Separating the cluster theory and models from Firecrown into an independent module.

examples/cluster_number_counts/cluster_redshift_richness.py

@@ -1,75 +1,73 @@
 """Likelihood factory function for cluster number counts."""
 
 import os
-from typing import Any, Dict
 
 import pyccl as ccl
 import sacc
 
-from firecrown.likelihood.gauss_family.statistic.cluster_number_counts import (
-    ClusterNumberCounts,
-)
+from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator
 from firecrown.likelihood.gauss_family.gaussian import ConstGaussian
+from firecrown.likelihood.gauss_family.statistic.binned_cluster_number_counts import (
+    BinnedClusterNumberCounts,
+)
+from firecrown.models.cluster.properties import ClusterProperty
 from firecrown.modeling_tools import ModelingTools
-from firecrown.models.cluster_abundance import ClusterAbundance
-from firecrown.models.cluster_mass_rich_proxy import ClusterMassRich
-from firecrown.models.cluster_redshift_spec import ClusterRedshiftSpec
-
+from firecrown.models.cluster.abundance import ClusterAbundance
+from firecrown.models.cluster.kernel import (
+    SpectroscopicRedshift,
+)
+from firecrown.models.cluster.mass_proxy import MurataBinned
+from firecrown.likelihood.likelihood import NamedParameters, Likelihood
+from typing import Tuple
 
-def build_likelihood(build_parameters):
-    """
-    Here we instantiate the number density (or mass function) object.
-    """
 
-    pivot_mass = 14.625862906
-    pivot_redshift = 0.6
+def get_cluster_abundance() -> ClusterAbundance:
+    hmf = ccl.halos.MassFuncBocquet16()
+    min_mass, max_mass = 13.0, 16.0
+    min_z, max_z = 0.2, 0.8
+    cluster_abundance = ClusterAbundance(min_mass, max_mass, min_z, max_z, hmf)
 
-    cluster_mass_r = ClusterMassRich(pivot_mass, pivot_redshift)
-    cluster_z = ClusterRedshiftSpec()
+    # Create and add the kernels you want in your cluster abundance
+    pivot_mass, pivot_redshift = 14.625862906, 0.6
+    mass_observable_kernel = MurataBinned(pivot_mass, pivot_redshift)
+    cluster_abundance.add_kernel(mass_observable_kernel)
 
-    # TODO: remove try/except when pyccl 3.0 is released
-    try:
-        hmd_200 = ccl.halos.MassDef200m()
-    except TypeError:
-        hmd_200 = ccl.halos.MassDef200m
+    redshift_proxy_kernel = SpectroscopicRedshift()
+    cluster_abundance.add_kernel(redshift_proxy_kernel)
 
-    hmf_args: Dict[str, Any] = {}
-    hmf_name = "Tinker08"
+    return cluster_abundance
 
-    cluster_abundance = ClusterAbundance(hmd_200, hmf_name, hmf_args)
 
