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

examples/cluster_number_counts/cluster_redshift_richness.py

@@ -1,75 +1,78 @@
 """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.abundance_data import AbundanceData
 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(sky_area: float) -> 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, sky_area
+    )
 
-    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
+    use_cluster_counts = build_parameters.get_bool("use_cluster_counts", True)
+    use_mean_log_mass = build_parameters.get_bool("use_mean_log_mass", False)
+    survey_name = "numcosmo_simulated_redshift_richness"
+    likelihood = ConstGaussian(
+        [
+            BinnedClusterNumberCounts(
+                use_cluster_counts, use_mean_log_mass, 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()
+    sacc_adapter = AbundanceData(
+        sacc_data, survey_name, use_cluster_counts, use_mean_log_mass
+    )
+    cluster_abundance = get_cluster_abundance(sacc_adapter.survey_tracer.sky_area)
+    modeling_tools = ModelingTools(cluster_abundance=cluster_abundance)
 
-    return lk, modeling
+    return likelihood, modeling_tools

examples/cluster_number_counts/compare_integration_methods.py

@@ -8,7 +8,7 @@
 import numpy as np
 from numcosmo_py import Ncm
 
-from firecrown.models.cluster_abundance import ClusterAbundance
+from firecrown.models.cluster_abundance_old import ClusterAbundance
 from firecrown.models.cluster_mass_rich_proxy import (
     ClusterMassRich,
     ClusterMassRichBinArgument,

firecrown/connector/cosmosis/likelihood.py

@@ -7,8 +7,6 @@
 likelihood abstract base class; it the implementation of a CosmoSIS module,
 not a specific likelihood.
 """
-
-
 import cosmosis.datablock
 from cosmosis.datablock import option_section
 from cosmosis.datablock import names as section_names
@@ -119,7 +117,7 @@ def execute(self, sample: cosmosis.datablock) -> int:
             msg = self.form_error_message(exc)
             raise RuntimeError(msg) from exc
 
-        self.tools.prepare(ccl_cosmo)
+        self.tools.prepare(ccl_cosmo, firecrown_params)
         loglike = self.likelihood.compute_loglike(self.tools)
 
         derived_params_collection = self.likelihood.get_derived_parameters()

firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py

@@ -0,0 +1,149 @@
+"""This module holds classes needed to predict the binned cluster number counts
+
+The binned cluster number counts statistic predicts the number of galaxy
+clusters within a single redshift and mass bin.
+"""
+from __future__ import annotations
+from typing import List, Optional, Tuple
+import sacc
+import numpy as np
+from firecrown.models.cluster.integrator.integrator import Integrator
+from firecrown.models.cluster.abundance_data import AbundanceData
+from firecrown.likelihood.gauss_family.statistic.statistic import (
+    Statistic,
+    DataVector,
+    TheoryVector,
+)
+from firecrown.likelihood.gauss_family.statistic.source.source import SourceSystematic
+from firecrown.modeling_tools import ModelingTools
+
+
+class BinnedClusterNumberCounts(Statistic):
+    """The Binned Cluster Number Counts statistic
+
+    This class will make a prediction for the number of clusters in a z, mass bin
+    and compare that prediction to the data provided in the sacc file.
+    """
+
+    def __init__(
+        self,
+        cluster_counts: bool,
+        mean_log_mass: bool,
+        survey_name: str,
+        integrator: Integrator,
+        systematics: Optional[List[SourceSystematic]] = None,
+    ):
+        super().__init__()
+        self.systematics = systematics or []
+        self.theory_vector: Optional[TheoryVector] = None
+        self.use_cluster_counts = cluster_counts
+        self.use_mean_log_mass = mean_log_mass
+        self.survey_name = survey_name
+        self.integrator = integrator
+        self.data_vector = DataVector.from_list([])
+        self.sky_area = 0.0
+        self.bin_limits: List[List[Tuple[float, float]]] = []
+
+    def read(self, sacc_data: sacc.Sacc) -> None:
+        # Build the data vector and indices needed for the likelihood
+
+        data_vector = []
+        sacc_indices = []
+
+        sacc_types = sacc.data_types.standard_types
+        sacc_adapter = AbundanceData(
+            sacc_data, self.survey_name, self.use_cluster_counts, self.use_mean_log_mass
+        )
+
+        if self.use_cluster_counts:
+            # pylint: disable=no-member
+            data, indices = sacc_adapter.get_data_and_indices(sacc_types.cluster_counts)
+            data_vector += data
+            sacc_indices += indices
+
+        if self.use_mean_log_mass:
+            # pylint: disable=no-member
+            data, indices = sacc_adapter.get_data_and_indices(
+                sacc_types.cluster_mean_log_mass
+            )
+            data_vector += data
+            sacc_indices += indices
+
+        self.sky_area = sacc_adapter.survey_tracer.sky_area
+        # Note - this is the same for both cl mass and cl counts... Why do we need to
+        # specify a data type?
+
+        # pylint: disable=no-member
+        self.bin_limits = sacc_adapter.get_bin_limits(sacc_types.cluster_mean_log_mass)
+        self.data_vector = DataVector.from_list(data_vector)
+
+        self.sacc_indices = np.array(sacc_indices)
+        super().read(sacc_data)
+
+    def get_data_vector(self) -> DataVector:
+        assert self.data_vector is not None
+        return self.data_vector
+
+    def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector:
+        assert tools.cluster_abundance is not None
+
+        theory_vector_list: List[float] = []
+        cluster_counts = []
+
+        if not self.use_cluster_counts and not self.use_mean_log_mass:
+            return TheoryVector.from_list(theory_vector_list)
+
+        cluster_counts = self.get_binned_cluster_counts(tools)
+
+        if self.use_cluster_counts:
+            theory_vector_list += cluster_counts
+
+        if self.use_mean_log_mass:
+            theory_vector_list += self.get_binned_cluster_masses(tools, cluster_counts)
+
+        return TheoryVector.from_list(theory_vector_list)
+
+    def get_binned_cluster_masses(
+        self, tools: ModelingTools, cluster_counts: List[float]
+    ) -> List[float]:
+        """Computes the mean mass of clusters in each bin
+
+        Using the data from the sacc file, this function evaluates the likelihood for
+        a single point of the parameter space, and returns the predicted mean mass of
+        the clusters in each bin."""
+        assert tools.cluster_abundance is not None
+
+        cluster_masses = []
+        for (z_proxy_limits, mass_proxy_limits), counts in zip(
+            self.bin_limits, cluster_counts
+        ):
+            integrand = tools.cluster_abundance.get_integrand(avg_mass=True)
+            self.integrator.set_integration_bounds(
+                tools.cluster_abundance, z_proxy_limits, mass_proxy_limits
+            )
+
+            total_mass = self.integrator.integrate(integrand)
+            mean_mass = total_mass / counts
+            cluster_masses.append(mean_mass)
+
+        return cluster_masses
+
+    def get_binned_cluster_counts(self, tools: ModelingTools) -> List[float]:
+        """Computes the number of clusters in each bin
+
+        Using the data from the sacc file, this function evaluates the likelihood for
+        a single point of the parameter space, and returns the predicted number of
+        clusters in each bin."""
+        assert tools.cluster_abundance is not None
+
+        cluster_counts = []
+        for z_proxy_limits, mass_proxy_limits in self.bin_limits:
+            self.integrator.set_integration_bounds(
+                tools.cluster_abundance, z_proxy_limits, mass_proxy_limits
+            )
+
+            integrand = tools.cluster_abundance.get_integrand()
+            counts = self.integrator.integrate(integrand)
+            cluster_counts.append(counts)
+
+        return cluster_counts