In [1]:
Copied!
# Run this notebook anywhere: install MeasureIA if it is not already available (e.g. on
# Google Colab). Everything else the notebook needs, including the mock catalogues and pyccl,
# comes with the package - no simulation data and no clone of the repository are required.
import importlib.util
import subprocess
import sys
if importlib.util.find_spec("measureia") is None:
subprocess.run([sys.executable, "-m", "pip", "install", "-q", "measureia"], check=True)
# Run this notebook anywhere: install MeasureIA if it is not already available (e.g. on
# Google Colab). Everything else the notebook needs, including the mock catalogues and pyccl,
# comes with the package - no simulation data and no clone of the repository are required.
import importlib.util
import subprocess
import sys
if importlib.util.find_spec("measureia") is None:
subprocess.run([sys.executable, "-m", "pip", "install", "-q", "measureia"], check=True)
In [2]:
Copied!
import matplotlib.pyplot as plt
import numpy as np
import pyccl as ccl
from measureia import MeasureIALightcone
from measureia import ReadData
from measureia.mocks import radial_alignment_lightcone_mock
import matplotlib.pyplot as plt
import numpy as np
import pyccl as ccl
from measureia import MeasureIALightcone
from measureia import ReadData
from measureia.mocks import radial_alignment_lightcone_mock
MeasureIA usage¶
MeasureIA can be used to measure intrinsic alignment correlation functions and clustering in lightcone data (observational or simulations). Both the w_gg, w_g+ and multipole moment estimators can be measured.
This notebook runs on a synthetic lightcone mock with a known radial-alignment signal, so it can be executed as-is without any survey data. Replace the data dictionaries below with your own catalogue to measure your own data.
Parameters for MeasureIA object¶
These are the parameters we need to enter to create the MeasureIA object. They govern measurement parameters, e.g. r(p) and pi bins (bounds, number of bins...).
In [3]:
Copied!
h = 0.7 # value of hubble parameter
num_bins_r = 8 # number of r or rp bins
num_bins_pi = 6 # number of pi bins. (trivial for wg+, not too little for multipoles)
separation_limits = [1.0, 20.0] # Mpc
pimax = 30. # pi max value
data_path_out = "./" # file path to output data
file_name = "example_IA_lightcone.hdf5" # datafile name (must be hdf5)
num_nodes = 1 # number of CPU cores available for calculation
IA_estimator = "galaxies" # type of estimator to be used. Choose "clusters" or "galaxies"
# Definition of the estimator will be printed when a method is called. "galaxies" uses the randoms (RR) in the
# denominator, "clusters" normalises by the data pairs instead; use the latter when the position sample is a
# cluster catalogue.
h = 0.7 # value of hubble parameter
num_bins_r = 8 # number of r or rp bins
num_bins_pi = 6 # number of pi bins. (trivial for wg+, not too little for multipoles)
separation_limits = [1.0, 20.0] # Mpc
pimax = 30. # pi max value
data_path_out = "./" # file path to output data
file_name = "example_IA_lightcone.hdf5" # datafile name (must be hdf5)
num_nodes = 1 # number of CPU cores available for calculation
IA_estimator = "galaxies" # type of estimator to be used. Choose "clusters" or "galaxies"
# Definition of the estimator will be printed when a method is called. "galaxies" uses the randoms (RR) in the
# denominator, "clusters" normalises by the data pairs instead; use the latter when the position sample is a
# cluster catalogue.
Parameters for methods¶
These are the parameters that are needed to configure the methods that are called. They govern what is measured and which methods are used. Check the method docstrings for extensive explanations.
In [4]:
Copied!
cosmology = ccl.Cosmology(Omega_c=0.27, Omega_b=0.049, h=h, sigma8=0.8, n_s=0.96) # pyccl cosmology to be used.
# If None, a default cosmology is used (Omega_c=0.225, Omega_b=0.045, sigma8=0.8, h=0.7, n_s=1.0).
over_h = False # if True, units are changed from Mpc -> Mpc/h
num_jk = 8 # number of jackknife patches. The covariance is measured whenever num_jk > 0 or
# jk_patches is given; use num_jk=0 to skip it.
corr_type = "both" # type of correlation to be calculated, choose g+, gg or both
masks = None # optional directory in form of data_dir containing masks to be placed over data in data_dir
masks_randoms = None # same as masks, but for the randoms
cosmology = ccl.Cosmology(Omega_c=0.27, Omega_b=0.049, h=h, sigma8=0.8, n_s=0.96) # pyccl cosmology to be used.
# If None, a default cosmology is used (Omega_c=0.225, Omega_b=0.045, sigma8=0.8, h=0.7, n_s=1.0).
over_h = False # if True, units are changed from Mpc -> Mpc/h
num_jk = 8 # number of jackknife patches. The covariance is measured whenever num_jk > 0 or
# jk_patches is given; use num_jk=0 to skip it.
corr_type = "both" # type of correlation to be calculated, choose g+, gg or both
masks = None # optional directory in form of data_dir containing masks to be placed over data in data_dir
masks_randoms = None # same as masks, but for the randoms
Create data dictionary¶
The data input is managed via a dictionary. These need to have specific keys (see below) and contain the position coordinates of both samples (RA, DEC, Redshift) and the shape information (e1, e2). The positional data is also needed for the random sample, which is stored in a separate dictionary.
Here both dictionaries come from the seeded mock catalogue in measureia.mocks: centrals uniform in comoving volume inside a cone section, with satellites scattered around them whose ellipticities point at their own central. The position (D) sample is the centrals and the shape (S) sample is the satellites. The mock is generated in comoving distance, so the distances are inverted to redshift with the same cosmology that is passed to the measurement methods.
In [5]:
Copied!
def redshift_of_comoving_distance(chi):
"""Invert the CCL comoving distance to get the redshift of each object."""
return 1. / ccl.scale_factor_of_chi(cosmology, chi) - 1.
# n_randoms_factor sets how many randoms per data object: enough randoms are needed to fill
# every rp/pi bin, or the estimator is undefined (NaN) there.
data, data_r, info = radial_alignment_lightcone_mock(n_centrals=400, n_sat=8, n_randoms_factor=10)
for catalogue in (data, data_r):
catalogue["Redshift"] = redshift_of_comoving_distance(catalogue.pop("r_com"))
catalogue["Redshift_shape_sample"] = redshift_of_comoving_distance(catalogue.pop("r_com_shape_sample"))
def redshift_of_comoving_distance(chi):
"""Invert the CCL comoving distance to get the redshift of each object."""
return 1. / ccl.scale_factor_of_chi(cosmology, chi) - 1.
# n_randoms_factor sets how many randoms per data object: enough randoms are needed to fill
# every rp/pi bin, or the estimator is undefined (NaN) there.
data, data_r, info = radial_alignment_lightcone_mock(n_centrals=400, n_sat=8, n_randoms_factor=10)
for catalogue in (data, data_r):
catalogue["Redshift"] = redshift_of_comoving_distance(catalogue.pop("r_com"))
catalogue["Redshift_shape_sample"] = redshift_of_comoving_distance(catalogue.pop("r_com_shape_sample"))
In [6]:
Copied!
# dictionary of randoms position data and shape-sample randoms (for the S+R term)
# data_r = {"Redshift": redshift of R_D sample,
# "Redshift_shape_sample": redshift of R_S sample (optional),
# "RA": RA of R_D sample, "RA_shape_sample": RA of R_S sample (optional),
# "DEC": DEC of R_D sample, "DEC_shape_sample": DEC of R_S sample (optional)}
# If only the Redshift, RA and DEC are provided, the code will use this random sample for both positions and shape clustering
# dictionary of the data (for the S+D term)
# data = {"Redshift": redshift of D sample, "Redshift_shape_sample": redshift of S sample,
# "RA": RA of D sample, "RA_shape_sample": RA of S sample,
# "DEC": DEC of D sample, "DEC_shape_sample": DEC of S sample,
# "e1": e1 of S sample, "e2": e2 of S sample}
print("data: ", {key: np.shape(value) for key, value in data.items()})
print("randoms:", {key: np.shape(value) for key, value in data_r.items()})
# dictionary of randoms position data and shape-sample randoms (for the S+R term)
# data_r = {"Redshift": redshift of R_D sample,
# "Redshift_shape_sample": redshift of R_S sample (optional),
# "RA": RA of R_D sample, "RA_shape_sample": RA of R_S sample (optional),
# "DEC": DEC of R_D sample, "DEC_shape_sample": DEC of R_S sample (optional)}
# If only the Redshift, RA and DEC are provided, the code will use this random sample for both positions and shape clustering
# dictionary of the data (for the S+D term)
# data = {"Redshift": redshift of D sample, "Redshift_shape_sample": redshift of S sample,
# "RA": RA of D sample, "RA_shape_sample": RA of S sample,
# "DEC": DEC of D sample, "DEC_shape_sample": DEC of S sample,
# "e1": e1 of S sample, "e2": e2 of S sample}
print("data: ", {key: np.shape(value) for key, value in data.items()})
print("randoms:", {key: np.shape(value) for key, value in data_r.items()})
data: {'RA': (400,), 'DEC': (400,), 'RA_shape_sample': (3200,), 'DEC_shape_sample': (3200,), 'e1': (3200,), 'e2': (3200,), 'weight': (400,), 'weight_shape_sample': (3200,), 'Redshift': (400,), 'Redshift_shape_sample': (3200,)}
randoms: {'RA': (4000,), 'DEC': (4000,), 'RA_shape_sample': (32000,), 'DEC_shape_sample': (32000,), 'weight': (4000,), 'weight_shape_sample': (32000,), 'Redshift': (4000,), 'Redshift_shape_sample': (32000,)}
Run with pre-assigned jackknife patches¶
If the jackknife patches are pre-assigned, they can be entered in the dictionary 'jk_patches' and entered in the input of the method.
Here they are generated with the assign_jackknife_patches method (kmeans on the sky), which is what the code does internally when num_jk is given instead; supply your own patch indices in the same format to control the regions yourself.
In [7]:
Copied!
MeasureIA_mock_obs = MeasureIALightcone(data, data_r,
separation_limits=separation_limits, pi_max=pimax,
num_bins_r=num_bins_r, num_nodes=num_nodes,
num_bins_pi=num_bins_pi, output_file_name=data_path_out + file_name)
dataset_name = "mock"
# patch index per object for the position (D), shape (S) and randoms (R_D, R_S) samples
patches = MeasureIA_mock_obs.assign_jackknife_patches(data, data_r, num_jk, seed=42)
# measure wgg, wg+
MeasureIA_mock_obs.measure_xi_w(IA_estimator, dataset_name, corr_type, jk_patches=patches,
masks=masks, masks_randoms=masks_randoms,
cosmology=cosmology, over_h=over_h)
# measure multipoles
MeasureIA_mock_obs.measure_xi_multipoles(IA_estimator, dataset_name, corr_type, jk_patches=patches,
masks=masks, masks_randoms=masks_randoms,
cosmology=cosmology, over_h=over_h)
MeasureIA_mock_obs = MeasureIALightcone(data, data_r,
separation_limits=separation_limits, pi_max=pimax,
num_bins_r=num_bins_r, num_nodes=num_nodes,
num_bins_pi=num_bins_pi, output_file_name=data_path_out + file_name)
dataset_name = "mock"
# patch index per object for the position (D), shape (S) and randoms (R_D, R_S) samples
patches = MeasureIA_mock_obs.assign_jackknife_patches(data, data_r, num_jk, seed=42)
# measure wgg, wg+
MeasureIA_mock_obs.measure_xi_w(IA_estimator, dataset_name, corr_type, jk_patches=patches,
masks=masks, masks_randoms=masks_randoms,
cosmology=cosmology, over_h=over_h)
# measure multipoles
MeasureIA_mock_obs.measure_xi_multipoles(IA_estimator, dataset_name, corr_type, jk_patches=patches,
masks=masks, masks_randoms=masks_randoms,
cosmology=cosmology, over_h=over_h)
MeasureIA object initialised with: observational data. There are 3200 galaxies in the shape sample and 400 galaxies in the position sample. The separation bin edges are given by [ 1. 1.45421543 2.11474253 3.07529122 4.47213595 6.50344913 9.45741609 13.75312044 20. ] Mpc. There are 8 r or r_p bins and 6 pi bins. The maximum pi used for binning is 30.0. The data will be written to ./example_IA_lightcone.hdf5 xi_g+ defined as (S+D - S+R)/RR, xi_gg as (SD - RD - SR)/RR + 1 There are 3200 galaxies in the shape sample and 400 galaxies in the position sample. There are 3200 galaxies in the shape sample and 4000 galaxies in the position sample.
There are 32000 galaxies in the shape sample and 400 galaxies in the position sample. There are 32000 galaxies in the shape sample and 4000 galaxies in the position sample.
{'D': 400, 'S': 3200, 'D_S': 0, 'R_D': 4000, 'R_S': 32000}
xi_g+ defined as (S+D - S+R)/RR, xi_gg as (SD - RD - SR)/RR + 1
There are 3200 galaxies in the shape sample and 400 galaxies in the position sample.
There are 3200 galaxies in the shape sample and 4000 galaxies in the position sample.
There are 32000 galaxies in the shape sample and 400 galaxies in the position sample. There are 32000 galaxies in the shape sample and 4000 galaxies in the position sample.
Run with internal jackknife patch assignment¶
If there are no pre-assigned jackknife patches, they can be assigned internally by entering the number of patches you want in num_jk. Note that for small samples, the covariance may differ because the patch assignment has a random seed.
In [8]:
Copied!
MeasureIA_mock_obs = MeasureIALightcone(data, data_r,
separation_limits=separation_limits, pi_max=pimax,
num_bins_r=num_bins_r, num_nodes=num_nodes,
num_bins_pi=num_bins_pi, output_file_name=data_path_out + file_name)
dataset_name_internal_jk = "mock_internal_jk"
# measure wgg, wg+
MeasureIA_mock_obs.measure_xi_w(IA_estimator, dataset_name_internal_jk, corr_type, num_jk=num_jk,
masks=masks, masks_randoms=masks_randoms,
cosmology=cosmology, over_h=over_h)
# measure multipoles
MeasureIA_mock_obs.measure_xi_multipoles(IA_estimator, dataset_name_internal_jk, corr_type, num_jk=num_jk,
masks=masks, masks_randoms=masks_randoms,
cosmology=cosmology, over_h=over_h)
MeasureIA_mock_obs = MeasureIALightcone(data, data_r,
separation_limits=separation_limits, pi_max=pimax,
num_bins_r=num_bins_r, num_nodes=num_nodes,
num_bins_pi=num_bins_pi, output_file_name=data_path_out + file_name)
dataset_name_internal_jk = "mock_internal_jk"
# measure wgg, wg+
MeasureIA_mock_obs.measure_xi_w(IA_estimator, dataset_name_internal_jk, corr_type, num_jk=num_jk,
masks=masks, masks_randoms=masks_randoms,
cosmology=cosmology, over_h=over_h)
# measure multipoles
MeasureIA_mock_obs.measure_xi_multipoles(IA_estimator, dataset_name_internal_jk, corr_type, num_jk=num_jk,
masks=masks, masks_randoms=masks_randoms,
cosmology=cosmology, over_h=over_h)
MeasureIA object initialised with: observational data. There are 3200 galaxies in the shape sample and 400 galaxies in the position sample. The separation bin edges are given by [ 1. 1.45421543 2.11474253 3.07529122 4.47213595 6.50344913 9.45741609 13.75312044 20. ] Mpc. There are 8 r or r_p bins and 6 pi bins. The maximum pi used for binning is 30.0. The data will be written to ./example_IA_lightcone.hdf5 xi_g+ defined as (S+D - S+R)/RR, xi_gg as (SD - RD - SR)/RR + 1 There are 3200 galaxies in the shape sample and 400 galaxies in the position sample. There are 3200 galaxies in the shape sample and 4000 galaxies in the position sample.
There are 32000 galaxies in the shape sample and 400 galaxies in the position sample. There are 32000 galaxies in the shape sample and 4000 galaxies in the position sample.
{'D': 400, 'S': 3200, 'D_S': 0, 'R_D': 4000, 'R_S': 32000}
xi_g+ defined as (S+D - S+R)/RR, xi_gg as (SD - RD - SR)/RR + 1
There are 3200 galaxies in the shape sample and 400 galaxies in the position sample.
There are 3200 galaxies in the shape sample and 4000 galaxies in the position sample.
There are 32000 galaxies in the shape sample and 400 galaxies in the position sample. There are 32000 galaxies in the shape sample and 4000 galaxies in the position sample.
Read data¶
The output file is read with the internal ReadData class, which fills the w_gg, w_g+ and multipole attributes
(and their jackknife errors) that are present in the file. See the box notebook for the alternative read_cat
route and the output file structure in the documentation.
The plots below mirror those of the box notebook, but the separations are in Mpc here because over_h = False.
Note that the largest-separation bins are consistent with zero for this small mock and can scatter negative;
such points simply fall off the logarithmic y-axis.
In [9]:
Copied!
IA_data = ReadData(None, file_name[:-len(".hdf5")], None, data_path=data_path_out)
IA_data.read_MeasureIA_output(dataset_name, num_jk)
IA_data = ReadData(None, file_name[:-len(".hdf5")], None, data_path=data_path_out)
IA_data.read_MeasureIA_output(dataset_name, num_jk)
Plot w_gg¶
In [10]:
Copied!
plt.figure()
plt.errorbar(IA_data.rp, IA_data.w_gg, yerr=IA_data.errors_w_gg,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$w_{gg}$")
plt.xlabel(r"$r_p [\mathrm{Mpc}]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.rp, IA_data.w_gg, yerr=IA_data.errors_w_gg,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$w_{gg}$")
plt.xlabel(r"$r_p [\mathrm{Mpc}]$")
plt.show()
Plot w_g+¶
In [11]:
Copied!
plt.figure()
plt.errorbar(IA_data.rp, IA_data.w_gp, yerr=IA_data.errors_w_gp,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$w_{g+}$")
plt.xlabel(r"$r_p [\mathrm{Mpc}]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.rp, IA_data.w_gp, yerr=IA_data.errors_w_gp,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$w_{g+}$")
plt.xlabel(r"$r_p [\mathrm{Mpc}]$")
plt.show()
Plot multipoles clustering¶
In [12]:
Copied!
plt.figure()
plt.errorbar(IA_data.r, IA_data.multipoles_gg, yerr=IA_data.errors_multipoles_gg,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$\tilde{\xi}_{gg,0}$")
plt.xlabel(r"$r [\mathrm{Mpc}]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.r, IA_data.multipoles_gg, yerr=IA_data.errors_multipoles_gg,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$\tilde{\xi}_{gg,0}$")
plt.xlabel(r"$r [\mathrm{Mpc}]$")
plt.show()
Plot multipoles g+¶
In [13]:
Copied!
plt.figure()
plt.errorbar(IA_data.r, IA_data.multipoles_gp, yerr=IA_data.errors_multipoles_gp,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$\tilde{\xi}_{g+,2}$")
plt.xlabel(r"$r [\mathrm{Mpc}]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.r, IA_data.multipoles_gp, yerr=IA_data.errors_multipoles_gp,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$\tilde{\xi}_{g+,2}$")
plt.xlabel(r"$r [\mathrm{Mpc}]$")
plt.show()
Plot r_p w_g+¶
Survey analyses often show the signal multiplied by r_p, which flattens the dynamic range.
In [14]:
Copied!
plt.figure()
plt.errorbar(IA_data.rp, IA_data.rp * IA_data.w_gp, yerr=IA_data.rp * IA_data.errors_w_gp,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$r_p w_{g+} [\mathrm{Mpc}^2]$")
plt.xlabel(r"$r_p [\mathrm{Mpc}]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.rp, IA_data.rp * IA_data.w_gp, yerr=IA_data.rp * IA_data.errors_w_gp,
marker="o", linestyle="none", capsize=3)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$r_p w_{g+} [\mathrm{Mpc}^2]$")
plt.xlabel(r"$r_p [\mathrm{Mpc}]$")
plt.show()
