Cacciato style CLF model

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

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

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

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)
+

docs/api/dsf.hod_cacciato.rst

@@ -0,0 +1,7 @@
+dsf.hod\_cacciato module
+========================
+
+.. automodule:: dsf.hod_cacciato
+   :members:
+   :show-inheritance:
+   :undoc-members:

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:

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
 ---------------

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

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: