Generate “realistic” mock data for a cluster ensemble ===================================================== Generate cluster ensemble with “realistic” mass and redshift distribution ------------------------------------------------------------------------- In this notebook, we generate a ‘realistic’ mock cluster ensemble, i.e. given a halo mass function, and mass-concentration relation (more complexity, such as miscentering, could also be added but it’s not done here). The notebook is structured as follows: - Implementation of the CCL and NumCosmo mass functions. - Some comparisons between the different backends. - Implementation of a accept-reject method to generate cluster samples from the mass function. - Generation of a cluster catalog for mass and redshift bins followed by a comparison with the theoretical prediction Then, the remaining of the notebook repeats what was done in the ``Example_cluster_ensemble.ipynb``. - Geration of mock galaxy data for each cluster - Creating ClusterEnsemble object and estimation of individual excess surface density profiles - An analysis of the covariance of the stack between radial bins Draw pairs of mass and redshift according to a mass function (computed with either CCL or NumCosmo) --------------------------------------------------------------------------------------------------- .. raw:: html In this notebook we use the Tinker halo mass function, and demonstrate the implementation via CCL or NumCosmo. .. code:: ipython3 import sys import numpy as np import pyccl as ccl import matplotlib.pyplot as plt import scipy import clmm from clmm import GalaxyCluster, ClusterEnsemble, GCData from clmm import Cosmology from clmm.support import mock_data as mock %matplotlib inline .. code:: ipython3 # NumCosmo imports try: import gi gi.require_version("NumCosmo", "1.0") gi.require_version("NumCosmoMath", "1.0") except: pass from gi.repository import GObject from gi.repository import NumCosmo as Nc from gi.repository import NumCosmoMath as Ncm Ncm.cfg_init() Ncm.cfg_set_log_handler(lambda msg: sys.stdout.write(msg) and sys.stdout.flush()) For reproducibility: .. code:: ipython3 np.random.seed(11) Mass function - CCL implementation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 cosmo = ccl.Cosmology( Omega_c=0.26, Omega_b=0.04, h=0.7, sigma8=0.8, n_s=0.96, Neff=3.04, m_nu=1.0e-05, mass_split="single", ) hmd_200c = ccl.halos.MassDef(200, "critical") # For a different multiplicty function, the user must change this function below def tinker08_ccl(logm, z): mass = 10 ** (logm) hmf_200c = ccl.halos.MassFuncTinker08(mass_def=hmd_200c) nm = hmf_200c(cosmo, mass, 1.0 / (1 + z)) return nm # dn/dlog10M # Computing the volume element def dV_over_dOmega_dz(z): a = 1.0 / (1.0 + z) da = ccl.background.angular_diameter_distance(cosmo, a) E = ccl.background.h_over_h0(cosmo, a) return ((1.0 + z) ** 2) * (da**2) * ccl.physical_constants.CLIGHT_HMPC / cosmo["h"] / E def pdf_tinker08_ccl(logm, z): return tinker08_ccl(logm, z) * dV_over_dOmega_dz(z) Mass function - Numcosmo implementation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 Ncm.cfg_init() # cosmo_nc = Nc.HICosmoDEXcdm() cosmo_nc = Nc.HICosmo.new_from_name(Nc.HICosmo, "NcHICosmoDECpl{'massnu-length':<1>}") cosmo_nc.omega_x2omega_k() cosmo_nc.param_set_by_name("w0", -1.0) cosmo_nc.param_set_by_name("w1", 0.0) cosmo_nc.param_set_by_name("Tgamma0", 2.725) cosmo_nc.param_set_by_name("massnu_0", 0.0) cosmo_nc.param_set_by_name("H0", 70.0) cosmo_nc.param_set_by_name("Omegab", 0.04) cosmo_nc.param_set_by_name("Omegac", 0.26) cosmo_nc.param_set_by_name("Omegak", 0.0) # ENnu = 3.046 - 3.0 * \ # cosmo_nc.E2Press_mnu(1.0e10) / (cosmo_nc.E2Omega_g(1.0e10) # * (7.0/8.0*(4.0/11.0)**(4.0/3.0))) ENnu = 3.046 cosmo_nc.param_set_by_name("ENnu", ENnu) reion = Nc.HIReionCamb.new() prim = Nc.HIPrimPowerLaw.new() cosmo_nc.add_submodel(reion) cosmo_nc.add_submodel(prim) dist = Nc.Distance.new(2.0) dist.prepare_if_needed(cosmo_nc) tf = Nc.TransferFunc.new_from_name("NcTransferFuncEH") psml = Nc.PowspecMLTransfer.new(tf) psml.require_kmin(1.0e-6) psml.require_kmax(1.0e3) psf = Ncm.PowspecFilter.new(psml, Ncm.PowspecFilterType.TOPHAT) psf.set_best_lnr0() prim.props.n_SA = 0.96 old_amplitude = np.exp(prim.props.ln10e10ASA) prim.props.ln10e10ASA = np.log((0.8 / cosmo_nc.sigma8(psf)) ** 2 * old_amplitude) print(0.8, cosmo_nc.sigma8(psf)) mulf = Nc.MultiplicityFuncTinker.new() mulf.set_mdef(Nc.MultiplicityFuncMassDef.CRITICAL) mulf.set_Delta(200.0) mf = Nc.HaloMassFunction.new(dist, psf, mulf) # # New mass function object using the objects defined above. # def tinker08_nc(logm, z): lnm = logm * np.log(10.0) # convert log10(M) to ln(M) res = mf.dn_dlnM(cosmo_nc, lnm, z) return res * np.log(10.0) # convert dn/dlnM to dn/dlog10M def pdf_tinker08_nc(logm, z): return tinker08_nc(logm, z) * mf.dv_dzdomega(cosmo_nc, z) .. parsed-literal:: 0.8 0.8 Sanity Checks - compare the CCL and NC implementations ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Volume element ^^^^^^^^^^^^^^ .. code:: ipython3 z_arr = np.linspace(0.01, 2, 100) v_nc = [mf.dv_dzdomega(cosmo_nc, z) for z in z_arr] v_ccl = dV_over_dOmega_dz(z_arr) plt.plot(z_arr, np.abs(v_nc / v_ccl - 1) * 100) plt.xlabel("Redshift", size=14) plt.ylabel("Relative difference [%]", size=14) plt.grid(which="major", lw=0.5) plt.grid(which="minor", lw=0.2) .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_13_0.png Mass function ^^^^^^^^^^^^^ The CCL and Numcosmo implementations give different results for the Tinker08 HMF, especially at high redshift (>1.5). Need to understand why, although it is not a regime we should find ourselves in with cluster studies at optical wavelength. Might be that the CCL and NC cosmologies defined above are not actually the same (neutrinos?). .. code:: ipython3 logm_arr = np.linspace(14.0, 15) z_arr = np.linspace(0.01, 1.7, 2) t08_vec_ccl = np.vectorize(tinker08_ccl) t08_vec_nc = np.vectorize(tinker08_nc) for z in z_arr: res_ccl = t08_vec_ccl(logm_arr, z) res_nc = t08_vec_nc(logm_arr, z) plt.plot(logm_arr, res_ccl, label=f"z={z}, CCL") plt.plot(logm_arr, res_nc, label=f"z={z}, NC") plt.yscale("log") plt.xlabel("log10(M200,c)", size=14) plt.ylabel("dn/d$\log_{10}$M", size=14) plt.legend() plt.grid(which="major", lw=0.5) plt.grid(which="minor", lw=0.2) .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_15_0.png .. code:: ipython3 z_arr = np.linspace(0.01, 1.7, 100) logm = 15.0 res_ccl = t08_vec_ccl(logm, z_arr) res_nc = t08_vec_nc(logm, z_arr) plt.plot(z_arr, np.abs(res_ccl / res_nc - 1) * 100) plt.xlabel("Redshift", size=14) plt.ylabel("Relative difference [%]", size=14) plt.grid(which="major", lw=0.5) plt.grid(which="minor", lw=0.2) .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_16_0.png Acceptance-rejection method to sample (M,z) from the mass function ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 def bivariate_draw( pdf, N=1000, logm_min=14.0, logm_max=15.0, zmin=0.01, zmax=1, Ngrid_m=30, Ngrid_z=30 ): """ Uses the rejection method for generating random numbers derived from an arbitrary probability distribution. Parameters ---------- pdf : func 2d distribution function to sample N : int number of points to generate log_min,logm_max : float log10 mass range zmin,zmax : float redshift range Returns: -------- ran_logm : list accepted logm values ran_z : list accepted redshift values acceptance : float acceptance ratio of the method """ # maximum value of the pdf over the mass and redshift space. # the pdf is not monotonous in mass and redshift, so we need # to find the maximum numerically. Here we scan the space # with a regular grid and use the maximum value. # Accuracy of the results depends on Ngrid_m and Ngrid_z # This should probably be improved logM_arr = np.linspace(logm_min, logm_max, Ngrid_m) z_arr = np.logspace(np.log10(zmin), np.log10(zmax), Ngrid_z) p = [] for logM in logM_arr: for z in z_arr: p.append(pdf(logM, z)) pmax = np.max(p) pmin = np.min(p) # Counters naccept = 0 ntrial = 0 # Keeps drawing until N points are accepted ran_logm = [] # output list of random numbers ran_z = [] # output list of random numbers while naccept < N: # draw (logm,z) from uniform distribution # draw p from uniform distribution logm = np.random.uniform(logm_min, logm_max) # x' z = np.random.uniform(zmin, zmax) # x' p = np.random.uniform(0.0, pmax) # y' if p < pdf(logm, z): # keep the point # print(logm, z, p, pdf(logm, z)) ran_logm.append(logm) ran_z.append(z) naccept = naccept + 1 ntrial = ntrial + 1 ran_logm = np.asarray(ran_logm) ran_z = np.asarray(ran_z) acceptance = float(N / ntrial) print(f"acceptance = {acceptance}") return ran_logm, ran_z, acceptance Draw :math:`10^4` (M,z) pairs and check consistency with mass function ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 # Define the numbers of pairs to draw and the mass and redshift ranges N = 10000 logm_min = 14.0 logm_max = 14.5 zmin = 0.2 zmax = 1.0 # Choose between the CCL and NC implementation (NC is faster) # pdf = pdf_tinker08_ccl pdf = pdf_tinker08_nc .. code:: ipython3 # Normalisation factor for the pdf to be used later norm = ( 1.0 / scipy.integrate.dblquad( pdf, zmin, zmax, lambda x: logm_min, lambda x: logm_max, epsrel=1.0e-4 )[0] ) .. code:: ipython3 # Random draw ran_logm, ran_z, acceptance = bivariate_draw( pdf, N=N, logm_min=logm_min, logm_max=logm_max, zmin=zmin, zmax=zmax ) .. parsed-literal:: acceptance = 0.32901230506020923 .. code:: ipython3 plt.hist2d(ran_logm, ran_z, bins=10) pass .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_23_0.png Check distribution in mass bins (integrate over the full redshift range) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 hist_m = np.histogram(ran_logm, bins=50) .. code:: ipython3 print(zmin, zmax) predicted_T08_m_quad = [ norm * N * scipy.integrate.dblquad(pdf, zmin, zmax, lambda x: lmmin, lambda x: lmmax, epsrel=1.0e-4)[0] for lmmin, lmmax in zip(hist_m[1], hist_m[1][1:]) ] .. parsed-literal:: 0.2 1.0 .. code:: ipython3 bin_centers_m = 0.5 * (hist_m[1][:-1] + hist_m[1][1:]) plt.hist( ran_logm, bins=hist_m[1], color="teal", alpha=0.3, label="drawn masses - integrated over redshift range", ) plt.plot( bin_centers_m, np.array(predicted_T08_m_quad), "-o", color="green", label="Expected from theory" ) plt.yscale("log") plt.xlabel("$\log_{10}$M", size=14) plt.ylabel("Number", size=14) plt.title(f"{N:,} clusters drawn in log10M={logm_min, logm_max}, z={zmin, zmax}") plt.legend() plt.grid(which="major", lw=0.5) plt.grid(which="minor", lw=0.2) .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_27_0.png Check distribution in redshift bins (integrate over the full mass range) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 hist_z = np.histogram(ran_z, bins=50) print(logm_min, logm_max) predicted_T08_z_quad = [ norm * N * scipy.integrate.dblquad( pdf, zzmin, zzmax, lambda x: logm_min, lambda x: logm_max, epsrel=1.0e-4 )[0] for zzmin, zzmax in zip(hist_z[1], hist_z[1][1:]) ] .. parsed-literal:: 14.0 14.5 .. code:: ipython3 bin_centers_z = 0.5 * (hist_z[1][:-1] + hist_z[1][1:]) plt.hist( ran_z, bins=hist_z[1], color="orange", alpha=0.3, label="drawn redshift - integrated over mass range", ) plt.plot( bin_centers_z, np.array(predicted_T08_z_quad), "-o", color="darkorange", label="Expected from theory", ) plt.yscale("log") plt.xlabel("Redshift", size=14) plt.ylabel("Number", size=14) plt.title(f"{N:,} clusters drawn in log10M={logm_min, logm_max}, z={zmin, zmax}") plt.legend(loc=4) plt.grid(which="major", lw=0.5) plt.grid(which="minor", lw=0.2) .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_30_0.png Generating a cluster catalog and associated source catalogs ----------------------------------------------------------- - First, set the cluster masses and redshifts given the generated mock data. Only keep 30 clusters from N generated above to avoid memory issues. - Then, instantiate a CCL concentration object to compute the concentration for each cluster from a mass-concentration realtion (Duffy et al. 2008). The actual concentration for each cluster is drawn from a lognormal distribution around the mean value - Last, we randomly generate the cluster center coordinates over the full sky (from 0 to 360 deg for ra and from -90 to 90 deg to dec) Cluster ensemble masses, redshifts and positions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. code:: ipython3 # Randomly keep 30 (M,z) pairs from the 1000s generated abve. Nstack = 30 indices = np.random.choice(len(ran_logm), Nstack) cluster_m = 10 ** ran_logm[indices] cluster_z = ran_z[indices] # Concentration CCL object to compute the theoretical concentration conc_obj = ccl.halos.ConcentrationDuffy08(mass_def=hmd_200c) conc_list = [] for number in range(0, len(cluster_m)): a = 1.0 / (1.0 + (cluster_z[number])) # mean value of the concentration for that cluster lnc_mean = np.log(conc_obj(cosmo, M=(cluster_m[number]), a=a)) # random draw of actual concentration from normal distribution around lnc_mean, with a 0.14 scatter lnc = np.random.normal(lnc_mean, 0.14) conc_list.append(np.exp(lnc)) conc_list = np.array(conc_list) # randomly draw cluster positions over the full sky ra = np.random.random(N) * 360 # from 0 to 360 deg sindec = np.random.random(N) * 2 - 1 dec = np.arcsin(sindec) * 180 / np.pi # from -90 to 90 deg .. code:: ipython3 len(cluster_m) Nstack .. parsed-literal:: 30 Background galaxy catalog generation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For each cluster of the ensemble, we use ``mock_data`` to generate a background galaxy catalog and store the results in a ``GalaxyCluster`` object. Note that: - The cluster density profiles follow the NFW parametrisation - The source redshifts follow the Chang et al. distribution and have associated pdfs - The shapes include shape noise and shape measurement errors - Background galaxy catalogs are independent, even if the clusters are close (i.e., no common galaxy between two catalogs). - For each cluster we then compute - the tangential and cross :math:`\Delta\Sigma` for each background galaxy - the weights ``w_ls`` to be used to compute the corresponding radial profiles (see ``demo_compute_deltasigma_weights.ipynb`` notebook for details) The cluster objects are then stored in ``gclist``. .. code:: ipython3 import warnings warnings.simplefilter( "ignore" ) # just to prevent warning print out when looping over the cluster ensemble below .. code:: ipython3 %%time gclist = [] # number of galaxies in each cluster field (alternatively, can use the galaxy density instead) n_gals = 10000 # ngal_density = 10 cosmo_clmm = Cosmology(H0=71.0, Omega_dm0=0.265 - 0.0448, Omega_b0=0.0448, Omega_k0=0.0) # cosmo_clmm.set_be_cosmo(cosmo) gal_args = dict( cosmo=cosmo_clmm, zsrc="chang13", delta_so=200, massdef="critical", halo_profile_model="nfw", zsrc_max=3.0, field_size=10.0, shapenoise=0.04, photoz_sigma_unscaled=0.02, ngals=n_gals, mean_e_err=0.1, ) for i in range(Nstack): # generate background galaxy catalog for cluster i cl = clmm.GalaxyCluster( f"mock_cluster_{i:04}", ra[i], dec[i], cluster_z[i], galcat=mock.generate_galaxy_catalog( cluster_m[i], cluster_z[i], conc_list[i], cluster_ra=ra[i], cluster_dec=dec[i], zsrc_min=cluster_z[i] + 0.1, **gal_args, ), ) # compute DeltaSigma for each background galaxy cl.compute_tangential_and_cross_components( shape_component1="e1", shape_component2="e2", tan_component="DS_t", cross_component="DS_x", cosmo=cosmo_clmm, is_deltasigma=True, use_pdz=True, ) # compute the weights to be used to bluid the DeltaSigma radial profiles cl.compute_galaxy_weights( use_pdz=True, use_shape_noise=True, shape_component1="e1", shape_component2="e2", use_shape_error=True, shape_component1_err="e_err", shape_component2_err="e_err", weight_name="w_ls", cosmo=cosmo_clmm, is_deltasigma=True, add=True, ) # append the cluster in the list gclist.append(cl) .. parsed-literal:: CPU times: user 1.68 s, sys: 280 ms, total: 1.96 s Wall time: 1.96 s .. code:: ipython3 cl.galcat.columns .. parsed-literal:: Creating ClusterEnsemble object and estimation of individual excess surface density profiles -------------------------------------------------------------------------------------------- From the galaxy cluster object list ``gclist``, we instantiate a cluster ensemble object ``clusterensemble``. This instantiation step uses - the individual galaxy input :math:`\Delta\Sigma_+` and :math:`\Delta\Sigma_{\times}` values (computed in the previous step, ``DS_{t,x}``) - the corresponding individual input weights ``w_ls`` (computed in the previous step) to compute - the output tangential ``DS_t`` and cross signal ``DS_x`` binned profiles (where the binning is controlled by ``bins``) - the associated profile weights ``W_l`` (that will be used to compute the stacked profile) .. code:: ipython3 %%time ensemble_id = 1 bins = np.logspace(np.log10(0.3), np.log10(5), 10) clusterensemble = ClusterEnsemble( ensemble_id, gclist, tan_component_in="DS_t", cross_component_in="DS_x", tan_component_out="DS_t", cross_component_out="DS_x", weights_in="w_ls", weights_out="W_l", bins=bins, bin_units="Mpc", cosmo=cosmo_clmm, ) .. parsed-literal:: CPU times: user 280 ms, sys: 0 ns, total: 280 ms Wall time: 280 ms There is also the option to create an ``ClusterEnsemble`` object without the clusters list. Then, the user may compute the individual profile for each wanted cluster and compute the radial profile once all the indvidual profiles have been computed. This method may be reccomended if there a large number of clusters to avoid excess of memory allocation, since the user can generate each cluster separately, compute its individual profile and then delete the cluster object. .. code:: ipython3 ensemble_id2 = 2 empty_cluster_ensemble = ClusterEnsemble(ensemble_id2) for cluster in gclist: empty_cluster_ensemble.make_individual_radial_profile( galaxycluster=cluster, tan_component_in="DS_t", cross_component_in="DS_x", tan_component_out="DS_t", cross_component_out="DS_x", weights_in="w_ls", weights_out="W_l", bins=bins, bin_units="Mpc", cosmo=cosmo_clmm, ) A third option is where all clusters already have the profile computed, and we just add those to the ``ClusterEnsemble`` object: .. code:: ipython3 # add profiles to gclist for cluster in gclist: cluster.make_radial_profile( tan_component_in="DS_t", cross_component_in="DS_x", tan_component_out="DS_t", cross_component_out="DS_x", weights_in="w_ls", weights_out="W_l", bins=bins, bin_units="Mpc", cosmo=cosmo_clmm, ) .. code:: ipython3 ensemble_id3 = 3 empty_cluster_ensemble = ClusterEnsemble(ensemble_id3) for cluster in gclist: empty_cluster_ensemble.add_individual_radial_profile( galaxycluster=cluster, profile_table=cluster.profile, tan_component="DS_t", cross_component="DS_x", weights="W_l", ) Stacked profile of the cluster ensemble ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The stacked radial profile of the ensemble is then obtained as .. code:: ipython3 clusterensemble.make_stacked_radial_profile(tan_component="DS_t", cross_component="DS_x") Covariance (Bootstrap, sample, Jackknife) of the stack between radial bins ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Radial bins may be correlated and the ``ClusterEnsemble`` class provides three methods to compute the covariance matrix of the stacked signal, from the data: - The Sample covariance directly computes the covariance between radial bins of the ``N`` individual cluster profiles of the stack. - The Bootstrap approach is a resampling technique generating ``n_bootstrap`` ensembles of ``N`` randomly drawn clusters from the original ensemble, allowing for duplication. For each new ensemble, the stacked profile is computed and the covariance computed over the ``n_bootstrap`` stacks. - The Jackknife approach is another resampling technique, that divides the sky in a given number of regions :math:`N_{\rm region}` and computes the covariance removing one region (i.e the clusters of the ensemble in that region) at a time. The stack is then computed using the remaining clusters and the covariance computed over the :math:`N_{\rm region}` number of stacks. The division of the sky is done using the Healpix pixelisation scheme and is controlled by the ``n_side`` parameter, with :math:`N_{\rm region}=12 N_{\rm side}^2`. NB: Approaches exist to compute the theoretical covariance of a stack but these are not (yet) available in ``CLMM`` .. code:: ipython3 clusterensemble.compute_sample_covariance(tan_component="DS_t", cross_component="DS_x") clusterensemble.compute_bootstrap_covariance( tan_component="DS_t", cross_component="DS_x", n_bootstrap=300 ) clusterensemble.compute_jackknife_covariance( n_side=16, tan_component="DS_t", cross_component="DS_x" ) .. code:: ipython3 plt.figure(figsize=(7, 7)) plt.rcParams["axes.linewidth"] = 2 plt.plot( clusterensemble.data["radius"][0], clusterensemble.cov["tan_sc"].diagonal() ** 0.5 / 1e13, "--", c="royalblue", label="Sample", linewidth=3, ) plt.plot( clusterensemble.data["radius"][0], clusterensemble.cov["tan_jk"].diagonal() ** 0.5 / 1e13, "-s", c="g", label="Bootstrap", linewidth=3, markersize=10, ) plt.plot( clusterensemble.data["radius"][0], clusterensemble.cov["tan_sc"].diagonal() ** 0.5 / 1e13, c="r", label="Jackknife", linewidth=3, ) plt.xlabel("R [Mpc]", fontsize=20) plt.ylabel(r"$\sigma_{\Delta\Sigma}\ (\times 10^{13} M_\odot /Mpc^2)$", fontsize=25) plt.tick_params(axis="both", which="major", labelsize=18) plt.legend(frameon=False, fontsize=20) plt.loglog() plt.grid(which="major", lw=0.5) plt.grid(which="minor", lw=0.2) .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_50_0.png .. code:: ipython3 fig, axes = plt.subplots(1, 3, figsize=(20, 5)) plt.rcParams["axes.linewidth"] = 2 fig.subplots_adjust(wspace=0.15, hspace=0) maximum = clusterensemble.cov["tan_sc"].max() for ax, cov, label in zip( axes, [clusterensemble.cov["tan_sc"], clusterensemble.cov["tan_jk"], clusterensemble.cov["tan_sc"]], ["Stack : Sample", "Stack : Bootstrap", "Stack : Jackknife"], ): ax.set_title(label, fontsize=20) ax.set_xlabel("radial bin index", fontsize=18) ax.set_ylabel("radial bin index", fontsize=18) ax.tick_params(axis="both", which="major", labelsize=12) im = ax.imshow(cov, cmap="Reds", vmin=0, vmax=maximum, origin="lower") plt.colorbar(im, ax=ax) .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_51_0.png Visualizing the stacked profiles and corresponding model ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ In the figure below, we plot: - the individual :math:`\Delta\Sigma` profiles of the clusters (light blue) - the stacked signal (red symbols) - the prediction computed using a NFW profile and the mean values of the mass, concentration and redshift in the stack (dashed black) .. code:: ipython3 moo = clmm.Modeling(massdef="critical", delta_mdef=200, halo_profile_model="nfw") moo.set_cosmo(cosmo_clmm) # Average values of mass and concentration of the ensemble to be used below # to overplot the model on the stacked profile moo.set_concentration(conc_list.mean()) moo.set_mass(cluster_m.mean()) .. code:: ipython3 r_stack, gt_stack, gx_stack = (clusterensemble.stacked_data[c] for c in ("radius", "DS_t", "DS_x")) plt.rcParams["axes.linewidth"] = 2 fig, axs = plt.subplots(1, 2, figsize=(17, 6)) err_gt = clusterensemble.cov["tan_sc"].diagonal() ** 0.5 / 1e13 err_gx = clusterensemble.cov["cross_sc"].diagonal() ** 0.5 / 1e13 axs[0].errorbar( r_stack, gt_stack / 1e13, err_gt, markersize=5, c="r", fmt="o", capsize=10, elinewidth=1, zorder=1000, alpha=1, label="Stack", ) axs[1].errorbar( r_stack, gx_stack / 1e13, err_gx, markersize=5, c="r", fmt="o", capsize=10, elinewidth=1, zorder=1000, alpha=1, label="Stack", ) # print(moo.eval_excess_surface_density(clusterensemble.data['radius'][0], cluster_z.mean())/1e13) axs[0].plot( clusterensemble.data["radius"][0], moo.eval_excess_surface_density(clusterensemble.data["radius"][0], cluster_z.mean()) / 1e13, "--k", linewidth=3, label="Prediction from stack mean cluster", zorder=100, ) axs[1].plot( clusterensemble.data["radius"][0], 0 * moo.eval_excess_surface_density(clusterensemble.data["radius"][0], cluster_z.mean()) / 1e13, "--k", linewidth=3, label="y=0", zorder=100, ) axs[0].set_xscale("log") axs[1].set_xscale("log") for i in range(Nstack): axs[0].plot( clusterensemble.data["radius"][i], clusterensemble.data["DS_t"][i] / 1e13, color="cyan", label="Individual", alpha=1, linewidth=1, ) axs[1].plot( clusterensemble.data["radius"][i], clusterensemble.data["DS_x"][i] / 1e13, color="cyan", label="Individual", alpha=1, linewidth=1, ) if i == 0: axs[0].legend(frameon=False, fontsize=15) axs[1].legend(frameon=False, fontsize=15) # axs[0].plot(np.average(clusterensemble.data['radius'], axis = 0), np.average(clusterensemble.data['gt'], weights = None, axis = 0)/1e13) axs[0].set_xlabel("R [Mpc]", fontsize=20) axs[1].set_xlabel("R [Mpc]", fontsize=20) axs[0].tick_params(axis="both", which="major", labelsize=18) axs[1].tick_params(axis="both", which="major", labelsize=18) axs[0].set_ylabel(r"$\Delta\Sigma_+$ $[\times 10^{13} M_\odot /Mpc^])$", fontsize=20) axs[1].set_ylabel(r"$\Delta\Sigma_\times$ $[\times 10^{13} M_\odot /Mpc^2]$", fontsize=20) axs[0].set_title(r"Tangential", fontsize=20) axs[1].set_title(r"Cross", fontsize=20) for ax in axs: ax.minorticks_on() ax.grid(lw=0.5) ax.grid(which="minor", lw=0.1) plt.show() .. image:: demo_mock_ensemble_realistic_files/demo_mock_ensemble_realistic_54_0.png Saving/Loading ClusterEnsemble ------------------------------ The ``ClusterEnsemble`` object also have an option for saving/loading: .. code:: ipython3 clusterensemble.save("ce.pkl") .. code:: ipython3 clusterensemble2 = ClusterEnsemble.load("ce.pkl")