-    stats = ClusterNumberCounts(
-        "numcosmo_simulated_redshift_richness",
-        cluster_abundance,
-        cluster_mass_r,
-        cluster_z,
-        use_cluster_counts=build_parameters.get_bool("use_cluster_counts", True),
-        use_mean_log_mass=build_parameters.get_bool("use_mean_log_mass", False),
+def build_likelihood(
+    build_parameters: NamedParameters,
+) -> Tuple[Likelihood, ModelingTools]:
+    """
+    Here we instantiate the number density (or mass function) object.
+    """
+    integrator = NumCosmoIntegrator()
+    # integrator = ScipyIntegrator()
+
+    # Pull params for the likelihood from build_parameters
+    average_properties = ClusterProperty.NONE
+    if build_parameters.get_bool("use_cluster_counts", True):
+        average_properties |= ClusterProperty.COUNTS
+    if build_parameters.get_bool("use_mean_log_mass", False):
+        average_properties |= ClusterProperty.MASS
+
+    survey_name = "numcosmo_simulated_redshift_richness"
+    likelihood = ConstGaussian(
+        [BinnedClusterNumberCounts(average_properties, survey_name, integrator)]
     )
-    stats_list = [stats]
 
-    # Here we instantiate the actual likelihood. The statistics argument carry
-    # the order of the data/theory vector.
-
-    lk = ConstGaussian(stats_list)
-
-    # We load the correct SACC file.
-    saccfile = os.path.expanduser(
-        os.path.expandvars(
-            "${FIRECROWN_DIR}/examples/cluster_number_counts/"
-            "cluster_redshift_richness_sacc_data.fits"
-        )
+    # Read in sacc data
+    sacc_file_nm = "cluster_redshift_richness_sacc_data.fits"
+    sacc_path = os.path.expanduser(
+        os.path.expandvars("${FIRECROWN_DIR}/examples/cluster_number_counts/")
     )
-    sacc_data = sacc.Sacc.load_fits(saccfile)
-
-    # The read likelihood method is called passing the loaded SACC file, the
-    # cluster number count functions will receive the appropriated sections
-    # of the SACC file.
-    lk.read(sacc_data)
+    sacc_data = sacc.Sacc.load_fits(os.path.join(sacc_path, sacc_file_nm))
+    likelihood.read(sacc_data)
 
-    # This script will be loaded by the appropriated connector. The framework
-    # then looks for the `likelihood` variable to find the instance that will
-    # be used to compute the likelihood.
-    modeling = ModelingTools()
+    cluster_abundance = get_cluster_abundance()
+    modeling_tools = ModelingTools(cluster_abundance=cluster_abundance)
 
-    return lk, modeling
+    return likelihood, modeling_tools

examples/cluster_number_counts/compare_integration_methods.py

@@ -1,26 +1,20 @@
 """Test integral methods for cluster abundance."""
-
-from typing import Any, Dict
 import time
 import itertools
 
-import pyccl as ccl
+import pyccl
 import numpy as np
-from numcosmo_py import Ncm
-
-from firecrown.models.cluster_abundance import ClusterAbundance
-from firecrown.models.cluster_mass_rich_proxy import (
-    ClusterMassRich,
-    ClusterMassRichBinArgument,
-)
-from firecrown.models.cluster_redshift_spec import (
-    ClusterRedshiftSpecArgument,
+from firecrown.models.cluster.mass_proxy import MurataBinned
+from firecrown.models.cluster.kernel import Kernel
+from firecrown.models.cluster.integrator.numcosmo_integrator import (
+    NumCosmoIntegrator,
+    NumCosmoIntegralMethod,
 )
+from firecrown.models.cluster.abundance import ClusterAbundance
+from firecrown.models.cluster.kernel import SpectroscopicRedshift
 
 
-def compare_integration():
-    """Compare integration methods."""
-
+def get_cosmology() -> pyccl.Cosmology:
     Omega_c = 0.262
     Omega_b = 0.049
     Omega_k = 0.0
@@ -32,7 +26,8 @@ def compare_integration():
     Neff = 3.046
     w0 = -1.0
     wa = 0.0
-    cosmo_ccl = ccl.Cosmology(
+
+    cosmo_ccl = pyccl.Cosmology(
         Omega_c=Omega_c,
         Omega_b=Omega_b,
         Neff=Neff,
@@ -45,56 +40,75 @@ def compare_integration():
         T_CMB=Tcmb0,
         m_nu=[0.00, 0.0, 0.0],
     )
+
+    return cosmo_ccl
+
+
+def get_mass_richness() -> Kernel:
     pivot_mass = 14.625862906
     pivot_redshift = 0.6
 
-    cluster_mass_r = ClusterMassRich(pivot_mass, pivot_redshift)
-
-    cluster_mass_r.mu_p0 = 3.0
-    cluster_mass_r.mu_p1 = 0.86
-    cluster_mass_r.mu_p2 = 0.0
-    cluster_mass_r.sigma_p0 = 3.0
-    cluster_mass_r.sigma_p1 = 0.7
-    cluster_mass_r.sigma_p2 = 0.0
-
-    # TODO: remove try/except when pyccl 3.0 is released
-    try:
-        hmd_200 = ccl.halos.MassDef200m()
-    except TypeError:
-        hmd_200 = ccl.halos.MassDef200m
-
-    hmf_args: Dict[str, Any] = {}
-    hmf_name = "Tinker08"
-    z_bins = np.linspace(0.0, 1.0, 4)
-    r_bins = np.linspace(1.0, 2.5, 5)
-
-    integ_options = [
-        Ncm.IntegralNDMethod.P,
-        Ncm.IntegralNDMethod.P_V,
-        Ncm.IntegralNDMethod.H,
-        Ncm.IntegralNDMethod.H_V,
-    ]
-    for integ_method in integ_options:
-        abundance_test = ClusterAbundance(
-            hmd_200, hmf_name, hmf_args, integ_method=integ_method
-        )
-        test_m_list = []
-        t1 = time.time()
+    mass_richness = MurataBinned(pivot_mass, pivot_redshift)
+
+    mass_richness.mu_p0 = 3.0
+    mass_richness.mu_p1 = 0.86
+    mass_richness.mu_p2 = 0.0
+    mass_richness.sigma_p0 = 3.0
+    mass_richness.sigma_p1 = 0.7
+    mass_richness.sigma_p2 = 0.0
+
+    return mass_richness
+
 
-        for i, j in itertools.product(range(3), range(4)):
-            cluster_mass_bin = ClusterMassRichBinArgument(
-                cluster_mass_r, 13, 15, r_bins[j], r_bins[j + 1]
+def compare_integration() -> None:
+    """Compare integration methods."""
+    hmf = pyccl.halos.MassFuncTinker08()
+    abundance = ClusterAbundance(13, 15, 0, 4, hmf)
+
+    mass_richness = get_mass_richness()
+    abundance.add_kernel(mass_richness)
+
+    redshift_proxy_kernel = SpectroscopicRedshift()
+    abundance.add_kernel(redshift_proxy_kernel)
+
+    cosmo = get_cosmology()
+    abundance.update_ingredients(cosmo)
+
+    sky_area = 360**2
+    integrand = abundance.get_integrand()
+
+    for method in NumCosmoIntegralMethod:
+        counts_list = []
+        t_start = time.time()
+
+        nc_integrator = NumCosmoIntegrator(method=method)
+        nc_integrator.set_integration_bounds(abundance, 496, (0, 4), (13, 15))
+
+        z_bins = np.linspace(0.0, 1.0, 4)
+        mass_proxy_bins = np.linspace(1.0, 2.5, 5)
+
+        for z_idx, mass_proxy_idx in itertools.product(range(3), range(4)):
+            z_proxy_limits = (z_bins[z_idx], z_bins[z_idx + 1])
+            mass_proxy_limits = (
+                mass_proxy_bins[mass_proxy_idx],
+                mass_proxy_bins[mass_proxy_idx + 1],
             )
-            cluster_z_bin = ClusterRedshiftSpecArgument(z_bins[i], z_bins[i + 1])
-            cluster_counts = abundance_test.compute(
-                cosmo_ccl, cluster_mass_bin, cluster_z_bin
+
+            nc_integrator.set_integration_bounds(
+                abundance,
+                sky_area,
+                z_proxy_limits,
+                mass_proxy_limits,
             )
-            test_m_list.append(cluster_counts)
-        t2 = time.time()
+
+            counts = nc_integrator.integrate(integrand)
+            counts_list.append(counts)
+
+        t_stop = time.time()
 
         print(
-            f"The time for {integ_method} is {t2-t1}\n\n"
-            f"The counts value is {test_m_list}\n\n"
+            f"The time for NumCosmo integration method {method} is {t_stop-t_start}\n\n"
+            f"The counts value is {counts_list}\n\n"
         )
 
 

examples/cluster_number_counts/generate_rich_mean_mass_sacc_data.py

@@ -19,11 +19,11 @@
 
 from astropy.io import fits
 from scipy import stats
-
+from typing import Any
 import sacc
 
 
-def generate_sacc_file():
+def generate_sacc_file() -> Any:
     """Generate a SACC file for cluster number counts."""
 
     H0 = 71.0
@@ -40,8 +40,7 @@ def generate_sacc_file():
     cosmo.add_submodel(prim)
 
     dist = Nc.Distance.new(2.0)
-
-    tf = Nc.TransferFunc.new_from_name("NcTransferFuncEH")
+    tf = Nc.TransferFuncEH.new()
 
     psml = Nc.PowspecMLTransfer.new(tf)
 
@@ -74,13 +73,9 @@ def generate_sacc_file():
     zu = 0.65
 
     # NumCosmo proxy model based on arxiv 1904.07524v2
-    cluster_z = Nc.ClusterRedshift.new_from_name(
-        f"NcClusterRedshiftNodist{{'z-min': <{zl:22.15e}>, 'z-max':<{zu:22.15e}>}}"
-    )
-
-    cluster_m = Nc.ClusterMass.new_from_name(
-        f"NcClusterMassAscaso{{'M0':<{3.0e14 / 0.71:22.15e}>,'z0':<0.6>, "
-        f"'lnRichness-min':<{lnRl:22.15e}>, 'lnRichness-max':<{lnRu:22.15e}>}}"
+    cluster_z = Nc.ClusterRedshiftNodist(z_max=zu, z_min=zl)
+    cluster_m = Nc.ClusterMassAscaso(
+        M0=3.0e14 / 0.71, z0=0.6, lnRichness_min=lnRl, lnRichness_max=lnRu
     )
     cluster_m.param_set_by_name("mup0", 3.19)
     cluster_m.param_set_by_name("mup1", 2 / np.log(10))

examples/des_y1_3x2pt/numcosmo_run.py

@@ -129,7 +129,7 @@ def setup_likelihood(
 def setup_fit(likelihood: Ncm.Likelihood, mset: Ncm.MSet) -> Ncm.Fit:
     """Setup the fit object."""
 
-    fit = Ncm.Fit.new(
+    fit = Ncm.Fit.factory(
         Ncm.FitType.NLOPT,
         "ln-neldermead",
         likelihood,
@@ -151,7 +151,7 @@ def setup_fit_all() -> Ncm.Fit:
     return fit
 
 
-def run_test():
+def run_test() -> None:
     """Run the fit."""
 
     fit = setup_fit_all()
@@ -162,7 +162,7 @@ def run_test():
     fit.run_restart(Ncm.FitRunMsgs.FULL, 1.0e-3, 0.0, None, None)
 
 
-def run_compute_best_fit():
+def run_compute_best_fit() -> None:
     """Run the fit."""
 
     fit = setup_fit_all()
@@ -178,7 +178,7 @@ def run_compute_best_fit():
     fit.run_restart(Ncm.FitRunMsgs.FULL, 1.0e-3, 0.0, None, None)
 
 
-def run_apes_sampler(ssize: int):
+def run_apes_sampler(ssize: int) -> None:
     """Run the fit."""
 
     fit = setup_fit_all()
@@ -192,7 +192,7 @@ def run_apes_sampler(ssize: int):
 
     nwalkers = mset.fparam_len * 100
     esmcmc = create_esmcmc(
-        fit.lh, mset, "des_y1_3x2pt_apes", nwalkers=nwalkers, nthreads=1
+        fit.props.likelihood, mset, "des_y1_3x2pt_apes", nwalkers=nwalkers, nthreads=1
     )
 
     esmcmc.start_run()

examples/des_y1_3x2pt/numcosmo_run_PT.py

@@ -77,7 +77,7 @@
 lh.priors_add_gauss_param_name(mset, "NcFirecrownPT:src0_delta_z", -0.001, 0.016)
 lh.priors_add_gauss_param_name(mset, "NcFirecrownPT:lens0_delta_z", +0.001, 0.008)
 
-fit = Ncm.Fit.new(
+fit = Ncm.Fit.factory(
     Ncm.FitType.NLOPT,
     "ln-neldermead",
     lh,

examples/srd_sn/numcosmo_run.py

@@ -81,7 +81,7 @@
 
 lh = Ncm.Likelihood(dataset=dset)
 
-fit = Ncm.Fit.new(
+fit = Ncm.Fit.factory(
     Ncm.FitType.NLOPT,
     "ln-neldermead",
     lh,