From 0752c3cf0be83e166a94d86d3a381aae6e229ef3 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 13:40:02 -0800 Subject: [PATCH 01/26] Removing "generic" logic from integrators. Keeping integrators simple, they accept a function to integrate and perform that integration. Adding the ability to pass extra arguments. --- .../models/cluster/integrator/integrator.py | 48 +-------- .../cluster/integrator/numcosmo_integrator.py | 101 +++--------------- .../cluster/integrator/scipy_integrator.py | 89 +-------------- 3 files changed, 19 insertions(+), 219 deletions(-) diff --git a/firecrown/models/cluster/integrator/integrator.py b/firecrown/models/cluster/integrator/integrator.py index 4503753b..59caa4b1 100644 --- a/firecrown/models/cluster/integrator/integrator.py +++ b/firecrown/models/cluster/integrator/integrator.py @@ -3,11 +3,8 @@ This module holds the classes that define the interface required to integrate an assembled cluster abundance. """ -import inspect from abc import ABC, abstractmethod -from typing import Tuple, Dict, get_args, List -from firecrown.models.cluster.abundance import ClusterAbundance, AbundanceIntegrand -from firecrown.models.cluster.kernel import KernelType +from typing import Tuple, List, Callable class Integrator(ABC): @@ -17,48 +14,9 @@ class Integrator(ABC): a specific set of methods to be used to integrate a cluster abundance integral.""" def __init__(self) -> None: - self.z_proxy_limits: Tuple[float, float] = (-1.0, -1.0) - self.mass_proxy_limits: Tuple[float, float] = (-1.0, -1.0) - self.sky_area: float = 360**2 self.integral_bounds: List[Tuple[float, float]] = [] - self.integral_args_lkp: Dict[KernelType, int] = self._default_integral_args() + self.extra_args: List[float] = [] @abstractmethod - def integrate( - self, - integrand: AbundanceIntegrand, - ) -> float: + def integrate(self, func_to_integrate: Callable) -> float: """Call this method to integrate the provided integrand argument.""" - - @abstractmethod - def set_integration_bounds( - self, - cl_abundance: ClusterAbundance, - sky_area: float, - z_proxy_limits: Tuple[float, float], - mass_proxy_limits: Tuple[float, float], - ) -> None: - """Sets the limits of integration used by the integration library. - - This method uses the mass and redshift proxy arguments, along with - the cluster abundance argument to construct the limits of integration - to be used in evaluating the cluster abundance.""" - - def _default_integral_args(self) -> Dict[KernelType, int]: - lkp: Dict[KernelType, int] = {} - lkp[KernelType.MASS] = 0 - lkp[KernelType.Z] = 1 - return lkp - - def _validate_integrand(self, integrand: AbundanceIntegrand) -> None: - expected_args, expected_return = get_args(AbundanceIntegrand) - - signature = inspect.signature(integrand) - params = signature.parameters.values() - param_types = [param.annotation for param in params] - - assert len(params) == len(expected_args) - - assert param_types == list(expected_args) - - assert signature.return_annotation == expected_return diff --git a/firecrown/models/cluster/integrator/numcosmo_integrator.py b/firecrown/models/cluster/integrator/numcosmo_integrator.py index d3105307..487a1935 100644 --- a/firecrown/models/cluster/integrator/numcosmo_integrator.py +++ b/firecrown/models/cluster/integrator/numcosmo_integrator.py @@ -2,14 +2,11 @@ This module holds the NumCosmo implementation of the integrator classes """ -from typing import Tuple, Callable, Sequence, Optional +from typing import Tuple, Callable, Sequence, Optional, List from enum import Enum import numpy as np import numpy.typing as npt from numcosmo_py import Ncm -from firecrown.models.cluster.kernel import KernelType -from firecrown.models.cluster.abundance import ClusterAbundance, AbundanceIntegrand -from firecrown.models.cluster.integrator.integrator import Integrator class NumCosmoIntegralMethod(Enum): @@ -21,7 +18,7 @@ class NumCosmoIntegralMethod(Enum): H_V = Ncm.IntegralNDMethod.H_V -class NumCosmoIntegrator(Integrator): +class NumCosmoIntegrator: """The NumCosmo implementation of the Integrator base class.""" def __init__( @@ -31,82 +28,18 @@ def __init__( absolute_tolerance: float = 1e-12, ) -> None: super().__init__() + self.integral_bounds: List[Tuple[float, float]] = [] + self.extra_args: List[float] = [] self._relative_tolerance = relative_tolerance self._absolute_tolerance = absolute_tolerance self._method = method or NumCosmoIntegralMethod.P_V - def _integral_wrapper( - self, - integrand: AbundanceIntegrand, - ) -> Callable[[npt.NDArray], Sequence[float]]: - self._validate_integrand(integrand) - - # mypy strict issue: npt.NDArray[npt.NDArray[np.float64]] not supported - def ncm_integrand(int_args: npt.NDArray) -> Sequence[float]: - default = np.ones_like(int_args[0]) * -1.0 - # pylint: disable=R0801 - mass = self._get_or_default(int_args, KernelType.MASS, default) - z = self._get_or_default(int_args, KernelType.Z, default) - mass_proxy = self._get_or_default(int_args, KernelType.MASS_PROXY, default) - z_proxy = self._get_or_default(int_args, KernelType.Z_PROXY, default) - - return_val = integrand( - mass, - z, - self.sky_area, - mass_proxy, - z_proxy, - self.mass_proxy_limits, - self.z_proxy_limits, - ).tolist() - assert isinstance(return_val, list) - return return_val - - return ncm_integrand - - def set_integration_bounds( - self, - cl_abundance: ClusterAbundance, - sky_area: float, - z_proxy_limits: Tuple[float, float], - mass_proxy_limits: Tuple[float, float], - ) -> None: - # pylint: disable=R0801 - self.integral_args_lkp = self._default_integral_args() - self.integral_bounds = [ - (cl_abundance.min_mass, cl_abundance.max_mass), - (cl_abundance.min_z, cl_abundance.max_z), - ] - - self.mass_proxy_limits = mass_proxy_limits - self.z_proxy_limits = z_proxy_limits - self.sky_area = sky_area - - for kernel in cl_abundance.dirac_delta_kernels: - if kernel.kernel_type == KernelType.Z_PROXY: - self.integral_bounds[1] = z_proxy_limits - - elif kernel.kernel_type == KernelType.MASS_PROXY: - self.integral_bounds[0] = mass_proxy_limits - - for kernel in cl_abundance.integrable_kernels: - idx = len(self.integral_bounds) - - if kernel.kernel_type == KernelType.Z_PROXY: - self.integral_bounds.append(z_proxy_limits) - self.integral_args_lkp[KernelType.Z_PROXY] = idx - - elif kernel.kernel_type == KernelType.MASS_PROXY: - self.integral_bounds.append(mass_proxy_limits) - self.integral_args_lkp[KernelType.MASS_PROXY] = idx - - def integrate( - self, - integrand: AbundanceIntegrand, - ) -> float: + def integrate(self, func_to_integrate: Callable) -> float: Ncm.cfg_init() - ncm_integrand = self._integral_wrapper(integrand) - int_nd = CountsIntegralND(len(self.integral_bounds), ncm_integrand) + + int_nd = CountsIntegralND( + len(self.integral_bounds), func_to_integrate, self.extra_args + ) int_nd.set_method(self._method.value) int_nd.set_reltol(self._relative_tolerance) int_nd.set_abstol(self._absolute_tolerance) @@ -117,18 +50,6 @@ def integrate( int_nd.eval(Ncm.Vector.new_array(bl), Ncm.Vector.new_array(bu), res, err) return res.get(0) - def _get_or_default( - self, - # mypy strict issue: npt.NDArray[npt.NDArray[np.float64]] not supported - int_args: npt.NDArray, - kernel_type: KernelType, - default: npt.NDArray[np.float64], - ) -> npt.NDArray[np.float64]: - try: - return int_args[:, self.integral_args_lkp[kernel_type]] - except KeyError: - return default - class CountsIntegralND(Ncm.IntegralND): """Integral subclass used to compute the integrals using NumCosmo.""" @@ -137,10 +58,12 @@ def __init__( self, dim: int, fun: Callable[[npt.NDArray], Sequence[float]], + args: List[float], ) -> None: super().__init__() self.dim = dim self.fun = fun + self.extra_args = args # pylint: disable-next=arguments-differ def do_get_dimensions(self) -> Tuple[int, int]: @@ -158,4 +81,4 @@ def do_integrand( ) -> None: """Called by NumCosmo to evaluate the integrand.""" x = np.array(x_vec.dup_array()).reshape(npoints, dim) - fval_vec.set_array(self.fun(x)) + fval_vec.set_array(self.fun(x, *self.extra_args)) diff --git a/firecrown/models/cluster/integrator/scipy_integrator.py b/firecrown/models/cluster/integrator/scipy_integrator.py index 87bc8420..0edd3bdf 100644 --- a/firecrown/models/cluster/integrator/scipy_integrator.py +++ b/firecrown/models/cluster/integrator/scipy_integrator.py @@ -2,13 +2,9 @@ This module holds the scipy implementation of the integrator classes """ -from typing import Callable, Dict, Tuple -import numpy as np -import numpy.typing as npt +from typing import Callable from scipy.integrate import nquad from firecrown.models.cluster.integrator.integrator import Integrator -from firecrown.models.cluster.kernel import KernelType -from firecrown.models.cluster.abundance import ClusterAbundance, AbundanceIntegrand class ScipyIntegrator(Integrator): @@ -21,80 +17,14 @@ def __init__( self._relative_tolerance = relative_tolerance self._absolute_tolerance = absolute_tolerance - self.integral_args_lkp: Dict[KernelType, int] = self._default_integral_args() - - def _integral_wrapper( - self, - integrand: AbundanceIntegrand, - ) -> Callable[..., float]: - self._validate_integrand(integrand) - - def scipy_integrand(*int_args: float) -> float: - default = -1.0 - # pylint: disable=R0801 - mass = self._get_or_default(int_args, KernelType.MASS, default) - z = self._get_or_default(int_args, KernelType.Z, default) - mass_proxy = self._get_or_default(int_args, KernelType.MASS_PROXY, default) - z_proxy = self._get_or_default(int_args, KernelType.Z_PROXY, default) - - return_val = integrand( - mass, - z, - self.sky_area, - mass_proxy, - z_proxy, - self.mass_proxy_limits, - self.z_proxy_limits, - )[0] - assert isinstance(return_val, float) - return return_val - - return scipy_integrand - - def set_integration_bounds( - self, - cl_abundance: ClusterAbundance, - sky_area: float, - z_proxy_limits: Tuple[float, float], - mass_proxy_limits: Tuple[float, float], - ) -> None: - # pylint: disable=R0801 - self.integral_args_lkp = self._default_integral_args() - self.integral_bounds = [ - (cl_abundance.min_mass, cl_abundance.max_mass), - (cl_abundance.min_z, cl_abundance.max_z), - ] - - self.mass_proxy_limits = mass_proxy_limits - self.z_proxy_limits = z_proxy_limits - self.sky_area = sky_area - - for kernel in cl_abundance.dirac_delta_kernels: - if kernel.kernel_type == KernelType.Z_PROXY: - self.integral_bounds[1] = z_proxy_limits - - elif kernel.kernel_type == KernelType.MASS_PROXY: - self.integral_bounds[0] = mass_proxy_limits - - for kernel in cl_abundance.integrable_kernels: - idx = len(self.integral_bounds) - - if kernel.kernel_type == KernelType.Z_PROXY: - self.integral_bounds.append(z_proxy_limits) - self.integral_args_lkp[KernelType.Z_PROXY] = idx - - elif kernel.kernel_type == KernelType.MASS_PROXY: - self.integral_bounds.append(mass_proxy_limits) - self.integral_args_lkp[KernelType.MASS_PROXY] = idx - def integrate( self, - integrand: AbundanceIntegrand, + func_to_integrate: Callable, ) -> float: - scipy_integrand = self._integral_wrapper(integrand) val = nquad( - scipy_integrand, + func_to_integrate, ranges=self.integral_bounds, + args=self.extra_args, opts={ "epsabs": self._absolute_tolerance, "epsrel": self._relative_tolerance, @@ -102,14 +32,3 @@ def integrate( )[0] assert isinstance(val, float) return val - - def _get_or_default( - self, - int_args: Tuple[float, ...], - kernel_type: KernelType, - default: float, - ) -> npt.NDArray[np.float64]: - try: - return np.atleast_1d(int_args[self.integral_args_lkp[kernel_type]]) - except KeyError: - return np.atleast_1d(default) From d165a708db757b2907a61d75e0751e98ac45b5c1 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 13:41:56 -0800 Subject: [PATCH 02/26] Removing the idea of a kernel from the cluster abundance class, and removing any generic construction of an integrand. Leaving the cluster abundance class as a collections of methods users can use to construct their own cluster "recipes" --- firecrown/models/cluster/abundance.py | 65 +-------------------------- 1 file changed, 1 insertion(+), 64 deletions(-) diff --git a/firecrown/models/cluster/abundance.py b/firecrown/models/cluster/abundance.py index 6d836ff1..dd97bbe2 100644 --- a/firecrown/models/cluster/abundance.py +++ b/firecrown/models/cluster/abundance.py @@ -4,15 +4,13 @@ and phenomenological predictions. This module contains the classes and functions that produce those predictions. """ -from typing import List, Callable, Optional, Dict, Tuple +from typing import Callable, Dict, Tuple import numpy as np import numpy.typing as npt import pyccl import pyccl.background as bkg from pyccl.cosmology import Cosmology from firecrown.updatable import Updatable, UpdatableCollection -from firecrown.models.cluster.kernel import Kernel -from firecrown.models.cluster.properties import ClusterProperty AbundanceIntegrand = Callable[ @@ -43,23 +41,6 @@ def cosmo(self) -> Cosmology: """The cosmology used to predict the cluster number count.""" return self._cosmo - @property - def analytic_kernels(self) -> List[Kernel]: - """The kernels in in the integrand which have an analytic solution.""" - return [x for x in self.kernels if x.has_analytic_sln] - - @property - def dirac_delta_kernels(self) -> List[Kernel]: - """The kernels in in the integrand which are dirac delta functions.""" - return [x for x in self.kernels if x.is_dirac_delta] - - @property - def integrable_kernels(self) -> List[Kernel]: - """The kernels in in the integrand which must be integrated.""" - return [ - x for x in self.kernels if not x.is_dirac_delta and not x.has_analytic_sln - ] - def __init__( self, min_mass: float, @@ -78,10 +59,6 @@ def __init__( self._hmf_cache: Dict[Tuple[float, float], float] = {} self._cosmo: Cosmology = None - def add_kernel(self, kernel: Kernel) -> None: - """Add a kernel to the cluster abundance integrand""" - self.kernels.append(kernel) - def update_ingredients(self, cosmo: Cosmology) -> None: """Update the cluster abundance calculation with a new cosmology.""" self._cosmo = cosmo @@ -126,43 +103,3 @@ def mass_function( return_vals.append(val) return np.asarray(return_vals, dtype=np.float64) - - def get_integrand( - self, *, average_properties: Optional[ClusterProperty] = None - ) -> AbundanceIntegrand: - """Returns a callable that evaluates the complete integrand.""" - - def integrand( - mass: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - sky_area: float, - mass_proxy: npt.NDArray[np.float64], - z_proxy: npt.NDArray[np.float64], - mass_proxy_limits: Tuple[float, float], - z_proxy_limits: Tuple[float, float], - ) -> npt.NDArray[np.float64]: - integrand = self.comoving_volume(z, sky_area) * self.mass_function(mass, z) - - for kernel in self.kernels: - assert isinstance(kernel, Kernel) - integrand *= kernel.distribution( - mass, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits - ) - - if average_properties is None: - return integrand - - for cluster_prop in ClusterProperty: - include_prop = cluster_prop & average_properties - if not include_prop: - continue - if cluster_prop == ClusterProperty.MASS: - integrand *= mass - elif cluster_prop == ClusterProperty.REDSHIFT: - integrand *= z - else: - raise NotImplementedError(f"Average for {cluster_prop}.") - - return integrand - - return integrand From 89828e88e5b4f0fd4f008d0de82714871c51acef Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 13:48:45 -0800 Subject: [PATCH 03/26] Removing Kernel, no longer needed. --- firecrown/models/cluster/kernel.py | 83 +++----------------------- firecrown/models/cluster/mass_proxy.py | 34 +++-------- 2 files changed, 16 insertions(+), 101 deletions(-) diff --git a/firecrown/models/cluster/kernel.py b/firecrown/models/cluster/kernel.py index 3e094269..68b3390d 100644 --- a/firecrown/models/cluster/kernel.py +++ b/firecrown/models/cluster/kernel.py @@ -2,12 +2,10 @@ This module holds the classes that define the kernels that can be included in the cluster abundance integrand.""" -from abc import ABC, abstractmethod from enum import Enum -from typing import List, Tuple, Optional +from typing import Tuple import numpy.typing as npt import numpy as np -from firecrown.updatable import Updatable class KernelType(Enum): @@ -21,55 +19,16 @@ class KernelType(Enum): PURITY = 6 -class Kernel(Updatable, ABC): - """The Kernel base class - - This class defines the common interface that any kernel must implement in order - to be included in the cluster abundance integrand.""" - - def __init__( - self, - kernel_type: KernelType, - is_dirac_delta: bool = False, - has_analytic_sln: bool = False, - integral_bounds: Optional[List[Tuple[float, float]]] = None, - ): - super().__init__() - self.integral_bounds = integral_bounds - self.is_dirac_delta = is_dirac_delta - self.kernel_type = kernel_type - self.has_analytic_sln = has_analytic_sln - - @abstractmethod - def distribution( - self, - mass: npt.NDArray[np.float64], - z: npt.NDArray[np.float64], - mass_proxy: npt.NDArray[np.float64], - z_proxy: npt.NDArray[np.float64], - mass_proxy_limits: Tuple[float, float], - z_proxy_limits: Tuple[float, float], - ) -> npt.NDArray[np.float64]: - """The functional form of the distribution or spread of this kernel""" - - -class Completeness(Kernel): +class Completeness: """The completeness kernel This kernel will affect the integrand by accounting for the incompleteness of a cluster selection.""" - def __init__(self) -> None: - super().__init__(KernelType.COMPLETENESS) - def distribution( self, mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], - _mass_proxy: npt.NDArray[np.float64], - _z_proxy: npt.NDArray[np.float64], - _mass_proxy_limits: Tuple[float, float], - _z_proxy_limits: Tuple[float, float], ) -> npt.NDArray[np.float64]: a_nc = 1.1321 b_nc = 0.7751 @@ -82,15 +41,12 @@ def distribution( return completeness -class Purity(Kernel): +class Purity: """The purity kernel This kernel will affect the integrand by accounting for the purity of a cluster selection.""" - def __init__(self) -> None: - super().__init__(KernelType.PURITY) - def _ln_rc(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: a_rc = 2.2183 b_rc = -0.6592 @@ -106,12 +62,9 @@ def _nc(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: def distribution( self, - _mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], mass_proxy: npt.NDArray[np.float64], - _z_proxy: npt.NDArray[np.float64], mass_proxy_limits: Tuple[float, float], - _z_proxy_limits: Tuple[float, float], ) -> npt.NDArray[np.float64]: if all(mass_proxy == -1.0): mean_mass = (mass_proxy_limits[0] + mass_proxy_limits[1]) / 2 @@ -126,42 +79,20 @@ def distribution( return purity -class TrueMass(Kernel): +class TrueMass: """The true mass kernel. Assuming we measure the true mass, this will always be 1.""" - def __init__(self) -> None: - super().__init__(KernelType.MASS_PROXY, True) - - def distribution( - self, - _mass: npt.NDArray[np.float64], - _z: npt.NDArray[np.float64], - _mass_proxy: npt.NDArray[np.float64], - _z_proxy: npt.NDArray[np.float64], - _mass_proxy_limits: Tuple[float, float], - _z_proxy_limits: Tuple[float, float], - ) -> npt.NDArray[np.float64]: + def distribution(self) -> npt.NDArray[np.float64]: return np.atleast_1d(1.0) -class SpectroscopicRedshift(Kernel): +class SpectroscopicRedshift: """The spec-z kernel. Assuming the spectroscopic redshift has no uncertainties, this is akin to multiplying by 1.""" - def __init__(self) -> None: - super().__init__(KernelType.Z_PROXY, True) - - def distribution( - self, - _mass: npt.NDArray[np.float64], - _z: npt.NDArray[np.float64], - _mass_proxy: npt.NDArray[np.float64], - _z_proxy: npt.NDArray[np.float64], - _mass_proxy_limits: Tuple[float, float], - _z_proxy_limits: Tuple[float, float], - ) -> npt.NDArray[np.float64]: + def distribution(self) -> npt.NDArray[np.float64]: return np.atleast_1d(1.0) diff --git a/firecrown/models/cluster/mass_proxy.py b/firecrown/models/cluster/mass_proxy.py index f64a4f6a..ae14d0eb 100644 --- a/firecrown/models/cluster/mass_proxy.py +++ b/firecrown/models/cluster/mass_proxy.py @@ -9,12 +9,15 @@ import numpy.typing as npt from scipy import special from firecrown import parameters -from firecrown.models.cluster.kernel import Kernel, KernelType +from firecrown.updatable import Updatable -class MassRichnessGaussian(Kernel): +class MassRichnessGaussian(Updatable): """The representation of mass richness relations that are of a gaussian form.""" + def __init__(self, parameter_prefix: str | None = None) -> None: + super().__init__(parameter_prefix) + @staticmethod def observed_value( p: Tuple[float, float, float], @@ -53,10 +56,7 @@ def _distribution_binned( self, mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], - _mass_proxy: npt.NDArray[np.float64], - _z_proxy: npt.NDArray[np.float64], mass_proxy_limits: Tuple[float, float], - _z_proxy_limits: Tuple[float, float], ) -> npt.NDArray[np.float64]: proxy_mean = self.get_proxy_mean(mass, z) proxy_sigma = self.get_proxy_sigma(mass, z) @@ -88,9 +88,6 @@ def _distribution_unbinned( mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], mass_proxy: npt.NDArray[np.float64], - _z_proxy: npt.NDArray[np.float64], - _mass_proxy_limits: Tuple[float, float], - _z_proxy_limits: Tuple[float, float], ) -> npt.NDArray[np.float64]: proxy_mean = self.get_proxy_mean(mass, z) proxy_sigma = self.get_proxy_sigma(mass, z) @@ -111,10 +108,8 @@ def __init__( self, pivot_mass: float, pivot_redshift: float, - integral_bounds: Optional[List[Tuple[float, float]]] = None, ): - super().__init__(KernelType.MASS_PROXY, False, True, integral_bounds) - + super().__init__() self.pivot_redshift = pivot_redshift self.pivot_mass = pivot_mass * np.log(10.0) # ln(M) self.log1p_pivot_redshift = np.log1p(self.pivot_redshift) @@ -159,14 +154,9 @@ def distribution( self, mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], - mass_proxy: npt.NDArray[np.float64], - z_proxy: npt.NDArray[np.float64], mass_proxy_limits: Tuple[float, float], - z_proxy_limits: Tuple[float, float], ) -> npt.NDArray[np.float64]: - return self._distribution_binned( - mass, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits - ) + return self._distribution_binned(mass, z, mass_proxy_limits) class MurataUnbinned(MassRichnessGaussian): @@ -178,8 +168,7 @@ def __init__( pivot_redshift: float, integral_bounds: Optional[List[Tuple[float, float]]] = None, ): - super().__init__(KernelType.MASS_PROXY, False, True, integral_bounds) - + super().__init__() self.pivot_redshift = pivot_redshift self.pivot_mass = pivot_mass * np.log(10.0) # ln(M) self.log1p_pivot_redshift = np.log1p(self.pivot_redshift) @@ -225,10 +214,5 @@ def distribution( mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], mass_proxy: npt.NDArray[np.float64], - z_proxy: npt.NDArray[np.float64], - mass_proxy_limits: Tuple[float, float], - z_proxy_limits: Tuple[float, float], ) -> npt.NDArray[np.float64]: - return self._distribution_unbinned( - mass, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits - ) + return self._distribution_unbinned(mass, z, mass_proxy) From dabad2ff49378d5e2a07b375b504ef571b7be68e Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 14:09:36 -0800 Subject: [PATCH 04/26] Adding new cluster "recipes" which are combinations of cluster ingredients that produce theoretical predictions. --- .../models/cluster/recipes/cluster_recipe.py | 23 +++++ .../cluster/recipes/murata_binned_spec_z.py | 83 +++++++++++++++++++ 2 files changed, 106 insertions(+) create mode 100644 firecrown/models/cluster/recipes/cluster_recipe.py create mode 100644 firecrown/models/cluster/recipes/murata_binned_spec_z.py diff --git a/firecrown/models/cluster/recipes/cluster_recipe.py b/firecrown/models/cluster/recipes/cluster_recipe.py new file mode 100644 index 00000000..05b606f9 --- /dev/null +++ b/firecrown/models/cluster/recipes/cluster_recipe.py @@ -0,0 +1,23 @@ +from abc import ABC, abstractmethod +from typing import Optional + +from firecrown.models.cluster.abundance import ClusterAbundance +from firecrown.models.cluster.binning import SaccBin +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.updatable import Updatable, UpdatableCollection + + +class ClusterRecipe(Updatable, ABC): + def __init__(self, parameter_prefix: str | None = None) -> None: + super().__init__(parameter_prefix) + self.my_updatables: UpdatableCollection = UpdatableCollection() + + @abstractmethod + def evaluate_theory_prediction( + self, + cluster_theory: ClusterAbundance, + bin: SaccBin, + sky_area: float, + average_on: Optional[ClusterProperty] = None, + ) -> float: + pass diff --git a/firecrown/models/cluster/recipes/murata_binned_spec_z.py b/firecrown/models/cluster/recipes/murata_binned_spec_z.py new file mode 100644 index 00000000..cc701da5 --- /dev/null +++ b/firecrown/models/cluster/recipes/murata_binned_spec_z.py @@ -0,0 +1,83 @@ +from typing import Callable, Optional + +from firecrown.models.cluster.abundance import ClusterAbundance +from firecrown.models.cluster.binning import SaccBin +from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator +from firecrown.models.cluster.kernel import SpectroscopicRedshift +from firecrown.models.cluster.mass_proxy import MurataBinned +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe + + +class MurataBinnedSpecZRecipe(ClusterRecipe): + def __init__(self) -> None: + super().__init__() + + self.integrator = NumCosmoIntegrator() + self.redshift_distribution = SpectroscopicRedshift() + pivot_mass, pivot_redshift = 14.625862906, 0.6 + self.mass_distribution = MurataBinned(pivot_mass, pivot_redshift) + self.my_updatables.append(self.mass_distribution) + + def get_theory_prediction( + self, + cluster_theory: ClusterAbundance, + average_on: Optional[ClusterProperty] = None, + ) -> Callable: + """_summary_""" + + def theory_prediction(mass, z, mass_proxy_limits, sky_area): + prediction = ( + cluster_theory.comoving_volume(z, sky_area) + * cluster_theory.mass_function(mass, z) + * self.redshift_distribution.distribution() + * self.mass_distribution.distribution(mass, z, mass_proxy_limits) + ) + + if average_on is None: + return prediction + + for cluster_prop in ClusterProperty: + include_prop = cluster_prop & average_on + if not include_prop: + continue + if cluster_prop == ClusterProperty.MASS: + prediction *= mass + elif cluster_prop == ClusterProperty.REDSHIFT: + prediction *= z + else: + raise NotImplementedError(f"Average for {cluster_prop}.") + + return prediction + + return theory_prediction + + def get_function_to_integrate(self, prediction: Callable) -> Callable: + def numcosmo_wrapper( + int_args, mass_proxy_low, mass_proxy_high, sky_area + ) -> float: + mass = int_args[:, 0] + z = int_args[:, 1] + return prediction(mass, z, (mass_proxy_low, mass_proxy_high), sky_area) + + return numcosmo_wrapper + + def evaluate_theory_prediction( + self, + cluster_theory: ClusterAbundance, + bin: SaccBin, + sky_area: float, + average_on: Optional[ClusterProperty] = None, + ) -> float: + self.integrator.integral_bounds = [ + (cluster_theory.min_mass, cluster_theory.max_mass), + bin.z_edges, + ] + self.integrator.extra_args = [*bin.mass_proxy_edges, sky_area] + + theory_prediction = self.get_theory_prediction(cluster_theory, average_on) + prediction_wrapper = self.get_function_to_integrate(theory_prediction) + + counts = self.integrator.integrate(prediction_wrapper) + + return counts From e5372caa3be938f38d1aeeec356cd6b574ee440e Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 14:10:54 -0800 Subject: [PATCH 05/26] Updating the likelihood script to just create the new recipe. --- .../cluster_redshift_richness.py | 30 ++++++------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/examples/cluster_number_counts/cluster_redshift_richness.py b/examples/cluster_number_counts/cluster_redshift_richness.py index a7acf0b3..2b2c5e02 100644 --- a/examples/cluster_number_counts/cluster_redshift_richness.py +++ b/examples/cluster_number_counts/cluster_redshift_richness.py @@ -1,24 +1,22 @@ """Likelihood factory function for cluster number counts.""" import os +from typing import Tuple import pyccl as ccl import sacc -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.likelihood.likelihood import Likelihood, NamedParameters from firecrown.modeling_tools import ModelingTools from firecrown.models.cluster.abundance import ClusterAbundance -from firecrown.models.cluster.kernel import ( - SpectroscopicRedshift, +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.models.cluster.recipes.murata_binned_spec_z import ( + MurataBinnedSpecZRecipe, ) -from firecrown.models.cluster.mass_proxy import MurataBinned -from firecrown.likelihood.likelihood import NamedParameters, Likelihood -from typing import Tuple def get_cluster_abundance() -> ClusterAbundance: @@ -27,14 +25,6 @@ def get_cluster_abundance() -> ClusterAbundance: min_z, max_z = 0.2, 0.8 cluster_abundance = ClusterAbundance(min_mass, max_mass, min_z, max_z, hmf) - # 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) - - redshift_proxy_kernel = SpectroscopicRedshift() - cluster_abundance.add_kernel(redshift_proxy_kernel) - return cluster_abundance @@ -44,19 +34,17 @@ def build_likelihood( """ 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 + average_on = ClusterProperty.NONE if build_parameters.get_bool("use_cluster_counts", True): - average_properties |= ClusterProperty.COUNTS + average_on |= ClusterProperty.COUNTS if build_parameters.get_bool("use_mean_log_mass", False): - average_properties |= ClusterProperty.MASS + average_on |= ClusterProperty.MASS survey_name = "numcosmo_simulated_redshift_richness" likelihood = ConstGaussian( - [BinnedClusterNumberCounts(average_properties, survey_name, integrator)] + [BinnedClusterNumberCounts(average_on, survey_name, MurataBinnedSpecZRecipe())] ) # Read in sacc data From bd6dd3da27fd457701a6f0770e9fcf28141a6ec6 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 14:12:01 -0800 Subject: [PATCH 06/26] Using cluster recipes in the statistic, simplifying the code. --- .../statistic/binned_cluster_number_counts.py | 55 ++++++++----------- 1 file changed, 23 insertions(+), 32 deletions(-) diff --git a/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py b/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py index 8d5e9e01..305f4b5d 100644 --- a/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py +++ b/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py @@ -4,20 +4,23 @@ clusters within a single redshift and mass bin. """ from __future__ import annotations + from typing import List, Optional -import sacc + import numpy as np -from firecrown.models.cluster.integrator.integrator import Integrator -from firecrown.models.cluster.abundance_data import AbundanceData -from firecrown.models.cluster.binning import SaccBin -from firecrown.models.cluster.properties import ClusterProperty +import sacc + +from firecrown.likelihood.gauss_family.statistic.source.source import SourceSystematic from firecrown.likelihood.gauss_family.statistic.statistic import ( - Statistic, DataVector, + Statistic, TheoryVector, ) -from firecrown.likelihood.gauss_family.statistic.source.source import SourceSystematic from firecrown.modeling_tools import ModelingTools +from firecrown.models.cluster.abundance_data import AbundanceData +from firecrown.models.cluster.binning import SaccBin +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe class BinnedClusterNumberCounts(Statistic): @@ -31,7 +34,7 @@ def __init__( self, cluster_properties: ClusterProperty, survey_name: str, - integrator: Integrator, + cluster_recipe: ClusterRecipe, systematics: Optional[List[SourceSystematic]] = None, ): super().__init__() @@ -39,7 +42,7 @@ def __init__( self.theory_vector: Optional[TheoryVector] = None self.cluster_properties = cluster_properties self.survey_name = survey_name - self.integrator = integrator + self.cluster_recipe = cluster_recipe self.data_vector = DataVector.from_list([]) self.sky_area = 0.0 self.bins: List[SaccBin] = [] @@ -110,23 +113,17 @@ def get_binned_cluster_property( the clusters in each bin.""" assert tools.cluster_abundance is not None - cluster_masses = [] - for bin_edge, counts in zip(self.bins, cluster_counts): - integrand = tools.cluster_abundance.get_integrand( - average_properties=cluster_properties - ) - self.integrator.set_integration_bounds( - tools.cluster_abundance, - self.sky_area, - bin_edge.z_edges, - bin_edge.mass_proxy_edges, + mean_values = [] + for this_bin, counts in zip(self.bins, cluster_counts): + total_observable = self.cluster_recipe.evaluate_theory_prediction( + tools.cluster_abundance, this_bin, self.sky_area, cluster_properties ) + cluster_counts.append(counts) - total_mass = self.integrator.integrate(integrand) - mean_mass = total_mass / counts - cluster_masses.append(mean_mass) + mean_observable = total_observable / counts + mean_values.append(mean_observable) - return cluster_masses + return mean_values def get_binned_cluster_counts(self, tools: ModelingTools) -> List[float]: """Computes the number of clusters in each bin @@ -137,16 +134,10 @@ def get_binned_cluster_counts(self, tools: ModelingTools) -> List[float]: assert tools.cluster_abundance is not None cluster_counts = [] - for bin_edge in self.bins: - self.integrator.set_integration_bounds( - tools.cluster_abundance, - self.sky_area, - bin_edge.z_edges, - bin_edge.mass_proxy_edges, + for this_bin in self.bins: + counts = self.cluster_recipe.evaluate_theory_prediction( + tools.cluster_abundance, this_bin, self.sky_area ) - - integrand = tools.cluster_abundance.get_integrand() - counts = self.integrator.integrate(integrand) cluster_counts.append(counts) return cluster_counts From 4918b63dcc794cf7ab0f556b67546e0638db9748 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 14:31:26 -0800 Subject: [PATCH 07/26] Updating the recipes to be able to handle the base NDimensional bin (needed for comparison script). --- firecrown/models/cluster/binning.py | 41 +++++++++++++++++++ .../cluster/recipes/murata_binned_spec_z.py | 4 +- 2 files changed, 43 insertions(+), 2 deletions(-) diff --git a/firecrown/models/cluster/binning.py b/firecrown/models/cluster/binning.py index f4183b71..6bb9a99b 100644 --- a/firecrown/models/cluster/binning.py +++ b/firecrown/models/cluster/binning.py @@ -84,3 +84,44 @@ def __hash__(self) -> int: """One bin's hash is determined by the dimension and lower/upper bound.""" bin_bounds = [(bin.lower, bin.upper) for bin in self.coordinate_bins] return hash((self.dimension, tuple(bin_bounds))) + + +class TupleBin(NDimensionalBin[Tuple]): + """An implementation of the N dimensional bin using sacc tracers.""" + + def __init__(self, bins: List[Tuple]): + super().__init__(bins) + + @property + def mass_proxy_edges(self) -> Tuple[float, float]: + mass_bin = self.coordinate_bins[0] + return mass_bin[0], mass_bin[1] + + @property + def z_edges(self) -> Tuple[float, float]: + z_bin = self.coordinate_bins[1] + return z_bin[0], z_bin[1] + + def __eq__(self, other: object) -> bool: + """Two bins are equal if they have the same lower/upper bound.""" + if not isinstance(other, TupleBin): + return False + + if self.dimension != other.dimension: + return False + + for i, my_bin in enumerate(self.coordinate_bins): + other_bin = other.coordinate_bins[i] + if len(my_bin) != len(other_bin): + return False + if my_bin[0] != other_bin[0]: + return False + if my_bin[1] != other_bin[1]: + return False + + return True + + def __hash__(self) -> int: + """One bin's hash is determined by the dimension and lower/upper bound.""" + bin_bounds = [(bin[0], bin[1]) for bin in self.coordinate_bins] + return hash((self.dimension, tuple(bin_bounds))) diff --git a/firecrown/models/cluster/recipes/murata_binned_spec_z.py b/firecrown/models/cluster/recipes/murata_binned_spec_z.py index cc701da5..881a09ce 100644 --- a/firecrown/models/cluster/recipes/murata_binned_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_binned_spec_z.py @@ -1,7 +1,7 @@ from typing import Callable, Optional from firecrown.models.cluster.abundance import ClusterAbundance -from firecrown.models.cluster.binning import SaccBin +from firecrown.models.cluster.binning import NDimensionalBin from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator from firecrown.models.cluster.kernel import SpectroscopicRedshift from firecrown.models.cluster.mass_proxy import MurataBinned @@ -65,7 +65,7 @@ def numcosmo_wrapper( def evaluate_theory_prediction( self, cluster_theory: ClusterAbundance, - bin: SaccBin, + bin: NDimensionalBin, sky_area: float, average_on: Optional[ClusterProperty] = None, ) -> float: From c852b159ab36ecd444cdf529ce49722345603d76 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 14:31:47 -0800 Subject: [PATCH 08/26] Fixing comparison script Note- H_V seems to perform best now, need to compare to old numbers. --- .../compare_integration_methods.py | 55 ++++++------------- .../cluster/integrator/numcosmo_integrator.py | 4 +- 2 files changed, 20 insertions(+), 39 deletions(-) diff --git a/examples/cluster_number_counts/compare_integration_methods.py b/examples/cluster_number_counts/compare_integration_methods.py index fb556340..b8d216bc 100644 --- a/examples/cluster_number_counts/compare_integration_methods.py +++ b/examples/cluster_number_counts/compare_integration_methods.py @@ -4,14 +4,14 @@ import pyccl import numpy as np -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 +from firecrown.models.cluster.recipes.murata_binned_spec_z import ( + MurataBinnedSpecZRecipe, +) +from firecrown.models.cluster.binning import TupleBin def get_cosmology() -> pyccl.Cosmology: @@ -44,64 +44,45 @@ def get_cosmology() -> pyccl.Cosmology: return cosmo_ccl -def get_mass_richness() -> Kernel: - pivot_mass = 14.625862906 - pivot_redshift = 0.6 - - 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 - - def compare_integration() -> None: """Compare integration methods.""" hmf = pyccl.halos.MassFuncTinker08() abundance = ClusterAbundance(13, 15, 0, 4, hmf) + cluster_recipe = MurataBinnedSpecZRecipe() - mass_richness = get_mass_richness() - abundance.add_kernel(mass_richness) - - redshift_proxy_kernel = SpectroscopicRedshift() - abundance.add_kernel(redshift_proxy_kernel) + cluster_recipe.mass_distribution.mu_p0 = 3.0 + cluster_recipe.mass_distribution.mu_p1 = 0.86 + cluster_recipe.mass_distribution.mu_p2 = 0.0 + cluster_recipe.mass_distribution.sigma_p0 = 3.0 + cluster_recipe.mass_distribution.sigma_p1 = 0.7 + cluster_recipe.mass_distribution.sigma_p2 = 0.0 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)) + cluster_recipe.integrator.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_limits = ( mass_proxy_bins[mass_proxy_idx], mass_proxy_bins[mass_proxy_idx + 1], ) + z_limits = (z_bins[z_idx], z_bins[z_idx + 1]) - nc_integrator.set_integration_bounds( - abundance, - sky_area, - z_proxy_limits, - mass_proxy_limits, - ) + bin = TupleBin([mass_limits, z_limits]) - counts = nc_integrator.integrate(integrand) + counts = cluster_recipe.evaluate_theory_prediction(abundance, bin, sky_area) counts_list.append(counts) t_stop = time.time() diff --git a/firecrown/models/cluster/integrator/numcosmo_integrator.py b/firecrown/models/cluster/integrator/numcosmo_integrator.py index 487a1935..7f80981b 100644 --- a/firecrown/models/cluster/integrator/numcosmo_integrator.py +++ b/firecrown/models/cluster/integrator/numcosmo_integrator.py @@ -28,11 +28,11 @@ def __init__( absolute_tolerance: float = 1e-12, ) -> None: super().__init__() + self.method = method or NumCosmoIntegralMethod.P_V self.integral_bounds: List[Tuple[float, float]] = [] self.extra_args: List[float] = [] self._relative_tolerance = relative_tolerance self._absolute_tolerance = absolute_tolerance - self._method = method or NumCosmoIntegralMethod.P_V def integrate(self, func_to_integrate: Callable) -> float: Ncm.cfg_init() @@ -40,7 +40,7 @@ def integrate(self, func_to_integrate: Callable) -> float: int_nd = CountsIntegralND( len(self.integral_bounds), func_to_integrate, self.extra_args ) - int_nd.set_method(self._method.value) + int_nd.set_method(self.method.value) int_nd.set_reltol(self._relative_tolerance) int_nd.set_abstol(self._absolute_tolerance) res = Ncm.Vector.new(1) From dc4582852bc9b78d925169f39627bd5d32a64581 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 14:58:25 -0800 Subject: [PATCH 09/26] Adding a true_mass_spec_z recipe for proof of concept. --- .../cluster_redshift_richness.py | 2 + .../cluster/recipes/true_mass_spec_z.py | 82 +++++++++++++++++++ 2 files changed, 84 insertions(+) create mode 100644 firecrown/models/cluster/recipes/true_mass_spec_z.py diff --git a/examples/cluster_number_counts/cluster_redshift_richness.py b/examples/cluster_number_counts/cluster_redshift_richness.py index 2b2c5e02..82ffb791 100644 --- a/examples/cluster_number_counts/cluster_redshift_richness.py +++ b/examples/cluster_number_counts/cluster_redshift_richness.py @@ -17,6 +17,7 @@ from firecrown.models.cluster.recipes.murata_binned_spec_z import ( MurataBinnedSpecZRecipe, ) +from firecrown.models.cluster.recipes.true_mass_spec_z import TrueMassSpecZRecipe def get_cluster_abundance() -> ClusterAbundance: @@ -45,6 +46,7 @@ def build_likelihood( survey_name = "numcosmo_simulated_redshift_richness" likelihood = ConstGaussian( [BinnedClusterNumberCounts(average_on, survey_name, MurataBinnedSpecZRecipe())] + # [BinnedClusterNumberCounts(average_on, survey_name, TrueMassSpecZRecipe())] ) # Read in sacc data diff --git a/firecrown/models/cluster/recipes/true_mass_spec_z.py b/firecrown/models/cluster/recipes/true_mass_spec_z.py new file mode 100644 index 00000000..35fd8714 --- /dev/null +++ b/firecrown/models/cluster/recipes/true_mass_spec_z.py @@ -0,0 +1,82 @@ +from typing import Callable, Optional + +from firecrown.models.cluster.abundance import ClusterAbundance +from firecrown.models.cluster.binning import NDimensionalBin +from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator +from firecrown.models.cluster.kernel import SpectroscopicRedshift, TrueMass +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe + + +class TrueMassSpecZRecipe(ClusterRecipe): + def __init__(self) -> None: + super().__init__() + + self.integrator = NumCosmoIntegrator() + self.redshift_distribution = SpectroscopicRedshift() + self.mass_distribution = TrueMass() + + def get_theory_prediction( + self, + cluster_theory: ClusterAbundance, + average_on: Optional[ClusterProperty] = None, + ) -> Callable: + """_summary_""" + + def theory_prediction(mass, z, sky_area): + prediction = ( + cluster_theory.comoving_volume(z, sky_area) + * cluster_theory.mass_function(mass, z) + * self.redshift_distribution.distribution() + * self.mass_distribution.distribution() + ) + + if average_on is None: + return prediction + + for cluster_prop in ClusterProperty: + include_prop = cluster_prop & average_on + if not include_prop: + continue + if cluster_prop == ClusterProperty.MASS: + prediction *= mass + elif cluster_prop == ClusterProperty.REDSHIFT: + prediction *= z + else: + raise NotImplementedError(f"Average for {cluster_prop}.") + + return prediction + + return theory_prediction + + def get_function_to_integrate(self, prediction: Callable) -> Callable: + def numcosmo_wrapper(int_args, sky_area) -> float: + mass = int_args[:, 0] + z = int_args[:, 1] + return prediction(mass, z, sky_area) + + return numcosmo_wrapper + + def evaluate_theory_prediction( + self, + cluster_theory: ClusterAbundance, + bin: NDimensionalBin, + sky_area: float, + average_on: Optional[ClusterProperty] = None, + ) -> float: + # Mock some fake mass to richness relation: + mass_low = bin.mass_proxy_edges[0] + 13 + mass_high = bin.mass_proxy_edges[1] + 13 + self.integrator.integral_bounds = [ + (mass_low, mass_high), + bin.z_edges, + ] + print(mass_low, mass_high) + self.integrator.extra_args = [sky_area] + + theory_prediction = self.get_theory_prediction(cluster_theory, average_on) + prediction_wrapper = self.get_function_to_integrate(theory_prediction) + + counts = self.integrator.integrate(prediction_wrapper) + + return counts From 5c65afc8283890efb105288bcbcf8ffcac91360c Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 15:09:59 -0800 Subject: [PATCH 10/26] Making mypy, flake8, and black happy. --- .../cluster_redshift_richness.py | 3 ++- firecrown/models/cluster/integrator/integrator.py | 6 +++--- .../cluster/integrator/numcosmo_integrator.py | 8 ++++++-- .../models/cluster/integrator/scipy_integrator.py | 1 + firecrown/models/cluster/kernel.py | 14 ++++++++++++-- firecrown/models/cluster/mass_proxy.py | 10 +++++----- 6 files changed, 29 insertions(+), 13 deletions(-) diff --git a/examples/cluster_number_counts/cluster_redshift_richness.py b/examples/cluster_number_counts/cluster_redshift_richness.py index 82ffb791..a259529c 100644 --- a/examples/cluster_number_counts/cluster_redshift_richness.py +++ b/examples/cluster_number_counts/cluster_redshift_richness.py @@ -17,7 +17,8 @@ from firecrown.models.cluster.recipes.murata_binned_spec_z import ( MurataBinnedSpecZRecipe, ) -from firecrown.models.cluster.recipes.true_mass_spec_z import TrueMassSpecZRecipe + +# from firecrown.models.cluster.recipes.true_mass_spec_z import TrueMassSpecZRecipe def get_cluster_abundance() -> ClusterAbundance: diff --git a/firecrown/models/cluster/integrator/integrator.py b/firecrown/models/cluster/integrator/integrator.py index 59caa4b1..a3031125 100644 --- a/firecrown/models/cluster/integrator/integrator.py +++ b/firecrown/models/cluster/integrator/integrator.py @@ -1,17 +1,17 @@ """The integrator module This module holds the classes that define the interface required to -integrate an assembled cluster abundance. +integrate a function. """ from abc import ABC, abstractmethod from typing import Tuple, List, Callable +# pylint: disable=too-few-public-methods class Integrator(ABC): """The integrator base class - This class acts as an adapter around an integration library, and must provides - a specific set of methods to be used to integrate a cluster abundance integral.""" + This class acts as an adapter around an integration library.""" def __init__(self) -> None: self.integral_bounds: List[Tuple[float, float]] = [] diff --git a/firecrown/models/cluster/integrator/numcosmo_integrator.py b/firecrown/models/cluster/integrator/numcosmo_integrator.py index 7f80981b..038ade7e 100644 --- a/firecrown/models/cluster/integrator/numcosmo_integrator.py +++ b/firecrown/models/cluster/integrator/numcosmo_integrator.py @@ -2,13 +2,17 @@ This module holds the NumCosmo implementation of the integrator classes """ -from typing import Tuple, Callable, Sequence, Optional, List from enum import Enum +from typing import Callable, List, Optional, Sequence, Tuple + import numpy as np import numpy.typing as npt from numcosmo_py import Ncm +from firecrown.models.cluster.integrator.integrator import Integrator + +# pylint: disable=too-few-public-methods class NumCosmoIntegralMethod(Enum): """The available NumCosmo integration methods.""" @@ -18,7 +22,7 @@ class NumCosmoIntegralMethod(Enum): H_V = Ncm.IntegralNDMethod.H_V -class NumCosmoIntegrator: +class NumCosmoIntegrator(Integrator): """The NumCosmo implementation of the Integrator base class.""" def __init__( diff --git a/firecrown/models/cluster/integrator/scipy_integrator.py b/firecrown/models/cluster/integrator/scipy_integrator.py index 0edd3bdf..40befb24 100644 --- a/firecrown/models/cluster/integrator/scipy_integrator.py +++ b/firecrown/models/cluster/integrator/scipy_integrator.py @@ -7,6 +7,7 @@ from firecrown.models.cluster.integrator.integrator import Integrator +# pylint: disable=too-few-public-methods class ScipyIntegrator(Integrator): """The scipy implementation of the Integrator base class using nquad.""" diff --git a/firecrown/models/cluster/kernel.py b/firecrown/models/cluster/kernel.py index 68b3390d..df3fcd08 100644 --- a/firecrown/models/cluster/kernel.py +++ b/firecrown/models/cluster/kernel.py @@ -19,8 +19,9 @@ class KernelType(Enum): PURITY = 6 +# pylint: disable=too-few-public-methods class Completeness: - """The completeness kernel + """The completeness kernel for the numcosmo simulated survey This kernel will affect the integrand by accounting for the incompleteness of a cluster selection.""" @@ -30,6 +31,7 @@ def distribution( mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: + """Evaluates and returns the completeness contribution to the integrand.""" a_nc = 1.1321 b_nc = 0.7751 a_mc = 13.31 @@ -41,8 +43,9 @@ def distribution( return completeness +# pylint: disable=too-few-public-methods class Purity: - """The purity kernel + """The purity kernel for the numcosmo simulated survey This kernel will affect the integrand by accounting for the purity of a cluster selection.""" @@ -66,6 +69,7 @@ def distribution( mass_proxy: npt.NDArray[np.float64], mass_proxy_limits: Tuple[float, float], ) -> npt.NDArray[np.float64]: + """Evaluates and returns the purity contribution to the integrand.""" if all(mass_proxy == -1.0): mean_mass = (mass_proxy_limits[0] + mass_proxy_limits[1]) / 2 ln_r = np.log(10**mean_mass) @@ -79,15 +83,19 @@ def distribution( return purity +# pylint: disable=too-few-public-methods class TrueMass: """The true mass kernel. Assuming we measure the true mass, this will always be 1.""" def distribution(self) -> npt.NDArray[np.float64]: + """Evaluates and returns the mass distribution contribution to the integrand. + We have set this to 1.0 (i.e. it does not affect the mass distribution)""" return np.atleast_1d(1.0) +# pylint: disable=too-few-public-methods class SpectroscopicRedshift: """The spec-z kernel. @@ -95,4 +103,6 @@ class SpectroscopicRedshift: multiplying by 1.""" def distribution(self) -> npt.NDArray[np.float64]: + """Evaluates and returns the z distribution contribution to the integrand. + We have set this to 1.0 (i.e. it does not affect the redshift distribution)""" return np.atleast_1d(1.0) diff --git a/firecrown/models/cluster/mass_proxy.py b/firecrown/models/cluster/mass_proxy.py index ae14d0eb..0125c83f 100644 --- a/firecrown/models/cluster/mass_proxy.py +++ b/firecrown/models/cluster/mass_proxy.py @@ -3,11 +3,13 @@ This module holds the classes that define the mass richness relations that can be included in the cluster abundance integrand. These are implementations of Kernels.""" -from typing import List, Tuple, Optional from abc import abstractmethod +from typing import Tuple + import numpy as np import numpy.typing as npt from scipy import special + from firecrown import parameters from firecrown.updatable import Updatable @@ -15,9 +17,6 @@ class MassRichnessGaussian(Updatable): """The representation of mass richness relations that are of a gaussian form.""" - def __init__(self, parameter_prefix: str | None = None) -> None: - super().__init__(parameter_prefix) - @staticmethod def observed_value( p: Tuple[float, float, float], @@ -156,6 +155,7 @@ def distribution( z: npt.NDArray[np.float64], mass_proxy_limits: Tuple[float, float], ) -> npt.NDArray[np.float64]: + """Evaluates and returns the mass-richness contribution to the integrand.""" return self._distribution_binned(mass, z, mass_proxy_limits) @@ -166,7 +166,6 @@ def __init__( self, pivot_mass: float, pivot_redshift: float, - integral_bounds: Optional[List[Tuple[float, float]]] = None, ): super().__init__() self.pivot_redshift = pivot_redshift @@ -215,4 +214,5 @@ def distribution( z: npt.NDArray[np.float64], mass_proxy: npt.NDArray[np.float64], ) -> npt.NDArray[np.float64]: + """Evaluates and returns the mass-richness contribution to the integrand.""" return self._distribution_unbinned(mass, z, mass_proxy) From 0f5e0e14253320cb6f154c78d52fb97aee6773fd Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 15:36:55 -0800 Subject: [PATCH 11/26] Removing old tests. --- tests/test_cluster_abundance.py | 228 ++++++++-------------- tests/test_cluster_integrators.py | 290 +++------------------------- tests/test_cluster_kernels.py | 83 ++------ tests/test_cluster_mass_richness.py | 32 +-- 4 files changed, 116 insertions(+), 517 deletions(-) diff --git a/tests/test_cluster_abundance.py b/tests/test_cluster_abundance.py index f811d873..24be2121 100644 --- a/tests/test_cluster_abundance.py +++ b/tests/test_cluster_abundance.py @@ -1,11 +1,8 @@ """Tests for the cluster abundance module.""" -from unittest.mock import Mock import pytest import pyccl import numpy as np from firecrown.models.cluster.abundance import ClusterAbundance -from firecrown.models.cluster.kernel import Kernel, KernelType -from firecrown.models.cluster.properties import ClusterProperty @pytest.fixture(name="cluster_abundance") @@ -16,42 +13,6 @@ def fixture_cluster_abundance(): return ca -@pytest.fixture(name="dirac_delta_kernel") -def fixture_dirac_delta_kernel(): - mk = Mock( - spec=Kernel, - kernel_type=KernelType.MASS, - is_dirac_delta=True, - has_analytic_sln=False, - ) - mk.distribution.return_value = np.atleast_1d(1.0) - return mk - - -@pytest.fixture(name="integrable_kernel") -def fixture_integrable_kernel(): - mk = Mock( - spec=Kernel, - kernel_type=KernelType.MASS_PROXY, - is_dirac_delta=False, - has_analytic_sln=False, - ) - mk.distribution.return_value = np.atleast_1d(1.0) - return mk - - -@pytest.fixture(name="analytic_sln_kernel") -def fixture_analytic_sln_solved_kernel(): - mk = Mock( - spec=Kernel, - kernel_type=KernelType.MASS, - is_dirac_delta=False, - has_analytic_sln=True, - ) - mk.distribution.return_value = np.atleast_1d(1.0) - return mk - - @pytest.fixture(name="integrand_args") def fixture_integrand_args(): sky_area = 439.78986 @@ -88,42 +49,6 @@ def test_cluster_update_ingredients(cluster_abundance: ClusterAbundance): assert cluster_abundance._hmf_cache == {} -def test_cluster_add_integrable_kernel( - cluster_abundance: ClusterAbundance, integrable_kernel: Kernel -): - cluster_abundance.add_kernel(integrable_kernel) - assert len(cluster_abundance.kernels) == 1 - assert isinstance(cluster_abundance.kernels[0], Kernel) - - assert len(cluster_abundance.analytic_kernels) == 0 - assert len(cluster_abundance.dirac_delta_kernels) == 0 - assert len(cluster_abundance.integrable_kernels) == 1 - - -def test_cluster_add_dirac_delta_kernel( - cluster_abundance: ClusterAbundance, dirac_delta_kernel: Kernel -): - cluster_abundance.add_kernel(dirac_delta_kernel) - assert len(cluster_abundance.kernels) == 1 - assert isinstance(cluster_abundance.kernels[0], Kernel) - - assert len(cluster_abundance.analytic_kernels) == 0 - assert len(cluster_abundance.dirac_delta_kernels) == 1 - assert len(cluster_abundance.integrable_kernels) == 0 - - -def test_cluster_add_analytic_sln_kernel( - cluster_abundance: ClusterAbundance, analytic_sln_kernel: Kernel -): - cluster_abundance.add_kernel(analytic_sln_kernel) - assert len(cluster_abundance.kernels) == 1 - assert isinstance(cluster_abundance.kernels[0], Kernel) - - assert len(cluster_abundance.analytic_kernels) == 1 - assert len(cluster_abundance.dirac_delta_kernels) == 0 - assert len(cluster_abundance.integrable_kernels) == 0 - - def test_abundance_comoving_returns_value(cluster_abundance: ClusterAbundance): cosmo = pyccl.CosmologyVanillaLCDM() cluster_abundance.update_ingredients(cosmo) @@ -149,79 +74,80 @@ def test_abundance_massfunc_returns_value(cluster_abundance: ClusterAbundance): assert np.all(result > 0) -@pytest.mark.slow -def test_abundance_get_integrand( - cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -): - cosmo = pyccl.CosmologyVanillaLCDM() - cluster_abundance.update_ingredients(cosmo) - cluster_abundance.add_kernel(integrable_kernel) - - integrand = cluster_abundance.get_integrand() - assert callable(integrand) - result = integrand(*integrand_args) - assert isinstance(result, np.ndarray) - assert np.issubdtype(result.dtype, np.float64) - - -@pytest.mark.slow -def test_abundance_get_integrand_avg_mass( - cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -): - cosmo = pyccl.CosmologyVanillaLCDM() - cluster_abundance.update_ingredients(cosmo) - cluster_abundance.add_kernel(integrable_kernel) - - integrand = cluster_abundance.get_integrand(average_properties=ClusterProperty.MASS) - assert callable(integrand) - result = integrand(*integrand_args) - assert isinstance(result, np.ndarray) - assert np.issubdtype(result.dtype, np.float64) - - -@pytest.mark.slow -def test_abundance_get_integrand_avg_redshift( - cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -): - cosmo = pyccl.CosmologyVanillaLCDM() - cluster_abundance.update_ingredients(cosmo) - cluster_abundance.add_kernel(integrable_kernel) - - integrand = cluster_abundance.get_integrand( - average_properties=ClusterProperty.REDSHIFT - ) - assert callable(integrand) - result = integrand(*integrand_args) - assert isinstance(result, np.ndarray) - assert np.issubdtype(result.dtype, np.float64) - - -@pytest.mark.slow -def test_abundance_get_integrand_avg_mass_and_redshift( - cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -): - cosmo = pyccl.CosmologyVanillaLCDM() - cluster_abundance.update_ingredients(cosmo) - cluster_abundance.add_kernel(integrable_kernel) - - average_properties = ClusterProperty.MASS | ClusterProperty.REDSHIFT - integrand = cluster_abundance.get_integrand(average_properties=average_properties) - assert callable(integrand) - result = integrand(*integrand_args) - assert isinstance(result, np.ndarray) - assert np.issubdtype(result.dtype, np.float64) - - -@pytest.mark.slow -def test_abundance_get_integrand_avg_not_implemented_throws( - cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -): - cosmo = pyccl.CosmologyVanillaLCDM() - cluster_abundance.update_ingredients(cosmo) - cluster_abundance.add_kernel(integrable_kernel) - - average_properties = ClusterProperty.SHEAR - integrand = cluster_abundance.get_integrand(average_properties=average_properties) - assert callable(integrand) - with pytest.raises(NotImplementedError): - _ = integrand(*integrand_args) +# @pytest.mark.slow +# def test_abundance_get_integrand( +# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args +# ): +# cosmo = pyccl.CosmologyVanillaLCDM() +# cluster_abundance.update_ingredients(cosmo) +# cluster_abundance.add_kernel(integrable_kernel) + +# integrand = cluster_abundance.get_integrand() +# assert callable(integrand) +# result = integrand(*integrand_args) +# assert isinstance(result, np.ndarray) +# assert np.issubdtype(result.dtype, np.float64) + + +# @pytest.mark.slow +# def test_abundance_get_integrand_avg_mass( +# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args +# ): +# cosmo = pyccl.CosmologyVanillaLCDM() +# cluster_abundance.update_ingredients(cosmo) +# cluster_abundance.add_kernel(integrable_kernel) + +# integrand = cluster_abundance.get_integrand(average_properties= +# ClusterProperty.MASS) +# assert callable(integrand) +# result = integrand(*integrand_args) +# assert isinstance(result, np.ndarray) +# assert np.issubdtype(result.dtype, np.float64) + + +# @pytest.mark.slow +# def test_abundance_get_integrand_avg_redshift( +# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args +# ): +# cosmo = pyccl.CosmologyVanillaLCDM() +# cluster_abundance.update_ingredients(cosmo) +# cluster_abundance.add_kernel(integrable_kernel) + +# integrand = cluster_abundance.get_integrand( +# average_properties=ClusterProperty.REDSHIFT +# ) +# assert callable(integrand) +# result = integrand(*integrand_args) +# assert isinstance(result, np.ndarray) +# assert np.issubdtype(result.dtype, np.float64) + + +# @pytest.mark.slow +# def test_abundance_get_integrand_avg_mass_and_redshift( +# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args +# ): +# cosmo = pyccl.CosmologyVanillaLCDM() +# cluster_abundance.update_ingredients(cosmo) +# cluster_abundance.add_kernel(integrable_kernel) + +# average_properties = ClusterProperty.MASS | ClusterProperty.REDSHIFT +# integrand = cluster_abundance.get_integrand(average_properties=average_properties) +# assert callable(integrand) +# result = integrand(*integrand_args) +# assert isinstance(result, np.ndarray) +# assert np.issubdtype(result.dtype, np.float64) + + +# @pytest.mark.slow +# def test_abundance_get_integrand_avg_not_implemented_throws( +# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args +# ): +# cosmo = pyccl.CosmologyVanillaLCDM() +# cluster_abundance.update_ingredients(cosmo) +# cluster_abundance.add_kernel(integrable_kernel) + +# average_properties = ClusterProperty.SHEAR +# integrand = cluster_abundance.get_integrand(average_properties=average_properties) +# assert callable(integrand) +# with pytest.raises(NotImplementedError): +# _ = integrand(*integrand_args) diff --git a/tests/test_cluster_integrators.py b/tests/test_cluster_integrators.py index acaf31a3..5ae7e512 100644 --- a/tests/test_cluster_integrators.py +++ b/tests/test_cluster_integrators.py @@ -1,291 +1,45 @@ """Tests for the integrator module.""" -from typing import Tuple -from unittest.mock import Mock + import numpy as np import numpy.typing as npt import pytest from firecrown.models.cluster.integrator.integrator import Integrator from firecrown.models.cluster.integrator.scipy_integrator import ScipyIntegrator from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator -from firecrown.models.cluster.kernel import KernelType, Kernel -from firecrown.models.cluster.abundance import ClusterAbundance - - -@pytest.fixture(name="integrator", params=[ScipyIntegrator, NumCosmoIntegrator]) -def fixture_integrator(request) -> Integrator: - return request.param() - - -def test_set_integration_bounds_no_kernels( - empty_cluster_abundance: ClusterAbundance, integrator: Integrator -): - z_array = np.linspace(0, 2, 10) - m_array = np.linspace(13, 17, 10) - z_bins = list(zip(z_array[:-1], z_array[1:])) - m_bins = list(zip(m_array[:-1], m_array[1:])) - sky_area = 100**2 - - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) - - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 2 - assert integrator.integral_bounds == [(13, 17), (0, 2)] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - } - - -def test_set_integration_bounds_dirac_delta( - empty_cluster_abundance: ClusterAbundance, integrator: Integrator -): - z_array = np.linspace(0, 2, 10) - m_array = np.linspace(13, 17, 10) - z_bins = list(zip(z_array[:-1], z_array[1:])) - m_bins = list(zip(m_array[:-1], m_array[1:])) - sky_area = 100**2 - - dd_kernel = Mock( - spec=Kernel, - kernel_type=KernelType.MASS_PROXY, - is_dirac_delta=True, - has_analytic_sln=False, - ) - dd_kernel.distribution.return_value = np.atleast_1d(1.0) - empty_cluster_abundance.add_kernel(dd_kernel) - - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) - - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 2 - assert integrator.integral_bounds == [mass_limits, (0, 2)] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - } - - empty_cluster_abundance.kernels.clear() - dd_kernel = Mock( - spec=Kernel, - kernel_type=KernelType.Z_PROXY, - has_analytic_sln=False, - is_dirac_delta=True, - ) - dd_kernel.distribution.return_value = np.atleast_1d(1.0) - empty_cluster_abundance.add_kernel(dd_kernel) - - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) - - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 2 - assert integrator.integral_bounds == [(13, 17), z_limits] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - } - dd_kernel2 = Mock( - spec=Kernel, - kernel_type=KernelType.MASS_PROXY, - has_analytic_sln=False, - is_dirac_delta=True, - ) - dd_kernel2.distribution.return_value = np.atleast_1d(1.0) - empty_cluster_abundance.add_kernel(dd_kernel2) - - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) - - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 2 - assert integrator.integral_bounds == [mass_limits, z_limits] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - } - -def test_set_integration_bounds_integrable_kernels( - empty_cluster_abundance: ClusterAbundance, integrator: Integrator -): - z_array = np.linspace(0, 2, 10) - m_array = np.linspace(13, 17, 10) - z_bins = list(zip(z_array[:-1], z_array[1:])) - m_bins = list(zip(m_array[:-1], m_array[1:])) - sky_area = 100**2 - ig_kernel = Mock( - spec=Kernel, - kernel_type=KernelType.MASS_PROXY, - has_analytic_sln=False, - is_dirac_delta=False, - ) - ig_kernel.distribution.return_value = np.atleast_1d(1.0) - empty_cluster_abundance.add_kernel(ig_kernel) +# @pytest.fixture(name="integrator", params=[ScipyIntegrator, NumCosmoIntegrator]) +# def fixture_integrator(request) -> Integrator: +# return request.param() - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 3 - assert integrator.integral_bounds == [(13, 17), (0, 2), mass_limits] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - KernelType.MASS_PROXY: 2, - } +def test_numcosmo_integrator_integrate(): + integrator = NumCosmoIntegrator() - empty_cluster_abundance.kernels.clear() - ig_kernel = Mock( - spec=Kernel, - kernel_type=KernelType.Z_PROXY, - has_analytic_sln=False, - is_dirac_delta=False, - ) - ig_kernel.distribution.return_value = np.atleast_1d(1.0) - empty_cluster_abundance.add_kernel(ig_kernel) - - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) - - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 3 - assert integrator.integral_bounds == [(13, 17), (0, 2), z_limits] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - KernelType.Z_PROXY: 2, - } - - ig_kernel2 = Mock( - spec=Kernel, - kernel_type=KernelType.MASS_PROXY, - has_analytic_sln=False, - is_dirac_delta=False, - ) - ig_kernel2.distribution.return_value = np.atleast_1d(1.0) - empty_cluster_abundance.add_kernel(ig_kernel2) - - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) - - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 4 - assert integrator.integral_bounds == [(13, 17), (0, 2), z_limits, mass_limits] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - KernelType.Z_PROXY: 2, - KernelType.MASS_PROXY: 3, - } - - -def test_set_integration_bounds_analytic_slns( - empty_cluster_abundance: ClusterAbundance, integrator: Integrator -): - z_array = np.linspace(0, 2, 10) - m_array = np.linspace(13, 17, 10) - z_bins = list(zip(z_array[:-1], z_array[1:])) - m_bins = list(zip(m_array[:-1], m_array[1:])) - sky_area = 100**2 - - a_kernel = Mock( - spec=Kernel, - kernel_type=KernelType.MASS_PROXY, - has_analytic_sln=True, - is_dirac_delta=False, - ) - a_kernel.distribution.return_value = np.atleast_1d(1.0) - empty_cluster_abundance.add_kernel(a_kernel) - - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) - - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 2 - assert integrator.integral_bounds == [(13, 17), (0, 2)] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - } - a_kernel2 = Mock( - spec=Kernel, - kernel_type=KernelType.Z_PROXY, - has_analytic_sln=True, - is_dirac_delta=False, - ) - a_kernel2.distribution.return_value = np.atleast_1d(1.0) - empty_cluster_abundance.add_kernel(a_kernel2) - - for z_limits, mass_limits in zip(z_bins, m_bins): - integrator.set_integration_bounds( - empty_cluster_abundance, sky_area, z_limits, mass_limits - ) + def integrand(int_args: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + # xy + a = int_args[:, 0] + b = int_args[:, 1] + result = a * b + return result - assert integrator.mass_proxy_limits == mass_limits - assert integrator.z_proxy_limits == z_limits - assert len(integrator.integral_bounds) == 2 - assert integrator.integral_bounds == [(13, 17), (0, 2)] - assert integrator.integral_args_lkp == { - KernelType.MASS: 0, - KernelType.Z: 1, - } + integrator.integral_bounds = [(0, 1), (0, 1)] + integrator.extra_args = [] + result = integrator.integrate(integrand) + # \int_0^1 \int_0^1 xy dx dy = 1/4 + assert result == pytest.approx(0.25, rel=1e-15, abs=0) -def test_integrator_integrate( - empty_cluster_abundance: ClusterAbundance, integrator: Integrator -): - empty_cluster_abundance.min_mass = 0 - empty_cluster_abundance.max_mass = 1 - empty_cluster_abundance.min_z = 0 - empty_cluster_abundance.max_z = 1 - sky_area = 100**2 +def test_scipy_integrator_integrate(): + integrator = ScipyIntegrator() - def integrand( - a: npt.NDArray[np.float64], - b: npt.NDArray[np.float64], - _c: float, - _d: npt.NDArray[np.float64], - _e: npt.NDArray[np.float64], - _f: Tuple[float, float], - _g: Tuple[float, float], - ) -> npt.NDArray[np.float64]: + def integrand(a: np.float64, b: np.float64) -> np.float64: # xy result = a * b return result - integrator.set_integration_bounds( - empty_cluster_abundance, - sky_area, - (0, 1), - (0, 1), - ) + integrator.integral_bounds = [(0, 1), (0, 1)] + integrator.extra_args = [] result = integrator.integrate(integrand) # \int_0^1 \int_0^1 xy dx dy = 1/4 assert result == pytest.approx(0.25, rel=1e-15, abs=0) diff --git a/tests/test_cluster_kernels.py b/tests/test_cluster_kernels.py index 1719bdbf..749426b4 100644 --- a/tests/test_cluster_kernels.py +++ b/tests/test_cluster_kernels.py @@ -1,11 +1,10 @@ """Tests for the cluster kernel module.""" -import pytest import numpy as np +import pytest + from firecrown.models.cluster.kernel import ( Completeness, Purity, - KernelType, - Kernel, SpectroscopicRedshift, TrueMass, ) @@ -13,84 +12,41 @@ def test_create_spectroscopic_redshift_kernel(): srk = SpectroscopicRedshift() - assert isinstance(srk, Kernel) - assert srk.kernel_type == KernelType.Z_PROXY - assert srk.is_dirac_delta is True - assert srk.integral_bounds is None - assert srk.has_analytic_sln is False + assert srk is not None def test_create_mass_kernel(): mk = TrueMass() - assert isinstance(mk, Kernel) - assert mk.kernel_type == KernelType.MASS_PROXY - assert mk.is_dirac_delta is True - assert mk.integral_bounds is None - assert mk.has_analytic_sln is False + assert mk is not None def test_create_completeness_kernel(): ck = Completeness() - assert isinstance(ck, Kernel) - assert ck.kernel_type == KernelType.COMPLETENESS - assert ck.is_dirac_delta is False - assert ck.integral_bounds is None - assert ck.has_analytic_sln is False + assert ck is not None def test_create_purity_kernel(): pk = Purity() - assert isinstance(pk, Kernel) - assert pk.kernel_type == KernelType.PURITY - assert pk.is_dirac_delta is False - assert pk.integral_bounds is None - assert pk.has_analytic_sln is False + assert pk is not None def test_spec_z_distribution(): srk = SpectroscopicRedshift() - - assert ( - srk.distribution( - _mass=np.linspace(13, 17, 5), - _z=np.linspace(0, 1, 5), - _mass_proxy=np.linspace(0, 5, 5), - _z_proxy=np.linspace(0, 1, 5), - _mass_proxy_limits=(0, 5), - _z_proxy_limits=(0, 1), - ) - == 1.0 - ) + assert srk.distribution() == 1.0 def test_true_mass_distribution(): tmk = TrueMass() - - assert ( - tmk.distribution( - _mass=np.linspace(13, 17, 5), - _z=np.linspace(0, 1, 5), - _mass_proxy=np.linspace(0, 5, 5), - _z_proxy=np.linspace(0, 1, 5), - _mass_proxy_limits=(0, 5), - _z_proxy_limits=(0, 1), - ) - == 1.0 - ) + assert tmk.distribution() == 1.0 @pytest.mark.precision_sensitive def test_purity_distribution(): pk = Purity() - - mass = np.linspace(13, 17, 10) mass_proxy = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) z = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) - z_proxy = np.linspace(0, 1, 10) - mass_proxy_limits = (1.0, 10.0) - z_proxy_limits = (0.1, 1.0) truth = np.array( [ @@ -107,9 +63,7 @@ def test_purity_distribution(): ] ) - purity = pk.distribution( - mass, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits - ) + purity = pk.distribution(z, mass_proxy, mass_proxy_limits) assert isinstance(purity, np.ndarray) for ref, true in zip(purity, truth): assert ref == pytest.approx(true, rel=1e-7, abs=0.0) @@ -118,15 +72,9 @@ def test_purity_distribution(): @pytest.mark.precision_sensitive def test_purity_distribution_uses_mean(): pk = Purity() - - mass = np.linspace(13, 17, 10) - z_proxy = np.linspace(0, 1, 10) - z = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) mass_proxy = np.ones_like(z) * -1.0 - mass_proxy_limits = (1.0, 10.0) - z_proxy_limits = (0.1, 1.0) truth = np.array( [ @@ -143,9 +91,7 @@ def test_purity_distribution_uses_mean(): ] ) - purity = pk.distribution( - mass, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits - ) + purity = pk.distribution(z, mass_proxy, mass_proxy_limits) assert isinstance(purity, np.ndarray) for ref, true in zip(purity, truth): assert ref == pytest.approx(true, rel=1e-7, abs=0.0) @@ -156,11 +102,6 @@ def test_completeness_distribution(): ck = Completeness() mass = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) z = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]) - mass_proxy = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) - z_proxy = np.linspace(0, 1, 10) - - mass_proxy_limits = (1.0, 10.0) - z_proxy_limits = (0.1, 1.0) truth = np.array( [ @@ -177,9 +118,7 @@ def test_completeness_distribution(): ] ) - comp = ck.distribution( - mass, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits - ) + comp = ck.distribution(mass, z) assert isinstance(comp, np.ndarray) for ref, true in zip(comp, truth): assert ref == pytest.approx(true, rel=1e-7, abs=0.0) diff --git a/tests/test_cluster_mass_richness.py b/tests/test_cluster_mass_richness.py index 7914cdcc..0cfeb801 100644 --- a/tests/test_cluster_mass_richness.py +++ b/tests/test_cluster_mass_richness.py @@ -7,10 +7,7 @@ MurataUnbinned, MassRichnessGaussian, ) -from firecrown.models.cluster.kernel import ( - KernelType, - Kernel, -) + PIVOT_Z = 0.6 PIVOT_MASS = 14.625862906 @@ -54,11 +51,6 @@ def fixture_murata_unbinned() -> MurataUnbinned: def test_create_musigma_kernel(): mb = MurataBinned(1, 1) - assert isinstance(mb, Kernel) - assert mb.kernel_type == KernelType.MASS_PROXY - assert mb.is_dirac_delta is False - assert mb.integral_bounds is None - assert mb.has_analytic_sln is True assert mb.pivot_mass == 1 * np.log(10) assert mb.pivot_redshift == 1 assert mb.log1p_pivot_redshift == np.log1p(1) @@ -91,7 +83,6 @@ def test_cluster_observed_mass(): def test_cluster_murata_binned_distribution(murata_binned_relation: MurataBinned): mass_array = np.linspace(7.0, 26.0, 20) mass_proxy_limits = (1.0, 5.0) - z_proxy_limits = (0.0, 1.0) for z in np.geomspace(1.0e-18, 2.0, 20): flip = False @@ -99,14 +90,12 @@ def test_cluster_murata_binned_distribution(murata_binned_relation: MurataBinned mass1 = np.atleast_1d(mass1) mass2 = np.atleast_1d(mass2) z = np.atleast_1d(z) - z_proxy = np.atleast_1d(0) - mass_proxy = np.atleast_1d(1) probability_0 = murata_binned_relation.distribution( - mass1, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits + mass1, z, mass_proxy_limits ) probability_1 = murata_binned_relation.distribution( - mass2, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits + mass2, z, mass_proxy_limits ) assert probability_0 >= 0 @@ -167,8 +156,6 @@ def test_cluster_murata_binned_variance(murata_binned_relation: MurataBinned): def test_cluster_murata_unbinned_distribution(murata_unbinned_relation: MurataUnbinned): mass_array = np.linspace(7.0, 26.0, 20) - mass_proxy_limits = (1.0, 5.0) - z_proxy_limits = (0.0, 1.0) for z in np.geomspace(1.0e-18, 2.0, 20): flip = False @@ -176,15 +163,10 @@ def test_cluster_murata_unbinned_distribution(murata_unbinned_relation: MurataUn mass1 = np.atleast_1d(mass1) mass2 = np.atleast_1d(mass2) z = np.atleast_1d(z) - z_proxy = np.atleast_1d(0) mass_proxy = np.atleast_1d(1) - probability_0 = murata_unbinned_relation.distribution( - mass1, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits - ) - probability_1 = murata_unbinned_relation.distribution( - mass2, z, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits - ) + probability_0 = murata_unbinned_relation.distribution(mass1, z, mass_proxy) + probability_1 = murata_unbinned_relation.distribution(mass2, z, mass_proxy) # Probability density should be initially monotonically increasing # and then monotonically decreasing. It should flip only once. @@ -218,9 +200,7 @@ def test_cluster_murata_unbinned_distribution_is_normalized( def integrand(mass_proxy): """Evaluate the unbinned distribution at fixed mass and redshift.""" # pylint: disable=cell-var-from-loop - return murata_unbinned_relation.distribution( - mass, z, mass_proxy, z, mass_proxy_limits, (0.0, 1.0) - ) + return murata_unbinned_relation.distribution(mass, z, mass_proxy) result, _ = quad( integrand, From 5ed95ac42ba3d5a4a68971f964451ba4712bf72f Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 15:54:29 -0800 Subject: [PATCH 12/26] Fixing a ton of warnings in test from a nquad warning. --- .../statistic/test_cluster_number_counts.py | 38 +++++++++---------- tests/test_cluster_integrators.py | 1 - tests/test_cluster_mass_richness.py | 4 +- 3 files changed, 19 insertions(+), 24 deletions(-) diff --git a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py index ff34e4d8..f4cdd461 100644 --- a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py +++ b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py @@ -4,6 +4,7 @@ import pytest import pyccl from firecrown.models.cluster.integrator.integrator import Integrator +from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe from firecrown.likelihood.gauss_family.statistic.source.source import SourceSystematic from firecrown.modeling_tools import ModelingTools from firecrown.models.cluster.properties import ClusterProperty @@ -15,9 +16,8 @@ def test_create_binned_number_counts(): - integrator = Mock(spec=Integrator) - integrator.integrate.return_value = 1.0 - bnc = BinnedClusterNumberCounts(ClusterProperty.NONE, "Test", integrator) + recipe = Mock(spec=ClusterRecipe) + bnc = BinnedClusterNumberCounts(ClusterProperty.NONE, "Test", recipe) assert bnc is not None assert bnc.cluster_properties == ClusterProperty.NONE assert bnc.survey_name == "Test" @@ -26,31 +26,27 @@ def test_create_binned_number_counts(): assert len(bnc.data_vector) == 0 bnc = BinnedClusterNumberCounts( - (ClusterProperty.COUNTS | ClusterProperty.MASS), "Test", integrator + (ClusterProperty.COUNTS | ClusterProperty.MASS), "Test", recipe ) assert ClusterProperty.COUNTS in bnc.cluster_properties assert ClusterProperty.MASS in bnc.cluster_properties systematics = [SourceSystematic("mock_systematic")] - bnc = BinnedClusterNumberCounts( - ClusterProperty.NONE, "Test", integrator, systematics - ) + bnc = BinnedClusterNumberCounts(ClusterProperty.NONE, "Test", recipe, systematics) assert bnc.systematics == systematics def test_get_data_vector(): - integrator = Mock(spec=Integrator) - integrator.integrate.return_value = 1.0 - bnc = BinnedClusterNumberCounts(ClusterProperty.NONE, "Test", integrator) + recipe = Mock(spec=ClusterRecipe) + bnc = BinnedClusterNumberCounts(ClusterProperty.NONE, "Test", recipe) dv = bnc.get_data_vector() assert dv is not None assert len(dv) == 0 def test_read(cluster_sacc_data: sacc.Sacc): - integrator = Mock(spec=Integrator) - integrator.integrate.return_value = 1.0 - bnc = BinnedClusterNumberCounts(ClusterProperty.NONE, "my_survey", integrator) + recipe = Mock(spec=ClusterRecipe) + bnc = BinnedClusterNumberCounts(ClusterProperty.NONE, "my_survey", recipe) with pytest.raises( ValueError, @@ -58,7 +54,7 @@ def test_read(cluster_sacc_data: sacc.Sacc): ): bnc.read(cluster_sacc_data) - bnc = BinnedClusterNumberCounts(ClusterProperty.COUNTS, "my_survey", integrator) + bnc = BinnedClusterNumberCounts(ClusterProperty.COUNTS, "my_survey", recipe) bnc.read(cluster_sacc_data) assert bnc.sky_area == 4000 assert len(bnc.bins) == 2 @@ -66,7 +62,7 @@ def test_read(cluster_sacc_data: sacc.Sacc): assert bnc.sacc_indices is not None assert len(bnc.sacc_indices) == 2 - bnc = BinnedClusterNumberCounts(ClusterProperty.MASS, "my_survey", integrator) + bnc = BinnedClusterNumberCounts(ClusterProperty.MASS, "my_survey", recipe) bnc.read(cluster_sacc_data) assert bnc.sky_area == 4000 assert len(bnc.bins) == 2 @@ -75,7 +71,7 @@ def test_read(cluster_sacc_data: sacc.Sacc): assert len(bnc.sacc_indices) == 2 bnc = BinnedClusterNumberCounts( - (ClusterProperty.COUNTS | ClusterProperty.MASS), "my_survey", integrator + (ClusterProperty.COUNTS | ClusterProperty.MASS), "my_survey", recipe ) bnc.read(cluster_sacc_data) assert bnc.sky_area == 4000 @@ -86,8 +82,8 @@ def test_read(cluster_sacc_data: sacc.Sacc): def test_compute_theory_vector(cluster_sacc_data: sacc.Sacc): - integrator = Mock(spec=Integrator) - integrator.integrate.return_value = 1.0 + recipe = Mock(spec=ClusterRecipe) + recipe.evaluate_theory_prediction.return_value = 1.0 tools = ModelingTools() hmf = pyccl.halos.MassFuncBocquet16() @@ -98,20 +94,20 @@ def test_compute_theory_vector(cluster_sacc_data: sacc.Sacc): tools.update(params) tools.prepare(cosmo) - bnc = BinnedClusterNumberCounts(ClusterProperty.COUNTS, "my_survey", integrator) + bnc = BinnedClusterNumberCounts(ClusterProperty.COUNTS, "my_survey", recipe) bnc.read(cluster_sacc_data) tv = bnc.compute_theory_vector(tools) assert tv is not None assert len(tv) == 2 - bnc = BinnedClusterNumberCounts(ClusterProperty.MASS, "my_survey", integrator) + bnc = BinnedClusterNumberCounts(ClusterProperty.MASS, "my_survey", recipe) bnc.read(cluster_sacc_data) tv = bnc.compute_theory_vector(tools) assert tv is not None assert len(tv) == 2 bnc = BinnedClusterNumberCounts( - (ClusterProperty.COUNTS | ClusterProperty.MASS), "my_survey", integrator + (ClusterProperty.COUNTS | ClusterProperty.MASS), "my_survey", recipe ) bnc.read(cluster_sacc_data) tv = bnc.compute_theory_vector(tools) diff --git a/tests/test_cluster_integrators.py b/tests/test_cluster_integrators.py index 5ae7e512..27764e5c 100644 --- a/tests/test_cluster_integrators.py +++ b/tests/test_cluster_integrators.py @@ -3,7 +3,6 @@ import numpy as np import numpy.typing as npt import pytest -from firecrown.models.cluster.integrator.integrator import Integrator from firecrown.models.cluster.integrator.scipy_integrator import ScipyIntegrator from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator diff --git a/tests/test_cluster_mass_richness.py b/tests/test_cluster_mass_richness.py index 0cfeb801..840c99ca 100644 --- a/tests/test_cluster_mass_richness.py +++ b/tests/test_cluster_mass_richness.py @@ -197,10 +197,10 @@ def test_cluster_murata_unbinned_distribution_is_normalized( sigma = murata_unbinned_relation.get_proxy_sigma(mass, z)[0] mass_proxy_limits = (mean - 5 * sigma, mean + 5 * sigma) - def integrand(mass_proxy): + def integrand(mass_proxy) -> float: """Evaluate the unbinned distribution at fixed mass and redshift.""" # pylint: disable=cell-var-from-loop - return murata_unbinned_relation.distribution(mass, z, mass_proxy) + return murata_unbinned_relation.distribution(mass, z, mass_proxy)[0] result, _ = quad( integrand, From 8069339e65dd865ac7f98c51ed94efc8d04cbaa5 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 16:02:34 -0800 Subject: [PATCH 13/26] Fixing unused import in tests. --- .../gauss_family/statistic/test_cluster_number_counts.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py index f4cdd461..d16c3cb5 100644 --- a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py +++ b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py @@ -3,7 +3,6 @@ import sacc import pytest import pyccl -from firecrown.models.cluster.integrator.integrator import Integrator from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe from firecrown.likelihood.gauss_family.statistic.source.source import SourceSystematic from firecrown.modeling_tools import ModelingTools From 0e79927e0a131be8d91b4babc683c39e5650f769 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 16:17:05 -0800 Subject: [PATCH 14/26] Removing verbose logging to see why mypy is failling on server when all CI commands succeed locally. --- .github/workflows/ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e1a6531a..0b485e76 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -86,9 +86,9 @@ jobs: - name: Running mypy shell: bash -l {0} run: | - mypy -p firecrown -vv - mypy -p examples -vv - mypy -p tests -vv + mypy -p firecrown + mypy -p examples + mypy -p tests - name: Running pylint shell: bash -l {0} run: | From 1a9daaa4853f101798211db6905a90aaff2dbb2f Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Mon, 8 Jan 2024 16:27:57 -0800 Subject: [PATCH 15/26] Python 3.11 feature was causing 3.9 tests to fail and then cancelling the 3.11 tests. --- .github/workflows/ci.yml | 6 +++--- firecrown/models/cluster/recipes/cluster_recipe.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0b485e76..e1a6531a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -86,9 +86,9 @@ jobs: - name: Running mypy shell: bash -l {0} run: | - mypy -p firecrown - mypy -p examples - mypy -p tests + mypy -p firecrown -vv + mypy -p examples -vv + mypy -p tests -vv - name: Running pylint shell: bash -l {0} run: | diff --git a/firecrown/models/cluster/recipes/cluster_recipe.py b/firecrown/models/cluster/recipes/cluster_recipe.py index 05b606f9..329925c6 100644 --- a/firecrown/models/cluster/recipes/cluster_recipe.py +++ b/firecrown/models/cluster/recipes/cluster_recipe.py @@ -8,7 +8,7 @@ class ClusterRecipe(Updatable, ABC): - def __init__(self, parameter_prefix: str | None = None) -> None: + def __init__(self, parameter_prefix: Optional[str] = None) -> None: super().__init__(parameter_prefix) self.my_updatables: UpdatableCollection = UpdatableCollection() From 58d0f0e02fc9ca85aa6adf1fd60a0075d3ec2a09 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Tue, 9 Jan 2024 13:23:28 -0800 Subject: [PATCH 16/26] * Disabling "too-few-public-methods" in models and global rc files. --- firecrown/models/pylintrc | 1 + pylintrc | 1 + 2 files changed, 2 insertions(+) diff --git a/firecrown/models/pylintrc b/firecrown/models/pylintrc index 83f475be..502ff948 100644 --- a/firecrown/models/pylintrc +++ b/firecrown/models/pylintrc @@ -19,6 +19,7 @@ recursive=yes # multiple time (only on the command line, not in the configuration file where # it should appear only once). See also the "--disable" option for examples. enable=c-extension-no-member +disable=too-few-public-methods [MISCELLANEOUS] diff --git a/pylintrc b/pylintrc index c7fefa57..fdc3102e 100644 --- a/pylintrc +++ b/pylintrc @@ -19,6 +19,7 @@ recursive=yes # multiple time (only on the command line, not in the configuration file where # it should appear only once). See also the "--disable" option for examples. enable=c-extension-no-member +disable=too-few-public-methods [MISCELLANEOUS] From 7d7ef67d5f5e88f8f8efda243124296ee4c42f96 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Tue, 9 Jan 2024 13:24:58 -0800 Subject: [PATCH 17/26] * Adding statically typed callables to integrators (base+concr impl) * Updating the way numcosmo calls the wrapped function in the integrator for clarity * Adding statically typed callables to the cluster recipes (forgot to do this one). --- .../models/cluster/integrator/integrator.py | 15 +++++-- .../cluster/integrator/numcosmo_integrator.py | 20 ++++++--- .../cluster/integrator/scipy_integrator.py | 9 +++- .../cluster/recipes/murata_binned_spec_z.py | 45 +++++++++++++++---- .../cluster/recipes/true_mass_spec_z.py | 29 +++++++++--- tests/test_cluster_integrators.py | 13 ++++-- 6 files changed, 100 insertions(+), 31 deletions(-) diff --git a/firecrown/models/cluster/integrator/integrator.py b/firecrown/models/cluster/integrator/integrator.py index a3031125..7620623f 100644 --- a/firecrown/models/cluster/integrator/integrator.py +++ b/firecrown/models/cluster/integrator/integrator.py @@ -4,10 +4,12 @@ integrate a function. """ from abc import ABC, abstractmethod -from typing import Tuple, List, Callable +from typing import Callable, List, Tuple + +import numpy as np +import numpy.typing as npt -# pylint: disable=too-few-public-methods class Integrator(ABC): """The integrator base class @@ -15,8 +17,13 @@ class Integrator(ABC): def __init__(self) -> None: self.integral_bounds: List[Tuple[float, float]] = [] - self.extra_args: List[float] = [] + self.extra_args: npt.NDArray[np.float64] = np.array([], dtype=np.float64) @abstractmethod - def integrate(self, func_to_integrate: Callable) -> float: + def integrate( + self, + func_to_integrate: Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.NDArray[np.float64] + ], + ) -> float: """Call this method to integrate the provided integrand argument.""" diff --git a/firecrown/models/cluster/integrator/numcosmo_integrator.py b/firecrown/models/cluster/integrator/numcosmo_integrator.py index 038ade7e..8eaed1a4 100644 --- a/firecrown/models/cluster/integrator/numcosmo_integrator.py +++ b/firecrown/models/cluster/integrator/numcosmo_integrator.py @@ -3,7 +3,7 @@ This module holds the NumCosmo implementation of the integrator classes """ from enum import Enum -from typing import Callable, List, Optional, Sequence, Tuple +from typing import Callable, List, Optional, Tuple import numpy as np import numpy.typing as npt @@ -12,7 +12,6 @@ from firecrown.models.cluster.integrator.integrator import Integrator -# pylint: disable=too-few-public-methods class NumCosmoIntegralMethod(Enum): """The available NumCosmo integration methods.""" @@ -34,11 +33,16 @@ def __init__( super().__init__() self.method = method or NumCosmoIntegralMethod.P_V self.integral_bounds: List[Tuple[float, float]] = [] - self.extra_args: List[float] = [] + self.extra_args: npt.NDArray[np.float64] = np.array([], dtype=np.float64) self._relative_tolerance = relative_tolerance self._absolute_tolerance = absolute_tolerance - def integrate(self, func_to_integrate: Callable) -> float: + def integrate( + self, + func_to_integrate: Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.NDArray[np.float64] + ], + ) -> float: Ncm.cfg_init() int_nd = CountsIntegralND( @@ -61,8 +65,10 @@ class CountsIntegralND(Ncm.IntegralND): def __init__( self, dim: int, - fun: Callable[[npt.NDArray], Sequence[float]], - args: List[float], + fun: Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.NDArray[np.float64] + ], + args: npt.NDArray[np.float64], ) -> None: super().__init__() self.dim = dim @@ -85,4 +91,4 @@ def do_integrand( ) -> None: """Called by NumCosmo to evaluate the integrand.""" x = np.array(x_vec.dup_array()).reshape(npoints, dim) - fval_vec.set_array(self.fun(x, *self.extra_args)) + fval_vec.set_array(list(self.fun(x, self.extra_args))) diff --git a/firecrown/models/cluster/integrator/scipy_integrator.py b/firecrown/models/cluster/integrator/scipy_integrator.py index 40befb24..4869f60c 100644 --- a/firecrown/models/cluster/integrator/scipy_integrator.py +++ b/firecrown/models/cluster/integrator/scipy_integrator.py @@ -3,11 +3,14 @@ This module holds the scipy implementation of the integrator classes """ from typing import Callable + +import numpy as np +import numpy.typing as npt from scipy.integrate import nquad + from firecrown.models.cluster.integrator.integrator import Integrator -# pylint: disable=too-few-public-methods class ScipyIntegrator(Integrator): """The scipy implementation of the Integrator base class using nquad.""" @@ -20,7 +23,9 @@ def __init__( def integrate( self, - func_to_integrate: Callable, + func_to_integrate: Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.NDArray[np.float64] + ], ) -> float: val = nquad( func_to_integrate, diff --git a/firecrown/models/cluster/recipes/murata_binned_spec_z.py b/firecrown/models/cluster/recipes/murata_binned_spec_z.py index 881a09ce..12d1a6ca 100644 --- a/firecrown/models/cluster/recipes/murata_binned_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_binned_spec_z.py @@ -1,4 +1,7 @@ -from typing import Callable, Optional +from typing import Callable, Optional, Tuple + +import numpy as np +import numpy.typing as npt from firecrown.models.cluster.abundance import ClusterAbundance from firecrown.models.cluster.binning import NDimensionalBin @@ -23,10 +26,18 @@ def get_theory_prediction( self, cluster_theory: ClusterAbundance, average_on: Optional[ClusterProperty] = None, - ) -> Callable: + ) -> Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64], Tuple[float, float], float], + npt.NDArray[np.float64], + ]: """_summary_""" - def theory_prediction(mass, z, mass_proxy_limits, sky_area): + def theory_prediction( + mass: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + mass_proxy_limits: Tuple[float, float], + sky_area: float, + ): prediction = ( cluster_theory.comoving_volume(z, sky_area) * cluster_theory.mass_function(mass, z) @@ -52,15 +63,31 @@ def theory_prediction(mass, z, mass_proxy_limits, sky_area): return theory_prediction - def get_function_to_integrate(self, prediction: Callable) -> Callable: - def numcosmo_wrapper( - int_args, mass_proxy_low, mass_proxy_high, sky_area - ) -> float: + def get_function_to_integrate( + self, + prediction: Callable[ + [ + npt.NDArray[np.float64], + npt.NDArray[np.float64], + Tuple[float, float], + float, + ], + npt.NDArray[np.float64], + ], + ) -> Callable[[npt.NDArray, npt.NDArray], npt.NDArray]: + def function_mapper( + int_args: npt.NDArray, extra_args: npt.NDArray + ) -> npt.NDArray[np.float64]: mass = int_args[:, 0] z = int_args[:, 1] + + mass_proxy_low = extra_args[0] + mass_proxy_high = extra_args[1] + sky_area = extra_args[2] + return prediction(mass, z, (mass_proxy_low, mass_proxy_high), sky_area) - return numcosmo_wrapper + return function_mapper def evaluate_theory_prediction( self, @@ -73,7 +100,7 @@ def evaluate_theory_prediction( (cluster_theory.min_mass, cluster_theory.max_mass), bin.z_edges, ] - self.integrator.extra_args = [*bin.mass_proxy_edges, sky_area] + self.integrator.extra_args = np.array([*bin.mass_proxy_edges, sky_area]) theory_prediction = self.get_theory_prediction(cluster_theory, average_on) prediction_wrapper = self.get_function_to_integrate(theory_prediction) diff --git a/firecrown/models/cluster/recipes/true_mass_spec_z.py b/firecrown/models/cluster/recipes/true_mass_spec_z.py index 35fd8714..d4653694 100644 --- a/firecrown/models/cluster/recipes/true_mass_spec_z.py +++ b/firecrown/models/cluster/recipes/true_mass_spec_z.py @@ -1,5 +1,8 @@ from typing import Callable, Optional +import numpy as np +import numpy.typing as npt + from firecrown.models.cluster.abundance import ClusterAbundance from firecrown.models.cluster.binning import NDimensionalBin from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator @@ -20,10 +23,15 @@ def get_theory_prediction( self, cluster_theory: ClusterAbundance, average_on: Optional[ClusterProperty] = None, - ) -> Callable: + ) -> Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64], float], + npt.NDArray[np.float64], + ]: """_summary_""" - def theory_prediction(mass, z, sky_area): + def theory_prediction( + mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], sky_area: float + ): prediction = ( cluster_theory.comoving_volume(z, sky_area) * cluster_theory.mass_function(mass, z) @@ -49,10 +57,21 @@ def theory_prediction(mass, z, sky_area): return theory_prediction - def get_function_to_integrate(self, prediction: Callable) -> Callable: - def numcosmo_wrapper(int_args, sky_area) -> float: + def get_function_to_integrate( + self, + prediction: Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64], float], + npt.NDArray[np.float64], + ], + ) -> Callable[ + [npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.NDArray[np.float64] + ]: + def numcosmo_wrapper( + int_args: npt.NDArray[np.float64], extra_args: npt.NDArray[np.float64] + ) -> npt.NDArray[np.float64]: mass = int_args[:, 0] z = int_args[:, 1] + sky_area = extra_args[0] return prediction(mass, z, sky_area) return numcosmo_wrapper @@ -72,7 +91,7 @@ def evaluate_theory_prediction( bin.z_edges, ] print(mass_low, mass_high) - self.integrator.extra_args = [sky_area] + self.integrator.extra_args = np.array([sky_area], dtype=np.float64) theory_prediction = self.get_theory_prediction(cluster_theory, average_on) prediction_wrapper = self.get_function_to_integrate(theory_prediction) diff --git a/tests/test_cluster_integrators.py b/tests/test_cluster_integrators.py index 27764e5c..fc1d6b93 100644 --- a/tests/test_cluster_integrators.py +++ b/tests/test_cluster_integrators.py @@ -15,7 +15,9 @@ def test_numcosmo_integrator_integrate(): integrator = NumCosmoIntegrator() - def integrand(int_args: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + def integrand( + int_args: npt.NDArray[np.float64], _extra_args: npt.NDArray[np.float64] + ) -> npt.NDArray[np.float64]: # xy a = int_args[:, 0] b = int_args[:, 1] @@ -23,7 +25,7 @@ def integrand(int_args: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: return result integrator.integral_bounds = [(0, 1), (0, 1)] - integrator.extra_args = [] + integrator.extra_args = np.array([], dtype=np.float64) result = integrator.integrate(integrand) # \int_0^1 \int_0^1 xy dx dy = 1/4 assert result == pytest.approx(0.25, rel=1e-15, abs=0) @@ -32,13 +34,16 @@ def integrand(int_args: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: def test_scipy_integrator_integrate(): integrator = ScipyIntegrator() - def integrand(a: np.float64, b: np.float64) -> np.float64: + # TODO: should we just remove this? + def integrand( + a: npt.NDArray[np.float64], b: npt.NDArray[np.float64] + ) -> npt.NDArray[np.float64]: # xy result = a * b return result integrator.integral_bounds = [(0, 1), (0, 1)] - integrator.extra_args = [] + integrator.extra_args = np.array([], dtype=np.float64) result = integrator.integrate(integrand) # \int_0^1 \int_0^1 xy dx dy = 1/4 assert result == pytest.approx(0.25, rel=1e-15, abs=0) From aa014be0e6c038d74db967234717d1a1b2ca1fbb Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Tue, 9 Jan 2024 14:07:18 -0800 Subject: [PATCH 18/26] * Removing bad logic in __eq__ that was causing a bug **FOUND using unit tests** * Adding TupleBin tests --- firecrown/models/cluster/binning.py | 3 +- .../models/cluster/recipes/cluster_recipe.py | 2 +- .../statistic/test_cluster_number_counts.py | 9 +- tests/test_cluster_binning.py | 89 ++++++++++++++++++- 4 files changed, 98 insertions(+), 5 deletions(-) diff --git a/firecrown/models/cluster/binning.py b/firecrown/models/cluster/binning.py index 6bb9a99b..272a060f 100644 --- a/firecrown/models/cluster/binning.py +++ b/firecrown/models/cluster/binning.py @@ -71,8 +71,7 @@ def __eq__(self, other: object) -> bool: other_bin = [ x for x in other.coordinate_bins if x.tracer_type == my_bin.tracer_type ] - if len(other_bin) != 1: - return False + if my_bin.lower != other_bin[0].lower: return False if my_bin.upper != other_bin[0].upper: diff --git a/firecrown/models/cluster/recipes/cluster_recipe.py b/firecrown/models/cluster/recipes/cluster_recipe.py index 329925c6..fd9b0fa9 100644 --- a/firecrown/models/cluster/recipes/cluster_recipe.py +++ b/firecrown/models/cluster/recipes/cluster_recipe.py @@ -20,4 +20,4 @@ def evaluate_theory_prediction( sky_area: float, average_on: Optional[ClusterProperty] = None, ) -> float: - pass + """Evaluate the theory prediction for this cluster recipe""" diff --git a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py index d16c3cb5..cf542d14 100644 --- a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py +++ b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py @@ -43,7 +43,7 @@ def test_get_data_vector(): assert len(dv) == 0 -def test_read(cluster_sacc_data: sacc.Sacc): +def test_read_throws_if_no_property(cluster_sacc_data: sacc.Sacc): recipe = Mock(spec=ClusterRecipe) bnc = BinnedClusterNumberCounts(ClusterProperty.NONE, "my_survey", recipe) @@ -53,6 +53,10 @@ def test_read(cluster_sacc_data: sacc.Sacc): ): bnc.read(cluster_sacc_data) + +def test_read_single_property(cluster_sacc_data: sacc.Sacc): + recipe = Mock(spec=ClusterRecipe) + bnc = BinnedClusterNumberCounts(ClusterProperty.COUNTS, "my_survey", recipe) bnc.read(cluster_sacc_data) assert bnc.sky_area == 4000 @@ -69,6 +73,9 @@ def test_read(cluster_sacc_data: sacc.Sacc): assert bnc.sacc_indices is not None assert len(bnc.sacc_indices) == 2 + +def test_read_multiple_properties(cluster_sacc_data: sacc.Sacc): + recipe = Mock(spec=ClusterRecipe) bnc = BinnedClusterNumberCounts( (ClusterProperty.COUNTS | ClusterProperty.MASS), "my_survey", recipe ) diff --git a/tests/test_cluster_binning.py b/tests/test_cluster_binning.py index 2304e1a6..4684ec26 100644 --- a/tests/test_cluster_binning.py +++ b/tests/test_cluster_binning.py @@ -1,7 +1,26 @@ """Tests for the cluster binning module""" from unittest.mock import Mock + +import pytest import sacc -from firecrown.models.cluster.binning import SaccBin, NDimensionalBin + +from firecrown.models.cluster.binning import NDimensionalBin, SaccBin, TupleBin + + +def test_bin_repr(): + tracer_z = sacc.tracers.BinZTracer("", 0, 1) + tracer_lambda = sacc.tracers.BinRichnessTracer("", 4, 5) + + sb = SaccBin([tracer_z, tracer_lambda]) + assert repr(sb) == "[(0, 1), (4, 5)]\n" + + +def test_bin_str(): + tracer_z = sacc.tracers.BinZTracer("", 0, 1) + tracer_lambda = sacc.tracers.BinRichnessTracer("", 4, 5) + + sb = SaccBin([tracer_z, tracer_lambda]) + assert str(sb) == "[(0, 1), (4, 5)]\n" def test_create_sacc_bin_with_correct_dimension(): @@ -25,6 +44,22 @@ def test_sacc_bin_z_edges(): assert sb.z_edges == (0, 1) +def test_sacc_bin_z_edges_throws_when_multiple_z_bins(): + tracer_z = sacc.tracers.BinZTracer("", 0, 1) + tracer_lambda = sacc.tracers.BinRichnessTracer("", 4, 5) + sb = SaccBin([tracer_z, tracer_z, tracer_lambda]) + with pytest.raises(ValueError, match="SaccBin must have exactly one z bin"): + print(sb.z_edges) + + +def test_sacc_bin_mass_proxy_edges_throws_when_multiple_mass_proxy_bins(): + tracer_z = sacc.tracers.BinZTracer("", 0, 1) + tracer_lambda = sacc.tracers.BinRichnessTracer("", 4, 5) + sb = SaccBin([tracer_z, tracer_lambda, tracer_lambda]) + with pytest.raises(ValueError, match="SaccBin must have exactly one richness bin"): + print(sb.mass_proxy_edges) + + def test_sacc_bin_richness_edges(): tracer_z = sacc.tracers.BinZTracer("", 0, 1) tracer_lambda = sacc.tracers.BinRichnessTracer("", 4, 5) @@ -95,3 +130,55 @@ def test_sacc_bin_must_be_equal_type(): sb = SaccBin([tracer_z, tracer_lambda]) assert sb != other_bin + + +def test_create_tuple_bin(): + tb = TupleBin([(1, 2), (3, 4)]) + assert tb is not None + assert tb.dimension == 2 + assert tb.mass_proxy_edges == (1, 2) + assert tb.z_edges == (3, 4) + + +def test_tuple_bins_are_equal(): + tb1 = TupleBin([(1, 2), (3, 4)]) + tb2 = TupleBin([(1, 2), (3, 4)]) + assert tb1 == tb2 + assert hash(tb1) == hash(tb2) + assert tb1 is not tb2 + + +def test_tuple_bin_neq_sacc_bin(): + tb = TupleBin([(1, 2), (3, 4)]) + tracer_z = sacc.tracers.BinZTracer("", 3, 4) + tracer_lambda = sacc.tracers.BinRichnessTracer("", 1, 2) + sb = SaccBin([tracer_z, tracer_lambda]) + assert tb != sb + + +def test_tuple_bin_different_dimensions_not_equal(): + tb1 = TupleBin([(1, 2), (3, 4), (5, 6)]) + tb2 = TupleBin([(1, 2), (3, 4)]) + assert tb1 != tb2 + + tb2 = TupleBin([(1, 2, 3), (3, 4)]) + tb1 = TupleBin([(1, 2), (3, 4)]) + assert tb1 != tb2 + + +def test_tuple_bin_different_bins_not_equal(): + tb1 = TupleBin([(1, 2), (3, 4)]) + tb2 = TupleBin([(1, 2), (3, 5)]) + assert tb1 != tb2 + + tb1 = TupleBin([(0, 2), (3, 4)]) + tb2 = TupleBin([(1, 2), (3, 4)]) + assert tb1 != tb2 + + tb1 = TupleBin([(0, 2), (0, 4)]) + tb2 = TupleBin([(1, 2), (3, 4)]) + assert tb1 != tb2 + + tb1 = TupleBin([(3, 4), (1, 2)]) + tb2 = TupleBin([(1, 2), (3, 4)]) + assert tb1 != tb2 From 4e8ecd8b2fbdfeb4ccc355c9291d3404889904a1 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Tue, 9 Jan 2024 14:54:28 -0800 Subject: [PATCH 19/26] * Removing unused fixture from conftest * Adding an init for cluster recipes * Writing (really, moving) tests for the cluster recipes --- firecrown/models/cluster/recipes/__init__.py | 0 .../test_murata_binned_spec_z.py | 184 ++++++++++++++++++ tests/conftest.py | 13 -- tests/test_cluster_abundance.py | 96 +-------- 4 files changed, 187 insertions(+), 106 deletions(-) create mode 100644 firecrown/models/cluster/recipes/__init__.py create mode 100644 tests/cluster_recipes/test_murata_binned_spec_z.py diff --git a/firecrown/models/cluster/recipes/__init__.py b/firecrown/models/cluster/recipes/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/cluster_recipes/test_murata_binned_spec_z.py b/tests/cluster_recipes/test_murata_binned_spec_z.py new file mode 100644 index 00000000..22ebb6d5 --- /dev/null +++ b/tests/cluster_recipes/test_murata_binned_spec_z.py @@ -0,0 +1,184 @@ +"""Tests for the cluster abundance module.""" +from unittest.mock import Mock + +import numpy as np +import pyccl +import pytest + +from firecrown.models.cluster.abundance import ClusterAbundance +from firecrown.models.cluster.binning import NDimensionalBin +from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator +from firecrown.models.cluster.kernel import SpectroscopicRedshift +from firecrown.models.cluster.mass_proxy import MurataBinned +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe +from firecrown.models.cluster.recipes.murata_binned_spec_z import ( + MurataBinnedSpecZRecipe, +) + + +@pytest.fixture(name="cluster_abundance") +def fixture_cluster_abundance() -> ClusterAbundance: + hmf = pyccl.halos.MassFuncBocquet16() + cl_abundance = ClusterAbundance( + min_z=0, + max_z=2, + min_mass=13, + max_mass=17, + halo_mass_function=hmf, + ) + cl_abundance.update_ingredients(pyccl.CosmologyVanillaLCDM()) + return cl_abundance + + +@pytest.fixture(name="murata_binned_spec_z") +def fixture_murata_binned_spec_z() -> MurataBinnedSpecZRecipe: + cluster_recipe = MurataBinnedSpecZRecipe() + cluster_recipe.mass_distribution.mu_p0 = 3.0 + cluster_recipe.mass_distribution.mu_p1 = 0.86 + cluster_recipe.mass_distribution.mu_p2 = 0.0 + cluster_recipe.mass_distribution.sigma_p0 = 3.0 + cluster_recipe.mass_distribution.sigma_p1 = 0.7 + cluster_recipe.mass_distribution.sigma_p2 = 0.0 + return cluster_recipe + + +def test_murata_binned_spec_z_init(): + recipe = MurataBinnedSpecZRecipe() + + assert recipe is not None + assert isinstance(recipe, ClusterRecipe) + assert recipe.integrator is not None + assert isinstance(recipe.integrator, NumCosmoIntegrator) + assert recipe.redshift_distribution is not None + assert isinstance(recipe.redshift_distribution, SpectroscopicRedshift) + assert recipe.mass_distribution is not None + assert isinstance(recipe.mass_distribution, MurataBinned) + assert recipe.my_updatables is not None + assert len(recipe.my_updatables) == 1 + assert recipe.my_updatables[0] is recipe.mass_distribution + + +def test_get_theory_prediction_returns_value( + cluster_abundance: ClusterAbundance, + murata_binned_spec_z: MurataBinnedSpecZRecipe, +): + prediction = murata_binned_spec_z.get_theory_prediction(cluster_abundance) + + assert prediction is not None + assert callable(prediction) + + mass = np.linspace(13, 17, 2) + z = np.linspace(0.1, 1, 2) + mass_proxy_limits = (0, 5) + sky_area = 360**2 + + result = prediction(mass, z, mass_proxy_limits, sky_area) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + +def test_get_theory_prediction_with_average_returns_value( + cluster_abundance: ClusterAbundance, + murata_binned_spec_z: MurataBinnedSpecZRecipe, +): + mass = np.linspace(13, 17, 2) + z = np.linspace(0.1, 1, 2) + mass_proxy_limits = (0, 5) + sky_area = 360**2 + + prediction = murata_binned_spec_z.get_theory_prediction( + cluster_abundance, average_on=ClusterProperty.MASS + ) + + assert prediction is not None + assert callable(prediction) + + result = prediction(mass, z, mass_proxy_limits, sky_area) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + prediction = murata_binned_spec_z.get_theory_prediction( + cluster_abundance, average_on=ClusterProperty.REDSHIFT + ) + + assert prediction is not None + assert callable(prediction) + + result = prediction(mass, z, mass_proxy_limits, sky_area) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + prediction = murata_binned_spec_z.get_theory_prediction( + cluster_abundance, average_on=(ClusterProperty.REDSHIFT | ClusterProperty.MASS) + ) + + assert prediction is not None + assert callable(prediction) + + result = prediction(mass, z, mass_proxy_limits, sky_area) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + +def test_get_theory_prediction_throws_with_nonimpl_average( + cluster_abundance: ClusterAbundance, + murata_binned_spec_z: MurataBinnedSpecZRecipe, +): + prediction = murata_binned_spec_z.get_theory_prediction( + cluster_abundance, average_on=ClusterProperty.SHEAR + ) + + assert prediction is not None + assert callable(prediction) + + mass = np.linspace(13, 17, 2) + z = np.linspace(0.1, 1, 2) + mass_proxy_limits = (0, 5) + sky_area = 360**2 + + with pytest.raises(NotImplementedError): + _ = prediction(mass, z, mass_proxy_limits, sky_area) + + +def test_get_function_to_integrate_returns_value( + cluster_abundance: ClusterAbundance, + murata_binned_spec_z: MurataBinnedSpecZRecipe, +): + prediction = murata_binned_spec_z.get_theory_prediction(cluster_abundance) + function_to_integrate = murata_binned_spec_z.get_function_to_integrate(prediction) + + assert function_to_integrate is not None + assert callable(function_to_integrate) + + int_args = np.array([[13.0, 0.1], [17.0, 1.0]]) + extra_args = np.array([0, 5, 360**2]) + + result = function_to_integrate(int_args, extra_args) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + +def test_evaluates_theory_prediction_returns_value( + cluster_abundance: ClusterAbundance, + murata_binned_spec_z: MurataBinnedSpecZRecipe, +): + mock_bin = Mock(spec=NDimensionalBin) + mock_bin.mass_proxy_edges = (0, 5) + mock_bin.z_edges = (0, 1) + + prediction = murata_binned_spec_z.evaluate_theory_prediction( + cluster_abundance, mock_bin, 360**2 + ) + + assert prediction > 0 diff --git a/tests/conftest.py b/tests/conftest.py index ea5b5c33..6944986e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -15,7 +15,6 @@ from firecrown.parameters import ParamsMap from firecrown.connector.mapping import MappingCosmoSIS, mapping_builder from firecrown.modeling_tools import ModelingTools -from firecrown.models.cluster.abundance import ClusterAbundance def pytest_addoption(parser): @@ -157,15 +156,3 @@ def fixture_cluster_sacc_data() -> sacc.Sacc: s.add_data_point(mlm, ("not_my_survey", "z_bin_tracer_2", "mass_bin_tracer_2"), 1) return s - - -@pytest.fixture(name="empty_cluster_abundance") -def fixture_empty_cluster_abundance() -> ClusterAbundance: - cl_abundance = ClusterAbundance( - min_z=0, - max_z=2, - min_mass=13, - max_mass=17, - halo_mass_function=None, - ) - return cl_abundance diff --git a/tests/test_cluster_abundance.py b/tests/test_cluster_abundance.py index 24be2121..7dd17178 100644 --- a/tests/test_cluster_abundance.py +++ b/tests/test_cluster_abundance.py @@ -1,7 +1,8 @@ """Tests for the cluster abundance module.""" -import pytest -import pyccl import numpy as np +import pyccl +import pytest + from firecrown.models.cluster.abundance import ClusterAbundance @@ -13,18 +14,6 @@ def fixture_cluster_abundance(): return ca -@pytest.fixture(name="integrand_args") -def fixture_integrand_args(): - sky_area = 439.78986 - mass = np.linspace(13, 17, 5) - z = np.linspace(0, 1, 5) - mass_proxy = np.linspace(0, 5, 5) - z_proxy = np.linspace(0, 1, 5) - mass_proxy_limits = (0, 5) - z_proxy_limits = (0, 1) - return (mass, z, sky_area, mass_proxy, z_proxy, mass_proxy_limits, z_proxy_limits) - - def test_cluster_abundance_init(cluster_abundance: ClusterAbundance): assert cluster_abundance is not None assert cluster_abundance.cosmo is None @@ -72,82 +61,3 @@ def test_abundance_massfunc_returns_value(cluster_abundance: ClusterAbundance): assert np.issubdtype(result.dtype, np.float64) assert len(result) == 5 assert np.all(result > 0) - - -# @pytest.mark.slow -# def test_abundance_get_integrand( -# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -# ): -# cosmo = pyccl.CosmologyVanillaLCDM() -# cluster_abundance.update_ingredients(cosmo) -# cluster_abundance.add_kernel(integrable_kernel) - -# integrand = cluster_abundance.get_integrand() -# assert callable(integrand) -# result = integrand(*integrand_args) -# assert isinstance(result, np.ndarray) -# assert np.issubdtype(result.dtype, np.float64) - - -# @pytest.mark.slow -# def test_abundance_get_integrand_avg_mass( -# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -# ): -# cosmo = pyccl.CosmologyVanillaLCDM() -# cluster_abundance.update_ingredients(cosmo) -# cluster_abundance.add_kernel(integrable_kernel) - -# integrand = cluster_abundance.get_integrand(average_properties= -# ClusterProperty.MASS) -# assert callable(integrand) -# result = integrand(*integrand_args) -# assert isinstance(result, np.ndarray) -# assert np.issubdtype(result.dtype, np.float64) - - -# @pytest.mark.slow -# def test_abundance_get_integrand_avg_redshift( -# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -# ): -# cosmo = pyccl.CosmologyVanillaLCDM() -# cluster_abundance.update_ingredients(cosmo) -# cluster_abundance.add_kernel(integrable_kernel) - -# integrand = cluster_abundance.get_integrand( -# average_properties=ClusterProperty.REDSHIFT -# ) -# assert callable(integrand) -# result = integrand(*integrand_args) -# assert isinstance(result, np.ndarray) -# assert np.issubdtype(result.dtype, np.float64) - - -# @pytest.mark.slow -# def test_abundance_get_integrand_avg_mass_and_redshift( -# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -# ): -# cosmo = pyccl.CosmologyVanillaLCDM() -# cluster_abundance.update_ingredients(cosmo) -# cluster_abundance.add_kernel(integrable_kernel) - -# average_properties = ClusterProperty.MASS | ClusterProperty.REDSHIFT -# integrand = cluster_abundance.get_integrand(average_properties=average_properties) -# assert callable(integrand) -# result = integrand(*integrand_args) -# assert isinstance(result, np.ndarray) -# assert np.issubdtype(result.dtype, np.float64) - - -# @pytest.mark.slow -# def test_abundance_get_integrand_avg_not_implemented_throws( -# cluster_abundance: ClusterAbundance, integrable_kernel: Kernel, integrand_args -# ): -# cosmo = pyccl.CosmologyVanillaLCDM() -# cluster_abundance.update_ingredients(cosmo) -# cluster_abundance.add_kernel(integrable_kernel) - -# average_properties = ClusterProperty.SHEAR -# integrand = cluster_abundance.get_integrand(average_properties=average_properties) -# assert callable(integrand) -# with pytest.raises(NotImplementedError): -# _ = integrand(*integrand_args) From 6716a99acea841312f0829f31a48a1efb7414cf3 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Tue, 9 Jan 2024 15:21:11 -0800 Subject: [PATCH 20/26] * Removing unused code * Adding docstrings * Adding an unbinned cluster recipe as proof of concept. --- .../cluster_redshift_richness.py | 3 - firecrown/models/cluster/abundance.py | 16 +----- .../cluster/integrator/numcosmo_integrator.py | 4 +- firecrown/models/cluster/mass_proxy.py | 2 - .../models/cluster/recipes/cluster_recipe.py | 8 ++- .../cluster/recipes/murata_binned_spec_z.py | 16 ++++-- ...ss_spec_z.py => murata_unbinned_spec_z.py} | 57 +++++++++++++------ 7 files changed, 60 insertions(+), 46 deletions(-) rename firecrown/models/cluster/recipes/{true_mass_spec_z.py => murata_unbinned_spec_z.py} (61%) diff --git a/examples/cluster_number_counts/cluster_redshift_richness.py b/examples/cluster_number_counts/cluster_redshift_richness.py index a259529c..2b2c5e02 100644 --- a/examples/cluster_number_counts/cluster_redshift_richness.py +++ b/examples/cluster_number_counts/cluster_redshift_richness.py @@ -18,8 +18,6 @@ MurataBinnedSpecZRecipe, ) -# from firecrown.models.cluster.recipes.true_mass_spec_z import TrueMassSpecZRecipe - def get_cluster_abundance() -> ClusterAbundance: hmf = ccl.halos.MassFuncBocquet16() @@ -47,7 +45,6 @@ def build_likelihood( survey_name = "numcosmo_simulated_redshift_richness" likelihood = ConstGaussian( [BinnedClusterNumberCounts(average_on, survey_name, MurataBinnedSpecZRecipe())] - # [BinnedClusterNumberCounts(average_on, survey_name, TrueMassSpecZRecipe())] ) # Read in sacc data diff --git a/firecrown/models/cluster/abundance.py b/firecrown/models/cluster/abundance.py index dd97bbe2..8c9b9a51 100644 --- a/firecrown/models/cluster/abundance.py +++ b/firecrown/models/cluster/abundance.py @@ -4,7 +4,7 @@ and phenomenological predictions. This module contains the classes and functions that produce those predictions. """ -from typing import Callable, Dict, Tuple +from typing import Dict, Tuple import numpy as np import numpy.typing as npt import pyccl @@ -13,20 +13,6 @@ from firecrown.updatable import Updatable, UpdatableCollection -AbundanceIntegrand = Callable[ - [ - npt.NDArray[np.float64], - npt.NDArray[np.float64], - float, - npt.NDArray[np.float64], - npt.NDArray[np.float64], - Tuple[float, float], - Tuple[float, float], - ], - npt.NDArray[np.float64], -] - - class ClusterAbundance(Updatable): """The class that calculates the predicted number counts of galaxy clusters diff --git a/firecrown/models/cluster/integrator/numcosmo_integrator.py b/firecrown/models/cluster/integrator/numcosmo_integrator.py index 8eaed1a4..e457620b 100644 --- a/firecrown/models/cluster/integrator/numcosmo_integrator.py +++ b/firecrown/models/cluster/integrator/numcosmo_integrator.py @@ -3,7 +3,7 @@ This module holds the NumCosmo implementation of the integrator classes """ from enum import Enum -from typing import Callable, List, Optional, Tuple +from typing import Callable, Optional, Tuple import numpy as np import numpy.typing as npt @@ -32,8 +32,6 @@ def __init__( ) -> None: super().__init__() self.method = method or NumCosmoIntegralMethod.P_V - self.integral_bounds: List[Tuple[float, float]] = [] - self.extra_args: npt.NDArray[np.float64] = np.array([], dtype=np.float64) self._relative_tolerance = relative_tolerance self._absolute_tolerance = absolute_tolerance diff --git a/firecrown/models/cluster/mass_proxy.py b/firecrown/models/cluster/mass_proxy.py index 0125c83f..54170f7d 100644 --- a/firecrown/models/cluster/mass_proxy.py +++ b/firecrown/models/cluster/mass_proxy.py @@ -180,8 +180,6 @@ def __init__( self.sigma_p1 = parameters.register_new_updatable_parameter() self.sigma_p2 = parameters.register_new_updatable_parameter() - # Verify this gets called last or first - def get_proxy_mean( self, mass: npt.NDArray[np.float64], diff --git a/firecrown/models/cluster/recipes/cluster_recipe.py b/firecrown/models/cluster/recipes/cluster_recipe.py index fd9b0fa9..4de2bf49 100644 --- a/firecrown/models/cluster/recipes/cluster_recipe.py +++ b/firecrown/models/cluster/recipes/cluster_recipe.py @@ -1,3 +1,4 @@ +"""Module for defining the ClusterRecipe class""" from abc import ABC, abstractmethod from typing import Optional @@ -8,6 +9,11 @@ class ClusterRecipe(Updatable, ABC): + """Abstract class defining a cluster recipe. + + A cluster recipe is a combination of different cluster theoretrical predictions + and models that produces a single prediction for an observable.""" + def __init__(self, parameter_prefix: Optional[str] = None) -> None: super().__init__(parameter_prefix) self.my_updatables: UpdatableCollection = UpdatableCollection() @@ -16,7 +22,7 @@ def __init__(self, parameter_prefix: Optional[str] = None) -> None: def evaluate_theory_prediction( self, cluster_theory: ClusterAbundance, - bin: SaccBin, + this_bin: SaccBin, sky_area: float, average_on: Optional[ClusterProperty] = None, ) -> float: diff --git a/firecrown/models/cluster/recipes/murata_binned_spec_z.py b/firecrown/models/cluster/recipes/murata_binned_spec_z.py index 12d1a6ca..9573c867 100644 --- a/firecrown/models/cluster/recipes/murata_binned_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_binned_spec_z.py @@ -1,3 +1,4 @@ +"""Module for defining the classes used in the MurataBinnedSpecZ cluster recipe.""" from typing import Callable, Optional, Tuple import numpy as np @@ -12,7 +13,11 @@ from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe +# pylint: disable=R0801 class MurataBinnedSpecZRecipe(ClusterRecipe): + """Cluster recipe using the Murata 2019 binned mass-richness relation and assuming + perfectly measured spec-zs.""" + def __init__(self) -> None: super().__init__() @@ -30,8 +35,6 @@ def get_theory_prediction( [npt.NDArray[np.float64], npt.NDArray[np.float64], Tuple[float, float], float], npt.NDArray[np.float64], ]: - """_summary_""" - def theory_prediction( mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], @@ -92,15 +95,18 @@ def function_mapper( def evaluate_theory_prediction( self, cluster_theory: ClusterAbundance, - bin: NDimensionalBin, + this_bin: NDimensionalBin, sky_area: float, average_on: Optional[ClusterProperty] = None, ) -> float: + """Evaluate the theoretical prediction for the observable in the provided bin + using the Murata 2019 binned mass-richness relation and assuming perfectly + measured redshifts.""" self.integrator.integral_bounds = [ (cluster_theory.min_mass, cluster_theory.max_mass), - bin.z_edges, + this_bin.z_edges, ] - self.integrator.extra_args = np.array([*bin.mass_proxy_edges, sky_area]) + self.integrator.extra_args = np.array([*this_bin.mass_proxy_edges, sky_area]) theory_prediction = self.get_theory_prediction(cluster_theory, average_on) prediction_wrapper = self.get_function_to_integrate(theory_prediction) diff --git a/firecrown/models/cluster/recipes/true_mass_spec_z.py b/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py similarity index 61% rename from firecrown/models/cluster/recipes/true_mass_spec_z.py rename to firecrown/models/cluster/recipes/murata_unbinned_spec_z.py index d4653694..fff0e3a4 100644 --- a/firecrown/models/cluster/recipes/true_mass_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py @@ -1,3 +1,4 @@ +"""Module for defining the classes used in the MurataUnbinnedSpecZ cluster recipe.""" from typing import Callable, Optional import numpy as np @@ -6,37 +7,53 @@ from firecrown.models.cluster.abundance import ClusterAbundance from firecrown.models.cluster.binning import NDimensionalBin from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator -from firecrown.models.cluster.kernel import SpectroscopicRedshift, TrueMass +from firecrown.models.cluster.kernel import SpectroscopicRedshift +from firecrown.models.cluster.mass_proxy import MurataUnbinned from firecrown.models.cluster.properties import ClusterProperty from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe -class TrueMassSpecZRecipe(ClusterRecipe): +# pylint: disable=R0801 +class MurataUnbinnedSpecZ(ClusterRecipe): + """Cluster recipe using the Murata 2019 unbinned mass-richness relation and assuming + perfectly measured spec-zs.""" + def __init__(self) -> None: super().__init__() self.integrator = NumCosmoIntegrator() self.redshift_distribution = SpectroscopicRedshift() - self.mass_distribution = TrueMass() + pivot_mass, pivot_redshift = 14.625862906, 0.6 + self.mass_distribution = MurataUnbinned(pivot_mass, pivot_redshift) + self.my_updatables.append(self.mass_distribution) def get_theory_prediction( self, cluster_theory: ClusterAbundance, average_on: Optional[ClusterProperty] = None, ) -> Callable[ - [npt.NDArray[np.float64], npt.NDArray[np.float64], float], + [ + npt.NDArray[np.float64], + npt.NDArray[np.float64], + npt.NDArray[np.float64], + float, + ], npt.NDArray[np.float64], ]: - """_summary_""" + import pdb def theory_prediction( - mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], sky_area: float + mass: npt.NDArray[np.float64], + z: npt.NDArray[np.float64], + mass_proxy: npt.NDArray[np.float64], + sky_area: float, ): + pdb.set_trace() prediction = ( cluster_theory.comoving_volume(z, sky_area) * cluster_theory.mass_function(mass, z) * self.redshift_distribution.distribution() - * self.mass_distribution.distribution() + * self.mass_distribution.distribution(mass, z, mass_proxy) ) if average_on is None: @@ -60,7 +77,12 @@ def theory_prediction( def get_function_to_integrate( self, prediction: Callable[ - [npt.NDArray[np.float64], npt.NDArray[np.float64], float], + [ + npt.NDArray[np.float64], + npt.NDArray[np.float64], + npt.NDArray[np.float64], + float, + ], npt.NDArray[np.float64], ], ) -> Callable[ @@ -71,27 +93,28 @@ def numcosmo_wrapper( ) -> npt.NDArray[np.float64]: mass = int_args[:, 0] z = int_args[:, 1] + mass_proxy = int_args[:, 2] sky_area = extra_args[0] - return prediction(mass, z, sky_area) + return prediction(mass, z, mass_proxy, sky_area) return numcosmo_wrapper def evaluate_theory_prediction( self, cluster_theory: ClusterAbundance, - bin: NDimensionalBin, + this_bin: NDimensionalBin, sky_area: float, average_on: Optional[ClusterProperty] = None, ) -> float: - # Mock some fake mass to richness relation: - mass_low = bin.mass_proxy_edges[0] + 13 - mass_high = bin.mass_proxy_edges[1] + 13 + """Evaluate the theoretical prediction for the observable in the provided bin + using the Murata 2019 unbinned mass-richness relation and assuming perfectly + measured redshifts.""" self.integrator.integral_bounds = [ - (mass_low, mass_high), - bin.z_edges, + (cluster_theory.min_mass, cluster_theory.max_mass), + this_bin.z_edges, + this_bin.mass_proxy_edges, ] - print(mass_low, mass_high) - self.integrator.extra_args = np.array([sky_area], dtype=np.float64) + self.integrator.extra_args = np.array([*this_bin.mass_proxy_edges, sky_area]) theory_prediction = self.get_theory_prediction(cluster_theory, average_on) prediction_wrapper = self.get_function_to_integrate(theory_prediction) From 4fd99b4233cb4f189d5b92f8c565d27c3290e51b Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Tue, 9 Jan 2024 15:31:38 -0800 Subject: [PATCH 21/26] Adding more docstrings --- .../models/cluster/recipes/murata_binned_spec_z.py | 9 +++++++++ .../cluster/recipes/murata_unbinned_spec_z.py | 14 ++++++++++---- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/firecrown/models/cluster/recipes/murata_binned_spec_z.py b/firecrown/models/cluster/recipes/murata_binned_spec_z.py index 9573c867..e18c897c 100644 --- a/firecrown/models/cluster/recipes/murata_binned_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_binned_spec_z.py @@ -35,6 +35,10 @@ def get_theory_prediction( [npt.NDArray[np.float64], npt.NDArray[np.float64], Tuple[float, float], float], npt.NDArray[np.float64], ]: + """Returns a callable function that accepts mass, redshift, mass proxy limits, + and the sky area of your survey and returns the theoretical prediction for the + expected number of clusters.""" + def theory_prediction( mass: npt.NDArray[np.float64], z: npt.NDArray[np.float64], @@ -78,6 +82,11 @@ def get_function_to_integrate( npt.NDArray[np.float64], ], ) -> Callable[[npt.NDArray, npt.NDArray], npt.NDArray]: + """Returns a callable function that can be evaluated by an integrator. + + This function is responsible for mapping arguments from the numerical integrator + to the arguments of the theoretical prediction function.""" + def function_mapper( int_args: npt.NDArray, extra_args: npt.NDArray ) -> npt.NDArray[np.float64]: diff --git a/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py b/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py index fff0e3a4..30210b14 100644 --- a/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py @@ -40,7 +40,9 @@ def get_theory_prediction( ], npt.NDArray[np.float64], ]: - import pdb + """Returns a callable function that accepts mass, redshift, mass proxy, + and the sky area of your survey and returns the theoretical prediction for the + expected number of clusters.""" def theory_prediction( mass: npt.NDArray[np.float64], @@ -48,7 +50,6 @@ def theory_prediction( mass_proxy: npt.NDArray[np.float64], sky_area: float, ): - pdb.set_trace() prediction = ( cluster_theory.comoving_volume(z, sky_area) * cluster_theory.mass_function(mass, z) @@ -88,7 +89,12 @@ def get_function_to_integrate( ) -> Callable[ [npt.NDArray[np.float64], npt.NDArray[np.float64]], npt.NDArray[np.float64] ]: - def numcosmo_wrapper( + """Returns a callable function that can be evaluated by an integrator. + + This function is responsible for mapping arguments from the numerical integrator + to the arguments of the theoretical prediction function.""" + + def function_mapper( int_args: npt.NDArray[np.float64], extra_args: npt.NDArray[np.float64] ) -> npt.NDArray[np.float64]: mass = int_args[:, 0] @@ -97,7 +103,7 @@ def numcosmo_wrapper( sky_area = extra_args[0] return prediction(mass, z, mass_proxy, sky_area) - return numcosmo_wrapper + return function_mapper def evaluate_theory_prediction( self, From 38cb0a251265adf65edfa068bf46d3080dfec41a Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Wed, 10 Jan 2024 12:11:16 -0800 Subject: [PATCH 22/26] * Discovered bug in unit tests! (unbinned recipe) * Adding unbinned unit tests. --- .../cluster/recipes/murata_unbinned_spec_z.py | 5 +- .../test_murata_unbinned_spec_z.py | 184 ++++++++++++++++++ 2 files changed, 187 insertions(+), 2 deletions(-) create mode 100644 tests/cluster_recipes/test_murata_unbinned_spec_z.py diff --git a/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py b/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py index 30210b14..16da5fec 100644 --- a/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py @@ -14,7 +14,7 @@ # pylint: disable=R0801 -class MurataUnbinnedSpecZ(ClusterRecipe): +class MurataUnbinnedSpecZRecipe(ClusterRecipe): """Cluster recipe using the Murata 2019 unbinned mass-richness relation and assuming perfectly measured spec-zs.""" @@ -101,6 +101,7 @@ def function_mapper( z = int_args[:, 1] mass_proxy = int_args[:, 2] sky_area = extra_args[0] + return prediction(mass, z, mass_proxy, sky_area) return function_mapper @@ -120,7 +121,7 @@ def evaluate_theory_prediction( this_bin.z_edges, this_bin.mass_proxy_edges, ] - self.integrator.extra_args = np.array([*this_bin.mass_proxy_edges, sky_area]) + self.integrator.extra_args = np.array([sky_area]) theory_prediction = self.get_theory_prediction(cluster_theory, average_on) prediction_wrapper = self.get_function_to_integrate(theory_prediction) diff --git a/tests/cluster_recipes/test_murata_unbinned_spec_z.py b/tests/cluster_recipes/test_murata_unbinned_spec_z.py new file mode 100644 index 00000000..76207c14 --- /dev/null +++ b/tests/cluster_recipes/test_murata_unbinned_spec_z.py @@ -0,0 +1,184 @@ +"""Tests for the cluster abundance module.""" +from unittest.mock import Mock + +import numpy as np +import pyccl +import pytest + +from firecrown.models.cluster.abundance import ClusterAbundance +from firecrown.models.cluster.binning import NDimensionalBin +from firecrown.models.cluster.integrator.numcosmo_integrator import NumCosmoIntegrator +from firecrown.models.cluster.kernel import SpectroscopicRedshift +from firecrown.models.cluster.mass_proxy import MurataUnbinned +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe +from firecrown.models.cluster.recipes.murata_unbinned_spec_z import ( + MurataUnbinnedSpecZRecipe, +) + + +@pytest.fixture(name="cluster_abundance") +def fixture_cluster_abundance() -> ClusterAbundance: + hmf = pyccl.halos.MassFuncBocquet16() + cl_abundance = ClusterAbundance( + min_z=0, + max_z=2, + min_mass=13, + max_mass=17, + halo_mass_function=hmf, + ) + cl_abundance.update_ingredients(pyccl.CosmologyVanillaLCDM()) + return cl_abundance + + +@pytest.fixture(name="murata_unbinned_spec_z") +def fixture_murata_unbinned_spec_z() -> MurataUnbinnedSpecZRecipe: + cluster_recipe = MurataUnbinnedSpecZRecipe() + cluster_recipe.mass_distribution.mu_p0 = 3.0 + cluster_recipe.mass_distribution.mu_p1 = 0.86 + cluster_recipe.mass_distribution.mu_p2 = 0.0 + cluster_recipe.mass_distribution.sigma_p0 = 3.0 + cluster_recipe.mass_distribution.sigma_p1 = 0.7 + cluster_recipe.mass_distribution.sigma_p2 = 0.0 + return cluster_recipe + + +def test_murata_unbinned_spec_z_init(): + recipe = MurataUnbinnedSpecZRecipe() + + assert recipe is not None + assert isinstance(recipe, ClusterRecipe) + assert recipe.integrator is not None + assert isinstance(recipe.integrator, NumCosmoIntegrator) + assert recipe.redshift_distribution is not None + assert isinstance(recipe.redshift_distribution, SpectroscopicRedshift) + assert recipe.mass_distribution is not None + assert isinstance(recipe.mass_distribution, MurataUnbinned) + assert recipe.my_updatables is not None + assert len(recipe.my_updatables) == 1 + assert recipe.my_updatables[0] is recipe.mass_distribution + + +def test_get_theory_prediction_returns_value( + cluster_abundance: ClusterAbundance, + murata_unbinned_spec_z: MurataUnbinnedSpecZRecipe, +): + prediction = murata_unbinned_spec_z.get_theory_prediction(cluster_abundance) + + assert prediction is not None + assert callable(prediction) + + mass = np.linspace(13, 17, 2) + z = np.linspace(0.1, 1, 2) + mass_proxy = np.linspace(0, 5, 2) + sky_area = 360**2 + + result = prediction(mass, z, mass_proxy, sky_area) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + +def test_get_theory_prediction_with_average_returns_value( + cluster_abundance: ClusterAbundance, + murata_unbinned_spec_z: MurataUnbinnedSpecZRecipe, +): + mass = np.linspace(13, 17, 2) + z = np.linspace(0.1, 1, 2) + mass_proxy = np.linspace(0, 5, 2) + sky_area = 360**2 + + prediction = murata_unbinned_spec_z.get_theory_prediction( + cluster_abundance, average_on=ClusterProperty.MASS + ) + + assert prediction is not None + assert callable(prediction) + + result = prediction(mass, z, mass_proxy, sky_area) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + prediction = murata_unbinned_spec_z.get_theory_prediction( + cluster_abundance, average_on=ClusterProperty.REDSHIFT + ) + + assert prediction is not None + assert callable(prediction) + + result = prediction(mass, z, mass_proxy, sky_area) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + prediction = murata_unbinned_spec_z.get_theory_prediction( + cluster_abundance, average_on=(ClusterProperty.REDSHIFT | ClusterProperty.MASS) + ) + + assert prediction is not None + assert callable(prediction) + + result = prediction(mass, z, mass_proxy, sky_area) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + +def test_get_theory_prediction_throws_with_nonimpl_average( + cluster_abundance: ClusterAbundance, + murata_unbinned_spec_z: MurataUnbinnedSpecZRecipe, +): + prediction = murata_unbinned_spec_z.get_theory_prediction( + cluster_abundance, average_on=ClusterProperty.SHEAR + ) + + assert prediction is not None + assert callable(prediction) + + mass = np.linspace(13, 17, 2) + z = np.linspace(0.1, 1, 2) + mass_proxy = np.linspace(0, 5, 2) + sky_area = 360**2 + + with pytest.raises(NotImplementedError): + _ = prediction(mass, z, mass_proxy, sky_area) + + +def test_get_function_to_integrate_returns_value( + cluster_abundance: ClusterAbundance, + murata_unbinned_spec_z: MurataUnbinnedSpecZRecipe, +): + prediction = murata_unbinned_spec_z.get_theory_prediction(cluster_abundance) + function_to_integrate = murata_unbinned_spec_z.get_function_to_integrate(prediction) + + assert function_to_integrate is not None + assert callable(function_to_integrate) + + int_args = np.array([[13.0, 0.1, 0.1], [17.0, 1.0, 4.9]]) + extra_args = np.array([360**2]) + + result = function_to_integrate(int_args, extra_args) + assert isinstance(result, np.ndarray) + assert np.issubdtype(result.dtype, np.float64) + assert len(result) == 2 + assert np.all(result > 0) + + +def test_evaluates_theory_prediction_returns_value( + cluster_abundance: ClusterAbundance, + murata_unbinned_spec_z: MurataUnbinnedSpecZRecipe, +): + mock_bin = Mock(spec=NDimensionalBin) + mock_bin.mass_proxy_edges = (0, 5) + mock_bin.z_edges = (0, 1) + + prediction = murata_unbinned_spec_z.evaluate_theory_prediction( + cluster_abundance, mock_bin, 360**2 + ) + + assert prediction > 0 From c71f2ccef63605d241f7e55aea448402be8165d2 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Thu, 18 Jan 2024 13:45:23 -0800 Subject: [PATCH 23/26] Disabling duplicate-code in tests pylint rc file. --- tests/pylintrc | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/pylintrc b/tests/pylintrc index 5c53d424..977e8043 100644 --- a/tests/pylintrc +++ b/tests/pylintrc @@ -19,6 +19,7 @@ recursive=yes # multiple time (only on the command line, not in the configuration file where # it should appear only once). See also the "--disable" option for examples. enable=c-extension-no-member +disable=duplicate-code [MISCELLANEOUS] From 68bd259e7702d66bfedd0426f160a58be7414c4a Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Thu, 18 Jan 2024 14:11:39 -0800 Subject: [PATCH 24/26] Deleting custom implementation of __repr__ since it wasn't implemented correctly --- firecrown/models/cluster/binning.py | 3 --- tests/test_cluster_binning.py | 8 -------- 2 files changed, 11 deletions(-) diff --git a/firecrown/models/cluster/binning.py b/firecrown/models/cluster/binning.py index 272a060f..7220fa1c 100644 --- a/firecrown/models/cluster/binning.py +++ b/firecrown/models/cluster/binning.py @@ -35,9 +35,6 @@ def mass_proxy_edges(self) -> Tuple[float, float]: def __str__(self) -> str: return f"[{self.z_edges}, {self.mass_proxy_edges}]\n" - def __repr__(self) -> str: - return f"[{self.z_edges}, {self.mass_proxy_edges}]\n" - class SaccBin(NDimensionalBin[sacc.BaseTracer]): """An implementation of the N dimensional bin using sacc tracers.""" diff --git a/tests/test_cluster_binning.py b/tests/test_cluster_binning.py index 4684ec26..ad665603 100644 --- a/tests/test_cluster_binning.py +++ b/tests/test_cluster_binning.py @@ -7,14 +7,6 @@ from firecrown.models.cluster.binning import NDimensionalBin, SaccBin, TupleBin -def test_bin_repr(): - tracer_z = sacc.tracers.BinZTracer("", 0, 1) - tracer_lambda = sacc.tracers.BinRichnessTracer("", 4, 5) - - sb = SaccBin([tracer_z, tracer_lambda]) - assert repr(sb) == "[(0, 1), (4, 5)]\n" - - def test_bin_str(): tracer_z = sacc.tracers.BinZTracer("", 0, 1) tracer_lambda = sacc.tracers.BinRichnessTracer("", 4, 5) From aa580b3bbd212aefabb7c56734cb8ae1546b9c8f Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Thu, 18 Jan 2024 15:09:14 -0800 Subject: [PATCH 25/26] * Added custom pylint plugin for duplicate-code * Added updated pylint to use the custom plugins through the rc file * Removing custom disable duplicate code warning --- .../cluster/recipes/murata_binned_spec_z.py | 1 - .../cluster/recipes/murata_unbinned_spec_z.py | 1 - firecrown/models/pylintrc | 3 +++ pylint_plugins/README.md | 10 +++++++++ pylint_plugins/duplicate_code.py | 21 +++++++++++++++++++ pylintrc | 3 +++ 6 files changed, 37 insertions(+), 2 deletions(-) create mode 100644 pylint_plugins/README.md create mode 100644 pylint_plugins/duplicate_code.py diff --git a/firecrown/models/cluster/recipes/murata_binned_spec_z.py b/firecrown/models/cluster/recipes/murata_binned_spec_z.py index e18c897c..bdfbf9ff 100644 --- a/firecrown/models/cluster/recipes/murata_binned_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_binned_spec_z.py @@ -13,7 +13,6 @@ from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe -# pylint: disable=R0801 class MurataBinnedSpecZRecipe(ClusterRecipe): """Cluster recipe using the Murata 2019 binned mass-richness relation and assuming perfectly measured spec-zs.""" diff --git a/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py b/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py index 16da5fec..472a9b34 100644 --- a/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py +++ b/firecrown/models/cluster/recipes/murata_unbinned_spec_z.py @@ -13,7 +13,6 @@ from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe -# pylint: disable=R0801 class MurataUnbinnedSpecZRecipe(ClusterRecipe): """Cluster recipe using the Murata 2019 unbinned mass-richness relation and assuming perfectly measured spec-zs.""" diff --git a/firecrown/models/pylintrc b/firecrown/models/pylintrc index 502ff948..abd8b0a6 100644 --- a/firecrown/models/pylintrc +++ b/firecrown/models/pylintrc @@ -12,6 +12,9 @@ py-version=3.9 # Discover python modules and packages in the file system subtree. recursive=yes +# Add custom pylint plugins +load-plugins=pylint_plugins.duplicate_code + [MESSAGES CONTROL] # Enable the message, report, category or checker with the given id(s). You can diff --git a/pylint_plugins/README.md b/pylint_plugins/README.md new file mode 100644 index 00000000..fe9bb619 --- /dev/null +++ b/pylint_plugins/README.md @@ -0,0 +1,10 @@ +### Custom Pylint Plugins +Pylint does not have any built in functionality to selectively apply rules to subdirectories. + +As an example, if we want to ignore a pylint warning in only one class/directory/module, we would need to create a custom `pylintrc` file and manually exclude this file from the parent pylint runs. + +Custom plugins allow us to selectively ignore warnings. This directory holds the custom plugins we load into pylint through the `pylintrc` files. + +### Table of Contents + +* `duplicate_code`: A custom plugin that will ignore the `duplicate-code` warning. Currently the only classes we allow to ignore this message is `firecrown.models.cluster.recipes`. \ No newline at end of file diff --git a/pylint_plugins/duplicate_code.py b/pylint_plugins/duplicate_code.py new file mode 100644 index 00000000..7191a921 --- /dev/null +++ b/pylint_plugins/duplicate_code.py @@ -0,0 +1,21 @@ +from astroid import MANAGER +from astroid import scoped_nodes +from pylint.lint import PyLinter + + +def register(linter: PyLinter): + pass + + +def transform(mod): + if "firecrown.models.cluster.recipes." not in mod.name: + return + print(mod.name) + c = mod.stream().read() + c = b"# pylint: disable=duplicate-code\n" + c + + # pylint will read from `.file_bytes` attribute later when tokenization + mod.file_bytes = c + + +MANAGER.register_transform(scoped_nodes.Module, transform) diff --git a/pylintrc b/pylintrc index fdc3102e..59b3b0b3 100644 --- a/pylintrc +++ b/pylintrc @@ -12,6 +12,9 @@ py-version=3.9 # Discover python modules and packages in the file system subtree. recursive=yes +# Add custom pylint plugins +load-plugins=pylint_plugins.duplicate_code + [MESSAGES CONTROL] # Enable the message, report, category or checker with the given id(s). You can From 516d37a1edf2fc1df0bdce997d325ce533f9ad87 Mon Sep 17 00:00:00 2001 From: Matt Kwiecien Date: Thu, 18 Jan 2024 15:13:31 -0800 Subject: [PATCH 26/26] Removing debug statement. --- pylint_plugins/duplicate_code.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylint_plugins/duplicate_code.py b/pylint_plugins/duplicate_code.py index 7191a921..3e95d734 100644 --- a/pylint_plugins/duplicate_code.py +++ b/pylint_plugins/duplicate_code.py @@ -10,7 +10,7 @@ def register(linter: PyLinter): def transform(mod): if "firecrown.models.cluster.recipes." not in mod.name: return - print(mod.name) + c = mod.stream().read() c = b"# pylint: disable=duplicate-code\n" + c