diff --git a/benchmarks/data_vector/figures/cacciato_esd.png b/benchmarks/data_vector/figures/cacciato_esd.png new file mode 100644 index 0000000..34b0b5c Binary files /dev/null and b/benchmarks/data_vector/figures/cacciato_esd.png differ diff --git a/benchmarks/data_vector/figures/cacciato_esd_with_digitized_reference_data.png b/benchmarks/data_vector/figures/cacciato_esd_with_digitized_reference_data.png new file mode 100644 index 0000000..33f3220 Binary files /dev/null and b/benchmarks/data_vector/figures/cacciato_esd_with_digitized_reference_data.png differ diff --git a/benchmarks/data_vector/reference/L4_esd_data_points_mandelbaum2006_webplotdigitizer.csv b/benchmarks/data_vector/reference/L4_esd_data_points_mandelbaum2006_webplotdigitizer.csv new file mode 100644 index 0000000..2ab4de5 --- /dev/null +++ b/benchmarks/data_vector/reference/L4_esd_data_points_mandelbaum2006_webplotdigitizer.csv @@ -0,0 +1,17 @@ +# Column units: r [kpc/h], esd [h_Msun_per_pc2] +r,esd,esd_low,esd_high +43.34071252214734,46.234926819148335,41.17684665471254,51.523542437639506 +79.91631166493562,24.003720150613397,20.123757294882612,27.639752007858743 +142.9544428982114,12.000341424045063,10.111412265219643,13.922933184041968 +231.11919574541656,9.441688216431537,7.995683046415844,10.871894822978975 +300.6339806722028,7.411654344783589,5.8344433678054495,8.997908852209086 +368.0321831972245,6.057746328679473,4.768648974582442,7.280515258917068 +449.02384012074515,5.06469455274024,4.037446486260684,6.117753743856538 +547.839070065059,4.298911778583965,3.435626176191919,5.140695340845657 +634.3699555661033,2.953284253529595,1.8862532035915682,4.015641390528406 +819.6272676708162,3.112822652309281,2.5447574615582442,3.694334049527693 +1005.0697296830173,3.0043020104561236,2.5441897784245957,3.494415754778914 +1224.186907861586,2.358524111110068,1.9772921234000844,2.764094304138051 +1410.39624169179,1.997004166040032,1.5023597735869665,2.486244875781691 +1658.1328618924867,2.1533647042176414,1.8513417876161165,2.436218456803349 +1903.9196189473482,1.8095737720679073,1.4461843692718677,2.163912156869191 diff --git a/benchmarks/data_vector/reference/fitavgsig.hh.all.L4.lowfdev.csv b/benchmarks/data_vector/reference/fitavgsig.hh.all.L4.lowfdev.csv new file mode 100644 index 0000000..73d9147 --- /dev/null +++ b/benchmarks/data_vector/reference/fitavgsig.hh.all.L4.lowfdev.csv @@ -0,0 +1,49 @@ +# Column units: r [Mpc/h], esd [h_Msun_per_pc2], ignore the rest two columns +# link https://github.com/rmandelb/mandelbaum-data/tree/master/halomass-2006 +r esd x y +2.110000e-02 5.248745e+01 6.724000e+01 1.000000e+20 +2.340000e-02 4.904119e+01 5.230000e+01 1.000000e+20 +2.580000e-02 4.580982e+01 1.840000e+00 1.000000e+20 +2.850000e-02 4.256520e+01 1.191000e+01 1.000000e+20 +3.150000e-02 3.937092e+01 1.296000e+01 3.365000e+01 +3.490000e-02 3.621436e+01 5.672000e+01 2.753000e+01 +3.850000e-02 3.330866e+01 1.969000e+01 2.673000e+01 +4.260000e-02 3.045324e+01 6.270000e+01 2.296000e+01 +4.710000e-02 2.777154e+01 4.061000e+01 2.097000e+01 +5.200000e-02 2.529141e+01 9.190000e+00 1.953000e+01 +5.750000e-02 2.294024e+01 5.525000e+01 1.878000e+01 +6.350000e-02 2.078655e+01 4.716000e+01 1.464000e+01 +7.020000e-02 1.878100e+01 2.965000e+01 1.399000e+01 +7.759999e-02 1.694533e+01 3.580000e+01 1.261000e+01 +8.570000e-02 1.528360e+01 1.122000e+01 1.237000e+01 +9.480000e-02 1.365408e+01 1.601000e+01 1.197000e+01 +1.047000e-01 1.221829e+01 1.651000e+01 1.039000e+01 +1.157000e-01 1.094171e+01 2.893000e+01 8.880000e+00 +1.279000e-01 9.810948e+00 -1.155000e+01 8.270000e+00 +1.414000e-01 8.813763e+00 -9.400000e-01 7.080000e+00 +1.562000e-01 7.944407e+00 4.680000e+00 6.440000e+00 +1.726000e-01 7.177272e+00 -9.800000e-01 5.900000e+00 +1.908000e-01 6.501606e+00 4.120000e+00 4.940000e+00 +2.109000e-01 5.909760e+00 9.220000e+00 4.500000e+00 +2.330000e-01 5.392006e+00 9.600000e-01 4.180000e+00 +2.576000e-01 4.933710e+00 5.540000e+00 4.130000e+00 +2.846000e-01 4.532073e+00 7.410000e+00 3.340000e+00 +3.146000e-01 4.173744e+00 3.800000e-01 3.150000e+00 +3.477000e-01 3.854030e+00 2.140000e+00 2.980000e+00 +3.842000e-01 3.566707e+00 -9.400000e-01 2.570000e+00 +4.246000e-01 3.305684e+00 4.120000e+00 2.160000e+00 +4.693000e-01 3.067446e+00 1.850000e+00 2.150000e+00 +5.187000e-01 2.849532e+00 1.480000e+00 1.820000e+00 +5.732000e-01 2.649164e+00 7.000000e-02 1.630000e+00 +6.335000e-01 2.462961e+00 1.570000e+00 1.700000e+00 +7.001000e-01 2.289478e+00 -1.960000e+00 1.450000e+00 +7.738000e-01 2.126962e+00 2.400000e-01 1.400000e+00 +8.551000e-01 1.984150e+00 2.990000e+00 1.200000e+00 +9.450999e-01 1.824450e+00 3.450000e+00 1.060000e+00 +1.044500e+00 1.670020e+00 1.030000e+00 9.300000e-01 +1.154300e+00 1.524806e+00 7.600000e-01 8.800000e-01 +1.275700e+00 1.387311e+00 3.080000e+00 7.900000e-01 +1.409900e+00 1.257159e+00 1.810000e+00 7.700000e-01 +1.558100e+00 1.134268e+00 2.270000e+00 6.500000e-01 +1.722000e+00 1.019302e+00 1.720000e+00 6.700000e-01 +1.903100e+00 9.131610e-01 1.030000e+00 5.400000e-01 diff --git a/benchmarks/data_vector/reference/rebin.lum.all.L4.lowfdev.csv b/benchmarks/data_vector/reference/rebin.lum.all.L4.lowfdev.csv new file mode 100644 index 0000000..e902fa1 --- /dev/null +++ b/benchmarks/data_vector/reference/rebin.lum.all.L4.lowfdev.csv @@ -0,0 +1,16 @@ +# Column units: r [kpc/h], esd [h_Msun_per_pc2], error on esd [h_Msun_per_pc2] +r esd err +22.3381501378219 59.0072536889055 32.3226437603873 +30.6230495594924 25.2314166405795 16.5566830686141 +45.6553692342722 31.3126760678966 11.0571116653852 +64.2592829703091 42.2616435597132 8.96193106301768 +86.7539448022708 19.7353330502735 7.15987601403994 +124.152052387709 5.88900918789102 4.21342273036555 +185.19276347883 4.78830009300953 2.64786296318038 +276.25061103256 3.51524033131159 1.81263762604785 +412.152526695816 1.90112421376724 1.20453545107694 +614.848765123912 0.0549364224469053 0.822038248063202 +917.280771503268 1.97240967455718 0.556025750361684 +1219.70607161488 2.03844521029045 0.587874807711166 +1489.74599664043 2.06348482617828 0.497792177106917 +1819.57041886367 1.33977276073258 0.423112667723181 diff --git a/benchmarks/test_cacciato.py b/benchmarks/test_cacciato.py new file mode 100644 index 0000000..e698dd1 --- /dev/null +++ b/benchmarks/test_cacciato.py @@ -0,0 +1,171 @@ +import sys + +sys.path.insert(0, + "/hpc/group/cosmology/nlc38/from_NERSC/dp2_proj/sham_proj/CCLX/my_work/cacciato_hod" +) +import cmasher as cmr +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +import pyccl as ccl +from cacciato_inputs import cacciato_med_pars, cosmo_pars, magnitude_to_luminosity + +from dsf.data_vector.delta_sigma_builder import DeltaSigmaCalculator +from dsf.pk2d_cacciato_hod import pk2d_cacciato_hod + +mpl.rcParams["text.usetex"] = "True" + +benchmark_indir = "data_vector/reference" +benchmark_figdir = "data_vector/figures" + +if __name__=="__main__": + + magfaint, magbright = -20,-21 + log_L1 = np.log10(magnitude_to_luminosity(magfaint)) #log10( L1 / [Lsun/h^2] ) + log_L2 = np.log10(magnitude_to_luminosity(magbright)) #log10( L2 / [Lsun/h^2] ) + hodparams = dict(log_L1=log_L1, log_L2=log_L2, h=cosmo_pars['h'], **cacciato_med_pars) + + print(cosmo_pars) + cosmo = ccl.Cosmology(**cosmo_pars) + + # prepare inputs for ESD calculation + r = np.geomspace(0.1, 10.0, 20) + + z_lens = 0.1 + a_lens = 1.0 / (1.0 + z_lens) + + k_array = np.geomspace(1.0e-3, 30.0, 64) + a_array = np.linspace(0.3, 1.0, 16) + + def pk2d_func(*, cosmo, hodpars=hodparams): + return pk2d_cacciato_hod( + cosmo, + k_array=k_array, + a_array=a_array, + **hodpars + ) + + calculator = DeltaSigmaCalculator(pk2d_func=pk2d_func) + delta_sigma = calculator.delta_sigma( + r=r, + a=a_lens, + cosmo=cosmo, + ) + + print("predicted ESD signal:\nr, ESD") + print(np.c_[r,delta_sigma]) + + # Now prepare for ESD calculation of ESD by splitting the full lens + # redshift range in multiple tomographic bins + z = np.linspace(0.05, 0.15, 9) + n_z = np.exp(-0.5 * ((z - z_lens) / 0.05) ** 2) + + ## Make this chose the same cosmology + #cosmo = make_ccl_cosmology( + # transfer_function="eisenstein_hu", + # matter_power_spectrum="halofit", + #) + + delta_sigma_bin = calculator.delta_sigma_lens_bin( + r=r, + lens_dndz=(z, n_z), + cosmo=cosmo, + z_min=0.05, + z_max=0.15, + ) + # compare the binned calculation with that at one redshift + ratio = delta_sigma_bin / delta_sigma + + # ----------------------------------------------- + # Get the benchmarking data: Mandelbaum 2006 + import pandas as pd + # note the data units for ESD:hMsun/pc^2, r:kpc/h + ddf = pd.read_csv(f"{benchmark_indir}/L4_esd_data_points_mandelbaum2006_webplotdigitizer.csv", comment="#") + err = ddf["esd_high"] - ddf["esd_low"] # data errorbar + + # theory prediction from Mandelbaum+2005 + df = pd.read_csv(f"{benchmark_indir}/rebin.lum.all.L4.lowfdev.csv", sep=" ", comment="#") + # real data + tdf = pd.read_csv(f"{benchmark_indir}/fitavgsig.hh.all.L4.lowfdev.csv", sep=" ", comment="#") + # ----------------------------------------------- + + colors = cmr.take_cmap_colors( + "viridis", + 3, + cmap_range=(0.10, 0.90), + return_fmt="hex", + ) + + fig, (ax, ax_res) = plt.subplots( + 2, + 1, + figsize=(7.0, 5.2), + sharex=True, + gridspec_kw={"height_ratios": [3.0, 1.0], "hspace": 0.06}, + ) + + ax.plot( + r*cosmo_pars['h'], # Mpc/h + delta_sigma/cosmo_pars['h'], # hMsun/pc^2 + color=colors[0], + marker="o", + markersize=6, + label=rf"single redshift, $z_l={z_lens:.2f}$", + ) + + ax.errorbar( + ddf["r"]/1000, # Mpc/h + ddf["esd"], # hMsun/pc^2 + yerr=err, + label="Mandelbaum+2006\n(Webplotdigitizer)" + ) + + ax.errorbar( + df.r/1000, # Mpc/h + df.esd, # hMsun/pc^2 + yerr=df.err, + label="L4 data\n(Mandelbaum+2006)" + ) + + ax.plot( + tdf.r, # Mpc/h + tdf.esd, # hMsun/pc^2 + label="fitavgsig L4\n(Mandelbaum+2006)" + ) + + ax.scatter( + r*cosmo_pars['h'], # Mpc/h + delta_sigma_bin/cosmo_pars['h'], # hMsun/pc^2 + color=colors[2], + s=50, + label="lens-bin averaged", + zorder=3, + ) + + ax.set_xscale("log") + ax.set_yscale("log") + + ax.set_ylabel(r"$\Delta\Sigma(R)\ [hM_\odot / {\rm pc}^2]$", fontsize=15) + ax.tick_params(axis="both", which="major", labelsize=13) + ax.tick_params(axis="both", which="minor", labelsize=11) + ax.legend(fontsize=13, frameon=False, title=r"$M_r-5\log h \in [-21,-20]$") + + ax_res.axhline(1.0, color="lightgray", linewidth=1.4, linestyle="--") + ax_res.plot( + r, + ratio, + color=colors[1], + ) + + ax_res.set_xscale("log") + ax_res.set_ylabel("ratio", fontsize=15) + ax_res.set_xlabel(r"$R\ [{\rm Mpc}/h]$", fontsize=15) + ax_res.tick_params(axis="both", which="major", labelsize=13) + ax_res.tick_params(axis="both", which="minor", labelsize=11) + + fig.subplots_adjust(left=0.18, right=0.97, bottom=0.15, top=0.94) + + outfname = f"{benchmark_figdir}/cacciato_esd.png" + plt.savefig(outfname, bbox_inches="tight", dpi=200) + print("saved figure at ", outfname) + diff --git a/docs/api/dsf.hod_cacciato.rst b/docs/api/dsf.hod_cacciato.rst new file mode 100644 index 0000000..38647b7 --- /dev/null +++ b/docs/api/dsf.hod_cacciato.rst @@ -0,0 +1,7 @@ +dsf.hod\_cacciato module +======================== + +.. automodule:: dsf.hod_cacciato + :members: + :show-inheritance: + :undoc-members: diff --git a/docs/api/dsf.pk2d_cacciato_hod.rst b/docs/api/dsf.pk2d_cacciato_hod.rst new file mode 100644 index 0000000..39211ca --- /dev/null +++ b/docs/api/dsf.pk2d_cacciato_hod.rst @@ -0,0 +1,7 @@ +dsf.pk2d\_cacciato\_hod module +============================== + +.. automodule:: dsf.pk2d_cacciato_hod + :members: + :show-inheritance: + :undoc-members: diff --git a/docs/api/dsf.rst b/docs/api/dsf.rst index f256734..6669a0f 100644 --- a/docs/api/dsf.rst +++ b/docs/api/dsf.rst @@ -19,7 +19,9 @@ Submodules :maxdepth: 2 dsf.delta_sigma_forecast_builder + dsf.hod_cacciato dsf.modelling + dsf.pk2d_cacciato_hod Module contents --------------- diff --git a/docs/api/dsf.utils.rst b/docs/api/dsf.utils.rst index 1211f4b..0fe5061 100644 --- a/docs/api/dsf.utils.rst +++ b/docs/api/dsf.utils.rst @@ -9,6 +9,7 @@ Submodules dsf.utils.converters dsf.utils.integrators + dsf.utils.special_func dsf.utils.thread_limits dsf.utils.types dsf.utils.validators diff --git a/docs/api/dsf.utils.special_func.rst b/docs/api/dsf.utils.special_func.rst new file mode 100644 index 0000000..4bc91e7 --- /dev/null +++ b/docs/api/dsf.utils.special_func.rst @@ -0,0 +1,7 @@ +dsf.utils.special\_func module +============================== + +.. automodule:: dsf.utils.special_func + :members: + :undoc-members: + :show-inheritance: diff --git a/src/dsf/hod_cacciato.py b/src/dsf/hod_cacciato.py new file mode 100644 index 0000000..966182b --- /dev/null +++ b/src/dsf/hod_cacciato.py @@ -0,0 +1,342 @@ +import numpy as np +import pyccl as ccl +from numpy.typing import NDArray +from scipy.special import erf + +from dsf.utils.special_func import safe_upper_inc_gamma + +__all__ = [ + "CacciatoHOD" +] + +class CacciatoHOD(ccl.halos.HaloProfileHOD): + r""" + This class implements the galaxy-halo connection in the HOD framework based + on the Conditional Luminosity Function (CLF) model of Cacciato et al. 2013. + This framework is designed to model galaxy samples defined by luminosity + bins. + + Args: + mass_def: Halo mass definition (e.g., '200m' or a MassDef object). + + concentration: Concentration-mass relation (e.g., 'Duffy08' or a + Concentration object). To apply the `eta` multiplicative factor, + `concentration` should be a class wrapped around + `pyccl.halos.concentration` that manages the multiplication. + + log_L1/log_L2: :math:`\log_{10}` (luminosity bin edge) of the lens + sample being modelled. The edges are defined s.t. :math:`(\log_{10}L_2 + > \log_{10}L_1)` + + HOD model based on Cacciato+2013 has 9 free parameters as listed below: + + .. rubric:: The CLF Parameters + + * **log_L0**: + :math:`\log_{10}` of the normalization factor :math:`(L_0)` + of the central galaxy luminosity-halo mass scaling relation, :math:`L_c(M)`. + + * **log_M1**: + :math:`\log_{10}` of the characteristic halo mass scale :math:`(M_1)`, + such that :math:`L_c(M) \propto M^{\gamma_1}` for :math:`{\rm halo-mass}, + M \ll M_2`. + + * **gamma_1**: + slope of the central galaxy luminosity-halo mass scaling relation, + :math:`L_c(M)` at the low-mass end :math:`(M \ll M1)`. + + * **gamma_2**: + slope of the central galaxy luminosity-halo mass scaling relation, + :math:`L_c(M)` at the high-mass end :math:`(M \gg M1)`. + + * **sigma_c**: + scatter in the luminosities of central galaxies in the sample populating + a fixed halo mass, :math:`M`. + + * **alpha_s**: + faint end of the slope of the satellite galaxy occupation number. + Currently, assumed independent of (M,z), but can be generalised. + + * **p_cs**: + The coupling parameter between the characteristic central luminosity :math:`L_c(M)` and + the cut-off satellite luminosity :math:`L_s^{*}(M)` for a halo of mass M + + * **b_0, b_1, b_2**: + parameters affecting the overall scaling of the the satellite galaxy + HOD/CLF/LF. + + Additional nusance parameters in the CLF model: + + * **eta**: + A multiplicative factor, as described in Cacciato+2013, which modifies + the parent halo concentration as follows: + + .. math:: + + c_\mathrm{halo}(M) = (1+\eta) \times c_\mathrm{baseline}(M) + + Here, :math:`c_\mathrm{baseline}` is the chosen c-M relation + which gets marginalised-over by `eta`. `eta` acts via the c-M wrapper class + passed above to the `concentration` parameter. + + * **R_s**: + A multiplicative factor, alters the concentration of the satellite + distribution. R_s=1 corresponds to fiducial case, where satellites + density distribution follows the NFW profile of the parent halo. + + Note: + HODs are unitless. But the characteristic masses and luminosities hold + the information of units and scalings involved. Here, we require the + masses and luminosities in :math:`M_\odot/h` and :math:`L_\odot/h^2` respectively to make + sure that the Cacciato+13 CLF parameters can be directly passed. + + CCL HMF is given as dn/dlog10(m[:math:`M_\odot`]), so make sure to work with + masses in :math:`M_\odot` units out of the HOD definition and use logarithmic base + 10 mass bins for integration. + + To use the Cacciato+2013 HOD parameters as is, all the methods in this class expect + the mass to be in :math:`M_\odot/h`. So an internal conversion of pyCCL's :math:`M_\odot` + mass is implemented where necessary. + + For Cacciato specific usecase, we do not vary free parameters in + redshift dependent manner. + + """ + + def __init__( + self, + *, + # halo mass def + mass_def: str | ccl.halos.MassDef, + # mass profile + concentration: str | ccl.halos.Concentration, + # sample bin def + log_L1, log_L2, + # hubble const + h=0.739, + # CLF pars + log_L0=9.95, log_M1=11.24, + gamma_1=3.18, gamma_2=0.245, p_cs = 0.562, + sigma_c=0.157, alpha_s=-1.18, + b_0=-1.17, b_1=1.53, b_2=-0.217, + # central prof pars + eta=0., + # satellite prof pars + R_s=1., + ): + super().__init__(mass_def=mass_def, concentration=concentration) + # disable (set to None) the vanilla HOD parameters not applicable to this CLF. + self._disable_vanilla_hod_parameters() + + # define inputs that don't change throughout the analysis. + self.h = h + #self.cM = concentration #use if needed + self.ns_independent = False + # sample luminosity bin definition + self.log_L1 = log_L1 + self.log_L2 = log_L2 + # update/define the CLF free parameters + self.R_s = R_s + bg_0 = 1./self.R_s + # followings are to be fixed for Cacciato like work + bg_p=0. + bmax_0=1. + bmax_p=0. + self.update_parameters( + log_L0=log_L0,log_M1=log_M1, gamma_1=gamma_1, gamma_2=gamma_2, + sigma_c=sigma_c, alpha_s=alpha_s, b_0=b_0, b_1=b_1, b_2=b_2, + bg_0=bg_0, bg_p=bg_p, bmax_0=bmax_0, bmax_p=bmax_p, eta=eta, + p_cs=p_cs, + ) + + def _disable_vanilla_hod_parameters(self): + """Nullifies the standard HOD parameters not utilized by the Cacciato CLF framework.""" + unused_attrs = [ + "log10Mmin_0", "log10Mmin_p", "log10M0_0", "log10M0_p", + "log10M1_0", "log10M1_p", "siglnM_0", "siglnM_p", "alpha_0", "alpha_p" + ] + for attr in unused_attrs: + setattr(self, attr, None) + + @property + def log_L1(self): + """float: The faint-end luminosity bin boundary in log10(Lsun/h^2).""" + return self._log_L1 + + @log_L1.setter + def log_L1(self, log_L1): + """Updates the faint-end luminosity bin boundary and triggers updates.""" + self._log_L1 = log_L1 + + @property + def log_L2(self): + """float: The bright-end luminosity bin boundary in log10(Lsun/h^2).""" + return self._log_L2 + + @log_L2.setter + def log_L2(self, log_L2): + """Updates the bright-end luminosity bin boundary and triggers updates.""" + self._log_L2 = log_L2 + + @property + def h(self): + """float: The Hubble parameter value""" + return self._h + @h.setter + def h(self, h): + """Updates the value of the Hubble parameter thoughout the analysis.""" + self._h = h + # cache log10(h) whenever h is altered. + self._log10_h = np.log10(h) + + def update_parameters(self, **kwargs): + r"""Updates a selective set of CLF parameters or resets them to defaults. + + If an empty dictionary or no arguments are provided, the parameters will fall + back to the default CLF parameters of Cacciato+13. Passing a subset of + parameters will only update those specific values, keeping the remaining + parameters at their previously instantiated values. + + This update does not affect the sample luminosity bin edges and + hubble-constant definitions. + + Args: + kwargs: a dict of CLF free parameters (e.g., log_L0, log_M1, gamma_1) + + Note: + To match the Cacciato parametrization exactly, set: + * bg = bg_0 (i.e., bg_p = 0) + * bmax = 1 (i.e., bmax_0 = 1, bmax_p = 0) + + And redefine: + * 1 / b_g = R_s (where R_s is a multiplicative factor to the + fiducial c-M relation as defined in Cacciato, not a CCL + definition). + """ + + # Luminosity thresholds (calibrated as h-dependent: log10(L / [h^-2 Lsun])) + self.log_L0 = kwargs.get('log_L0', getattr(self, 'log_L0', 9.95)) + # Mass scales (calibrated as h-dependent: log10(M / [h^-1 Msun])) + self.log_M1 = kwargs.get('log_M1', getattr(self, 'log_M1', 11.24)) + # Dimensionless polynomial and shape parameters + self.gamma_1 = kwargs.get('gamma_1', getattr(self, 'gamma_1', 3.18)) + self.gamma_2 = kwargs.get('gamma_2', getattr(self, 'gamma_2', 0.245)) + self.sigma_c = kwargs.get('sigma_c', getattr(self, 'sigma_c', 0.157)) + self.alpha_s = kwargs.get('alpha_s', getattr(self, 'alpha_s', -1.18)) + #self.b_0 = kwargs.get('b_0', getattr(self, 'b_0', -1.17)) + self.b_0 = kwargs.get('b_0', getattr(self, 'b_0', -1.17)) + self.b_1 = kwargs.get('b_1', getattr(self, 'b_1', 1.53)) + self.b_2 = kwargs.get('b_2', getattr(self, 'b_2', -0.217)) + self.p_cs = kwargs.get('p_cs', getattr(self, 'p_cs', 0.562)) + # satellite concentration bias parameters + self.R_s = kwargs.get('R_s', getattr(self, 'R_s', 1)) + self.bg_0 = 1.0 / self.R_s + # central incompleteness modeling, as a constant value + self.fc_0 = kwargs.get('fc_0', getattr(self, 'fc_0', 1.)) + self.fc_p = kwargs.get('fc_p', getattr(self, 'fc_p', 0.)) + self.a_pivot = kwargs.get('a_pivot', getattr(self, 'a_pivot', 0.)) + # parent halo concentration scaling injected into concentration-Mass relation + self.concentration.eta = kwargs.get('eta', getattr(self, 'eta', 0.)) + + def _log_Lc(self, M_h: NDArray[np.float64]) -> NDArray[np.float64]: + r"""Calculates \log_{10} characteristic central luminosity as a + function of halo mass M_h + + Args: + M_h: Halo mass in units of Msun/h. + + Returns: + The $\log_{10}$ of the characteristic central luminosity in units of + Lsun/h^2. + """ + + return ( + self.log_L0 + + self.gamma_1 * (np.log10(M_h) - self.log_M1) + - (self.gamma_1 - self.gamma_2) * np.log10(1 + M_h / 10 ** self.log_M1) + ) + + def _log_L_star_s(self, M_h: NDArray[np.float64]) -> NDArray[np.float64]: + r"""Calculates \log_{10} of characteristic satellite luminosity as + a function of halo mass M_h + + This represents the cut-off luminosity of the satellite galaxies + populating a halo of mass M_h. + + Args: + M_h: Halo mass in Msun/h + + Returns: + The $\log_{10}$ of the characteristic satellite luminosity in units + of Lsun/h^2. + """ + return np.log10(self.p_cs) + self._log_Lc(M_h) + + def _log_phi_star_s(self, M_h: NDArray[np.float64]) -> NDArray[np.float64]: + r"""Calculates \log_{10} of the amplitude of satellite CLF as a function of halo + mass M_h + + Args: + M_h: Mass in Msun/h; Msun/h unit is required to extract the correct + response from this function using Cacciato+2013 parameters + + Returns: + \log_{10} of the amplitude of satellite CLF + """ + # the dimensionless Cacciato scaled-mass parameter: log10( M / [10^12 * Msun/h] ) + log_M12 = np.log10(M_h) - 12.0 + return self.b_0 + self.b_1 * log_M12 + self.b_2 * log_M12 ** 2 + + def _Nc(self, M: NDArray[np.float64], a=None) -> NDArray[np.float64]: + """Calculates the central HOD: The number of central galaxies with + luminosities between [L-dL, L+dL] populated in a halo of mass M. + + Args: + M: Mass in Msun; internally converted to Msun/h + + Note: + Because of the internal conversion of Msun to Msun/h units, + constraints on the mass parameters will still be in units of + Msun/h. + + Returns: + Central HOD for the given CLF parameters + """ + # internal conversion to Msun/h + M = M * self.h + + log_Lc = self._log_Lc(M) + denom = np.sqrt(2) * self.sigma_c + term_max = (self.log_L2 - log_Lc) / denom + term_min = (self.log_L1 - log_Lc) / denom + return 0.5 * (erf(term_max) - erf(term_min)) + + def _Ns(self, M: NDArray[np.float64], a=None) -> NDArray[np.float64]: + # Here, adding`a` as a kwarg, but are not needed for this work. Putting + # just to be able to run the delta_sigma_builder. Remove it in the + # future. + """Calculates the satellite HOD: The number of satellite galaxies with + luminosities between [L-dL, L+dL] populated in a halo of mass M. + + Args: + M: Halo mass in Msun; internally converted to Msun/h + + Note: + Because of the internal conversion of Msun to Msun/h units, + constraints on the mass parameters will still be in units of + Msun/h. + + Returns: + Satellite HOD for the given CLF parameters + """ + # internal conversion to Msun/h + M = M * self.h + + phi_star_s = 10 ** self._log_phi_star_s(M) + log_Ls = self._log_L_star_s(M) #Lsun/h^2 + shape = (self.alpha_s + 1) / 2.0 + x_min = 2.0 * (self.log_L1 - log_Ls) + x_max = 2.0 * (self.log_L2 - log_Ls) + integral_min = safe_upper_inc_gamma(shape, 10 ** x_min) + integral_max = safe_upper_inc_gamma(shape, 10 ** x_max) + return (phi_star_s / 2.0) * (integral_min - integral_max) diff --git a/src/dsf/pk2d_cacciato_hod.py b/src/dsf/pk2d_cacciato_hod.py new file mode 100644 index 0000000..8a4a585 --- /dev/null +++ b/src/dsf/pk2d_cacciato_hod.py @@ -0,0 +1,135 @@ +import numpy as np +import pyccl as ccl +from numpy.typing import NDArray + +from dsf.data_vector.profiles import density_weighted_power_spectrum +from dsf.hod_cacciato import CacciatoHOD +from dsf.modelling import ( + _validate_pk2d_grids, + make_ccl_cosmology, +) + +__all__ = [ + "MASS_DEF", + "CONCENTRATION", + "MASS_FUNCTION", + "HALO_BIAS", + "MATTER_PROFILE", + "HM_CALCULATOR", + "make_ccl_cosmology", + "pk2d_cacciato_hod", + "ScaledConcentration" +] + + +class ScaledConcentration(ccl.halos.Concentration): + """ A CCL concentration class wrapper that scales a baseline + concentration-Mass relation via a free parameter `eta` as (1+eta).""" + + name = 'Scaled_cM_Cacciato2013' + + def __init__(self, baseline_cM, *, eta=1.0, rho_type="matter", **baseline_cM_kwargs): + """ + baseline_cM: a pyccl.halos.concentration class + baseline_cM_kwargs: Arguments for a pyccl.halos.concentration class + eta: This multiplies the baseline c-M relation in the form (1+eta) + """ + self._eta = eta + self.mass_def = baseline_cM_kwargs.pop("mass_def") + + assert self.mass_def.rho_type == rho_type, "Inconsistent mass definition in concentration" + + print(f"Extra baseline concentration parameters: -- {baseline_cM_kwargs} --") + self.baseline_relation = baseline_cM(mass_def=self.mass_def, **baseline_cM_kwargs) + assert rho_type == self.baseline_relation.mass_def.rho_type + + super().__init__(mass_def=self.mass_def) + + def _check_mass_def_strict(self, mass_def): + """ + To exactly work with Cacciato+2013 conventions, mandate 200m + """ + return mass_def.name not in ["200m"] + + def _concentration(self, cosmo, M, a): + return (1.0+self.eta) * self.baseline_relation(cosmo, M, a) + + @property + def eta(self): + """float: The concentration-Mass relation modifier parameter that + scales it by an overall factor of (1+eta)""" + return self._eta + @eta.setter + def eta(self, eta): + """Updates the scaling factor (1+eta) of the c-M relation.""" + self._eta = eta + +MASS_DEF = ccl.halos.MassDef200m + +CONCENTRATION = ScaledConcentration( + ccl.halos.ConcentrationDuffy08, + mass_def=MASS_DEF, +) + +MASS_FUNCTION = ccl.halos.MassFuncTinker10( + mass_def=MASS_DEF, +) + +HALO_BIAS = ccl.halos.HaloBiasTinker10( + mass_def=MASS_DEF, +) + +MATTER_PROFILE = ccl.halos.HaloProfileNFW( + mass_def=MASS_DEF, + concentration=CONCENTRATION, + fourier_analytic=True, +) + +HM_CALCULATOR = ccl.halos.HMCalculator( + mass_function=MASS_FUNCTION, + halo_bias=HALO_BIAS, + mass_def=MASS_DEF, +) + +def pk2d_cacciato_hod( + cosmo: ccl.Cosmology, + *, + k_array: NDArray[np.float64], + a_array: NDArray[np.float64], + **hod_kwargs, +) -> ccl.Pk2D: + """Return the HOD galaxy-matter power spectrum for Delta Sigma. + + Args: + cosmo: CCL cosmology object. + k_array: Wavenumber grid used for the halo-model ``Pk2D`` spline. + a_array: Scale-factor grid used for the halo-model ``Pk2D`` spline. + **hod_kwargs: Keyword arguments passed to ``CacciatoHOD``. + + Returns: + Density-weighted galaxy-matter ``Pk2D`` object. + """ + k_arr, a_arr = _validate_pk2d_grids(k_array, a_array) + + galaxy_profile = CacciatoHOD( + mass_def=MASS_DEF, + concentration=CONCENTRATION, + **hod_kwargs, + ) + + pk_gm = ccl.halos.halomod_Pk2D( + cosmo, + HM_CALCULATOR, + galaxy_profile, + prof2=MATTER_PROFILE, + lk_arr=np.log(k_arr), + a_arr=a_arr, + p_of_k_a=cosmo.get_nonlin_power(), + ) + + def pk_gm_power(cosmo, k, a): + return pk_gm(k, a, cosmo=cosmo) + + return density_weighted_power_spectrum(cosmo, pk_gm_power) + + diff --git a/src/dsf/utils/special_func.py b/src/dsf/utils/special_func.py new file mode 100644 index 0000000..e3932ae --- /dev/null +++ b/src/dsf/utils/special_func.py @@ -0,0 +1,42 @@ +"""Special functions and their integrals""" + +from __future__ import annotations + +import numpy as np +from scipy.special import gamma, gammaincc + +__all__ = [ + "upper_inc_gamma", + "safe_upper_inc_gamma", +] + +def upper_inc_gamma(a, x): + r"""Computes the upper incomplete gamma function :math:`\Gamma(a, x)` + :math:`\forall a \in \mathbb{R}` and :math:`x \ge 0`. + + .. math:: + + \Gamma(a, x) = \int_x^\infty t^{a - 1}e^{-t} dt. + + This function uses the recurrence relation, + + .. math:: + + \Gamma(a, x) = \frac{\Gamma(a+1, x) - x^a \cdot e^{-x}}{a} + + to handle :math:`a \lt 0` cases, where scipy's implementation returns nan. + + Parameters: + a (float): The shape parameter, can be any real number. + x (float): The lower limit of upper-gamma function integral. + + Returns: + a float + """ + if x <= 0 and a <= 0: + return np.inf + if a > 0: + return gammaincc(a, x) * gamma(a) + return (upper_inc_gamma(a + 1, x) - (x ** a) * np.exp(-x)) / a + +safe_upper_inc_gamma = np.vectorize(upper_inc_gamma, excluded=["a"])