From 0e7f39a518d1d0ca9d1c93cc338bd044829acc03 Mon Sep 17 00:00:00 2001
From: eduardojsbarroso <eduardojsbarroso@gmail.com>
Date: Fri, 31 Jan 2025 17:04:11 +0100
Subject: [PATCH 1/7] Modified selection function to be updatable
---
firecrown/models/cluster/kernel.py | 102 +++++++--
.../test_murata_binned_spec_z_selection.py | 199 ++++++++++++++++++
tests/test_cluster_kernels.py | 104 +++++----
3 files changed, 346 insertions(+), 59 deletions(-)
create mode 100644 tests/cluster_recipes/test_murata_binned_spec_z_selection.py
diff --git a/firecrown/models/cluster/kernel.py b/firecrown/models/cluster/kernel.py
index 3b7b6d78..cf10453f 100644
--- a/firecrown/models/cluster/kernel.py
+++ b/firecrown/models/cluster/kernel.py
@@ -8,6 +8,15 @@
import numpy.typing as npt
import numpy as np
+from firecrown import parameters
+from firecrown.updatable import Updatable
+
+
+REDMAPPER_DEFAULT_AC_NC = 0.38
+REDMAPPER_DEFAULT_BC_NC = 1.2634
+REDMAPPER_DEFAULT_AC_MC = 13.31
+REDMAPPER_DEFAULT_BC_MC = 0.2025
+
class KernelType(Enum):
"""The kernels that can be included in the cluster abundance integrand."""
@@ -20,47 +29,98 @@ class KernelType(Enum):
PURITY = 6
-class Completeness:
+class Completeness(Updatable):
"""The completeness kernel for the numcosmo simulated survey.
This kernel will affect the integrand by accounting for the incompleteness
of a cluster selection.
"""
+ def __init__(
+ self,
+ ):
+ super().__init__()
+ # Updatable parameters
+ self.ac_nc = parameters.register_new_updatable_parameter(
+ default_value=REDMAPPER_DEFAULT_AC_NC
+ )
+ self.bc_nc = parameters.register_new_updatable_parameter(
+ default_value=REDMAPPER_DEFAULT_BC_NC
+ )
+ self.ac_mc = parameters.register_new_updatable_parameter(
+ default_value=REDMAPPER_DEFAULT_AC_MC
+ )
+ self.bc_mc = parameters.register_new_updatable_parameter(
+ default_value=REDMAPPER_DEFAULT_BC_MC
+ )
+
+ def _mc(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
+ ac_mc = self.ac_mc
+ bc_mc = self.bc_mc
+ log_mc = ac_mc + bc_mc * (1.0 + z)
+ mc = 10.0**log_mc
+ return mc.astype(np.float64)
+
+ def _nc(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
+ ac_nc = self.ac_nc
+ bc_nc = self.bc_nc
+ nc = ac_nc + bc_nc * (1.0 + z)
+ assert isinstance(nc, np.ndarray)
+ return nc
+
def distribution(
self,
- mass: npt.NDArray[np.float64],
+ log_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
- b_mc = 0.2025
- log_mc = a_mc + b_mc * (1.0 + z)
- nc = a_nc + b_nc * (1.0 + z)
- completeness = (mass / log_mc) ** nc / ((mass / log_mc) ** nc + 1.0)
+ mc = self._mc(z)
+ mass = 10.0**log_mass
+ nc = self._nc(z)
+ completeness = (mass / mc) ** nc / ((mass / mc) ** nc + 1.0)
assert isinstance(completeness, np.ndarray)
return completeness
-class Purity:
+REDMAPPER_DEFAULT_AP_NC = 3.9193
+REDMAPPER_DEFAULT_BP_NC = -0.3323
+REDMAPPER_DEFAULT_AP_RC = 1.1839
+REDMAPPER_DEFAULT_BP_RC = -0.4077
+
+
+class Purity(Updatable):
"""The purity kernel for the numcosmo simulated survey.
This kernel will affect the integrand by accounting for the purity
of a cluster selection.
"""
- def _ln_rc(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
- a_rc = 2.2183
- b_rc = -0.6592
- ln_rc = a_rc + b_rc * (1.0 + z)
- return ln_rc.astype(np.float64)
+ def __init__(self):
+ super().__init__()
+ self.ap_nc = parameters.register_new_updatable_parameter(
+ default_value=REDMAPPER_DEFAULT_AP_NC
+ )
+ self.bp_nc = parameters.register_new_updatable_parameter(
+ default_value=REDMAPPER_DEFAULT_BP_NC
+ )
+ self.ap_rc = parameters.register_new_updatable_parameter(
+ default_value=REDMAPPER_DEFAULT_AP_RC
+ )
+ self.bp_rc = parameters.register_new_updatable_parameter(
+ default_value=REDMAPPER_DEFAULT_BP_RC
+ )
+
+ def _rc(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
+ ap_rc = self.ap_rc
+ bp_rc = self.bp_rc
+ log_rc = ap_rc + bp_rc * (1.0 + z)
+ rc = 10**log_rc
+ return rc.astype(np.float64)
def _nc(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
- b_nc = np.log(10) * 0.3527
- a_nc = np.log(10) * 0.8612
- nc = a_nc + b_nc * (1.0 + z)
+ bp_nc = self.bp_nc
+ ap_nc = self.ap_nc
+ nc = ap_nc + bp_nc * (1.0 + z)
assert isinstance(nc, np.ndarray)
return nc
@@ -73,11 +133,11 @@ def distribution(
"""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)
+ r = 10**mean_mass
else:
- ln_r = np.log(10**mass_proxy)
+ r = 10**mass_proxy
- r_over_rc = ln_r / self._ln_rc(z)
+ r_over_rc = r / self._rc(z)
purity = (r_over_rc) ** self._nc(z) / (r_over_rc ** self._nc(z) + 1.0)
assert isinstance(purity, np.ndarray)
diff --git a/tests/cluster_recipes/test_murata_binned_spec_z_selection.py b/tests/cluster_recipes/test_murata_binned_spec_z_selection.py
new file mode 100644
index 00000000..e92ebc41
--- /dev/null
+++ b/tests/cluster_recipes/test_murata_binned_spec_z_selection.py
@@ -0,0 +1,199 @@
+"""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_selection import (
+ MurataBinnedSpecZSelectionRecipe,
+)
+
+
+@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_selection")
+def fixture_murata_binned_spec_z() -> MurataBinnedSpecZSelectionRecipe:
+ cluster_recipe = MurataBinnedSpecZSelectionRecipe()
+ 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
+ cluster_recipe.purity_distribution.ap_rc = 1.1839
+ cluster_recipe.purity_distribution.bp_rc = -0.4077
+ cluster_recipe.purity_distribution.ap_nc = 3.9193
+ cluster_recipe.purity_distribution.bp_nc = -0.3323
+ cluster_recipe.completeness_distribution.ac_mc = 13.31
+ cluster_recipe.completeness_distribution.bc_mc = 0.2025
+ cluster_recipe.completeness_distribution.ac_nc = 0.38
+ cluster_recipe.completeness_distribution.bc_nc = 1.2634
+ return cluster_recipe
+
+
+def test_murata_binned_spec_z_init():
+ recipe = MurataBinnedSpecZSelectionRecipe()
+
+ 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.completeness_distribution is not None
+ assert recipe.purity_distribution is not None
+ assert recipe.my_updatables is not None
+ assert len(recipe.my_updatables) == 3
+ assert recipe.my_updatables[0] is recipe.mass_distribution
+ assert recipe.my_updatables[1] is recipe.completeness_distribution
+ assert recipe.my_updatables[2] is recipe.purity_distribution
+
+
+def test_get_theory_prediction_returns_value(
+ cluster_abundance: ClusterAbundance,
+ murata_binned_spec_z_selection: MurataBinnedSpecZSelectionRecipe,
+):
+ prediction = murata_binned_spec_z_selection.get_theory_prediction(cluster_abundance)
+
+ assert prediction is not None
+ assert callable(prediction)
+
+ mass = np.linspace(13, 17, 2, dtype=np.float64)
+ z = np.linspace(0.1, 1, 2, dtype=np.float64)
+ 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_selection: MurataBinnedSpecZSelectionRecipe,
+):
+ mass = np.linspace(13, 17, 2, dtype=np.float64)
+ z = np.linspace(0.1, 1, 2, dtype=np.float64)
+ mass_proxy_limits = (0, 5)
+ sky_area = 360**2
+
+ prediction = murata_binned_spec_z_selection.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_selection.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_selection.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_selection: MurataBinnedSpecZSelectionRecipe,
+):
+ prediction = murata_binned_spec_z_selection.get_theory_prediction(
+ cluster_abundance, average_on=ClusterProperty.SHEAR
+ )
+
+ assert prediction is not None
+ assert callable(prediction)
+
+ mass = np.linspace(13, 17, 2, dtype=np.float64)
+ z = np.linspace(0.1, 1, 2, dtype=np.float64)
+ 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_selection: MurataBinnedSpecZSelectionRecipe,
+):
+ prediction = murata_binned_spec_z_selection.get_theory_prediction(cluster_abundance)
+ function_to_integrate = murata_binned_spec_z_selection.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_selection: MurataBinnedSpecZSelectionRecipe,
+):
+ mock_bin = Mock(spec=NDimensionalBin)
+ mock_bin.mass_proxy_edges = (0, 5)
+ mock_bin.z_edges = (0, 1)
+
+ prediction = murata_binned_spec_z_selection.evaluate_theory_prediction(
+ cluster_abundance, mock_bin, 360**2
+ )
+
+ assert prediction > 0
diff --git a/tests/test_cluster_kernels.py b/tests/test_cluster_kernels.py
index 36e83d4e..5740d88c 100644
--- a/tests/test_cluster_kernels.py
+++ b/tests/test_cluster_kernels.py
@@ -23,12 +23,28 @@ def test_create_mass_kernel():
def test_create_completeness_kernel():
ck = Completeness()
+ ck.ac_mc = 13.31
+ ck.bc_mc = 0.2025
+ ck.ac_nc = 0.38
+ ck.bc_nc = 1.2634
assert ck is not None
+ assert ck.ac_mc == 13.31
+ assert ck.bc_mc == 0.2025
+ assert ck.ac_nc == 0.38
+ assert ck.bc_nc == 1.2634
def test_create_purity_kernel():
pk = Purity()
+ pk.ap_nc = 3.9193
+ pk.bp_nc = -0.3323
+ pk.ap_rc = 1.1839
+ pk.bp_rc = -0.4077
assert pk is not None
+ assert pk.ap_nc == 3.9193
+ assert pk.bp_nc == -0.3323
+ assert pk.ap_rc == 1.1839
+ assert pk.bp_rc == -0.4077
def test_spec_z_distribution():
@@ -44,56 +60,64 @@ def test_true_mass_distribution():
@pytest.mark.precision_sensitive
def test_purity_distribution():
pk = Purity()
- mass_proxy = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
+ pk.ap_nc = 3.9193
+ pk.bp_nc = -0.3323
+ pk.ap_rc = 1.1839
+ pk.bp_rc = -0.4077
+ mass_proxy = np.linspace(0.0, 2.5, 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_limits = (1.0, 10.0)
truth = np.array(
[
- 0.77657274,
- 0.96966127,
- 0.99286409,
- 0.99780586,
- 0.999224,
- 0.99970302,
- 0.99988111,
- 0.99995125,
- 0.99997982,
- 0.99999166,
- ]
+ 0.00242882,
+ 0.03294582,
+ 0.3122527,
+ 0.85213252,
+ 0.98584893,
+ 0.99875485,
+ 0.99988632,
+ 0.99998911,
+ 0.99999891,
+ 0.99999988,
+ ],
+ dtype=np.float64,
)
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)
+ assert ref == pytest.approx(true, rel=1e-5, abs=0.0)
@pytest.mark.precision_sensitive
def test_purity_distribution_uses_mean():
pk = Purity()
+ pk.ap_nc = 3.9193
+ pk.bp_nc = -0.3323
+ pk.ap_rc = 1.1839
+ pk.bp_rc = -0.4077
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)
+ mass_proxy = np.array([-1.0])
+ mass_proxy_limits = (0.0, 2.0)
truth = np.array(
[
- 0.9978693724040568,
- 0.9984319673134954,
- 0.9988620014089232,
- 0.9991864843696077,
- 0.9994279315032029,
- 0.999604893383804,
- 0.9997324678841709,
- 0.9998227843987537,
- 0.9998854531462606,
- 0.9999279749997235,
+ 0.89705651,
+ 0.92238419,
+ 0.94154163,
+ 0.95593305,
+ 0.96670586,
+ 0.97476117,
+ 0.98078884,
+ 0.98530847,
+ 0.98870753,
+ 0.99127329,
],
dtype=np.float64,
)
-
- purity = pk.distribution(z, mass_proxy.astype(np.float64), mass_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)
@@ -102,21 +126,25 @@ def test_purity_distribution_uses_mean():
@pytest.mark.precision_sensitive
def test_completeness_distribution():
ck = Completeness()
- mass = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
+ ck.ac_mc = 13.31
+ ck.bc_mc = 0.2025
+ ck.ac_nc = 0.38
+ ck.bc_nc = 1.2634
+ mass = np.linspace(13.0, 15.0, 10)
z = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
truth = np.array(
[
- 0.0056502277493542,
- 0.01896566878380423,
- 0.03805597500308377,
- 0.06224888967250564,
- 0.09124569979282898,
- 0.12486247682690908,
- 0.16290218589569144,
- 0.20507815091349266,
- 0.2509673905442634,
- 0.2999886170051561,
+ 0.10239024,
+ 0.19090539,
+ 0.35438466,
+ 0.58952617,
+ 0.80866296,
+ 0.93327968,
+ 0.98115635,
+ 0.99543348,
+ 0.99902667,
+ 0.99981606,
]
)
From 189885e958083b19bff24d6e8a801203249592ab Mon Sep 17 00:00:00 2001
From: eduardojsbarroso <eduardojsbarroso@gmail.com>
Date: Fri, 31 Jan 2025 17:10:03 +0100
Subject: [PATCH 2/7] Fixed mypy error
---
firecrown/models/cluster/kernel.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/firecrown/models/cluster/kernel.py b/firecrown/models/cluster/kernel.py
index cf10453f..40802bd2 100644
--- a/firecrown/models/cluster/kernel.py
+++ b/firecrown/models/cluster/kernel.py
@@ -135,7 +135,7 @@ def distribution(
mean_mass = (mass_proxy_limits[0] + mass_proxy_limits[1]) / 2
r = 10**mean_mass
else:
- r = 10**mass_proxy
+ r = np.power(10, mass_proxy)
r_over_rc = r / self._rc(z)
From 381ed01183b858e2774a15ae1b383772f2449601 Mon Sep 17 00:00:00 2001
From: eduardojsbarroso <eduardojsbarroso@gmail.com>
Date: Fri, 31 Jan 2025 17:17:44 +0100
Subject: [PATCH 3/7] typo
---
firecrown/models/cluster/kernel.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/firecrown/models/cluster/kernel.py b/firecrown/models/cluster/kernel.py
index 40802bd2..c2f6ee7e 100644
--- a/firecrown/models/cluster/kernel.py
+++ b/firecrown/models/cluster/kernel.py
@@ -135,7 +135,7 @@ def distribution(
mean_mass = (mass_proxy_limits[0] + mass_proxy_limits[1]) / 2
r = 10**mean_mass
else:
- r = np.power(10, mass_proxy)
+ r = np.power(10., mass_proxy).item()
r_over_rc = r / self._rc(z)
From 3e31a31b9423e4de5bafdee20c356563d9ede7fa Mon Sep 17 00:00:00 2001
From: eduardojsbarroso <eduardojsbarroso@gmail.com>
Date: Fri, 31 Jan 2025 17:20:09 +0100
Subject: [PATCH 4/7] Black
---
firecrown/models/cluster/kernel.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/firecrown/models/cluster/kernel.py b/firecrown/models/cluster/kernel.py
index c2f6ee7e..e614de2f 100644
--- a/firecrown/models/cluster/kernel.py
+++ b/firecrown/models/cluster/kernel.py
@@ -135,7 +135,7 @@ def distribution(
mean_mass = (mass_proxy_limits[0] + mass_proxy_limits[1]) / 2
r = 10**mean_mass
else:
- r = np.power(10., mass_proxy).item()
+ r = np.power(10.0, mass_proxy).item()
r_over_rc = r / self._rc(z)
From aef072a5c94fd6359c3ec64310f29ac5ce9d213f Mon Sep 17 00:00:00 2001
From: eduardojsbarroso <eduardojsbarroso@gmail.com>
Date: Fri, 31 Jan 2025 17:31:12 +0100
Subject: [PATCH 5/7] Test fix
---
tests/cluster_recipes/test_murata_binned_spec_z_selection.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/tests/cluster_recipes/test_murata_binned_spec_z_selection.py b/tests/cluster_recipes/test_murata_binned_spec_z_selection.py
index e92ebc41..07a1d804 100644
--- a/tests/cluster_recipes/test_murata_binned_spec_z_selection.py
+++ b/tests/cluster_recipes/test_murata_binned_spec_z_selection.py
@@ -33,7 +33,7 @@ def fixture_cluster_abundance() -> ClusterAbundance:
@pytest.fixture(name="murata_binned_spec_z_selection")
-def fixture_murata_binned_spec_z() -> MurataBinnedSpecZSelectionRecipe:
+def fixture_murata_binned_spec_z_selection() -> MurataBinnedSpecZSelectionRecipe:
cluster_recipe = MurataBinnedSpecZSelectionRecipe()
cluster_recipe.mass_distribution.mu_p0 = 3.0
cluster_recipe.mass_distribution.mu_p1 = 0.86
@@ -52,7 +52,7 @@ def fixture_murata_binned_spec_z() -> MurataBinnedSpecZSelectionRecipe:
return cluster_recipe
-def test_murata_binned_spec_z_init():
+def test_murata_binned_spec_z_selection_init():
recipe = MurataBinnedSpecZSelectionRecipe()
assert recipe is not None
From 08a647680714d5e39ed9e835ba0e9113dcde4560 Mon Sep 17 00:00:00 2001
From: eduardojsbarroso <eduardojsbarroso@gmail.com>
Date: Tue, 11 Feb 2025 15:49:46 +0100
Subject: [PATCH 6/7] MInor fixes and recipe creation
---
firecrown/models/cluster/kernel.py | 4 +-
firecrown/models/cluster/mass_proxy.py | 2 +-
.../recipes/murata_binned_spec_z_selection.py | 143 ++++++++++++++++++
3 files changed, 146 insertions(+), 3 deletions(-)
create mode 100644 firecrown/models/cluster/recipes/murata_binned_spec_z_selection.py
diff --git a/firecrown/models/cluster/kernel.py b/firecrown/models/cluster/kernel.py
index e614de2f..3e0763b1 100644
--- a/firecrown/models/cluster/kernel.py
+++ b/firecrown/models/cluster/kernel.py
@@ -128,14 +128,14 @@ def distribution(
self,
z: npt.NDArray[np.float64],
mass_proxy: npt.NDArray[np.float64],
- mass_proxy_limits: tuple[float, float],
+ mass_proxy_limits: tuple[float, float] = None,
) -> 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
r = 10**mean_mass
else:
- r = np.power(10.0, mass_proxy).item()
+ r = np.power(10.0, mass_proxy)
r_over_rc = r / self._rc(z)
diff --git a/firecrown/models/cluster/mass_proxy.py b/firecrown/models/cluster/mass_proxy.py
index 3fc105cf..4a0a1e1b 100644
--- a/firecrown/models/cluster/mass_proxy.py
+++ b/firecrown/models/cluster/mass_proxy.py
@@ -93,7 +93,7 @@ def _distribution_unbinned(
normalization = 1 / np.sqrt(2 * np.pi * proxy_sigma**2)
result = normalization * np.exp(
- -0.5 * ((mass_proxy - proxy_mean) / proxy_sigma) ** 2
+ -0.5 * ((mass_proxy * np.log(10) - proxy_mean) / proxy_sigma) ** 2
)
assert isinstance(result, np.ndarray)
diff --git a/firecrown/models/cluster/recipes/murata_binned_spec_z_selection.py b/firecrown/models/cluster/recipes/murata_binned_spec_z_selection.py
new file mode 100644
index 00000000..3e5e34e4
--- /dev/null
+++ b/firecrown/models/cluster/recipes/murata_binned_spec_z_selection.py
@@ -0,0 +1,143 @@
+"""Module for defining the classes used in the MurataBinnedSpecZ cluster recipe."""
+
+from typing import Callable
+
+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
+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.kernel import Completeness, Purity
+
+
+class MurataBinnedSpecZSelectionRecipe(ClusterRecipe):
+ """Cluster recipe with Murata19 mass-richness and spec-zs.
+
+ This recipe uses the Murata 2019 binned mass-richness relation and assumes
+ perfectly measured spec-zs.
+ """
+
+ def __init__(self) -> None:
+ super().__init__()
+
+ self.integrator = NumCosmoIntegrator()
+ self.redshift_distribution = SpectroscopicRedshift()
+ pivot_mass, pivot_redshift = 14.3, 0.5
+ self.mass_distribution = MurataUnbinned(pivot_mass, pivot_redshift)
+ self.completeness_distribution = Completeness()
+ self.purity_distribution = Purity()
+ self.my_updatables.append(self.mass_distribution)
+ self.my_updatables.append(self.completeness_distribution)
+ self.my_updatables.append(self.purity_distribution)
+
+ def get_theory_prediction(
+ self,
+ cluster_theory: ClusterAbundance,
+ average_on: None | ClusterProperty = None,
+ ) -> Callable[
+ [npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64], float],
+ npt.NDArray[np.float64],
+ ]:
+ """Get a callable that evaluates a cluster theory prediction.
+
+ 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],
+ mass_proxy: npt.NDArray[np.float64],
+ sky_area: float,
+ ):
+ 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 / np.log(10))
+ * self.completeness_distribution.distribution(mass, z)
+ / self.purity_distribution.distribution(
+ z, mass_proxy=mass_proxy
+ )
+ )
+
+ 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[
+ [
+ npt.NDArray[np.float64],
+ npt.NDArray[np.float64],
+ npt.NDArray[np.float64],
+ float,
+ ],
+ 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]:
+ mass = int_args[:, 0]
+ 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
+
+ def evaluate_theory_prediction(
+ self,
+ cluster_theory: ClusterAbundance,
+ this_bin: NDimensionalBin,
+ sky_area: float,
+ average_on: None | ClusterProperty = None,
+ ) -> float:
+ """Evaluate the theory prediction for this cluster recipe.
+
+ 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),
+ this_bin.z_edges, np.log(10) * np.array(this_bin.mass_proxy_edges),
+ ]
+ 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)
+
+ counts = self.integrator.integrate(prediction_wrapper)
+
+ return counts
From 6ef1ca4a6bb1d3164b296d46d09070e01345eff8 Mon Sep 17 00:00:00 2001
From: Marc Paterno <paterno@fnal.gov>
Date: Fri, 14 Feb 2025 14:36:29 -0600
Subject: [PATCH 7/7] Apply black
---
.../recipes/murata_binned_spec_z_selection.py | 18 +++++++++++-------
setup.py | 3 +--
2 files changed, 12 insertions(+), 9 deletions(-)
diff --git a/firecrown/models/cluster/recipes/murata_binned_spec_z_selection.py b/firecrown/models/cluster/recipes/murata_binned_spec_z_selection.py
index 3e5e34e4..4fff475b 100644
--- a/firecrown/models/cluster/recipes/murata_binned_spec_z_selection.py
+++ b/firecrown/models/cluster/recipes/murata_binned_spec_z_selection.py
@@ -40,7 +40,12 @@ def get_theory_prediction(
cluster_theory: ClusterAbundance,
average_on: None | ClusterProperty = None,
) -> Callable[
- [npt.NDArray[np.float64], 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],
]:
"""Get a callable that evaluates a cluster theory prediction.
@@ -62,9 +67,7 @@ def theory_prediction(
* self.redshift_distribution.distribution()
* self.mass_distribution.distribution(mass, z, mass_proxy / np.log(10))
* self.completeness_distribution.distribution(mass, z)
- / self.purity_distribution.distribution(
- z, mass_proxy=mass_proxy
- )
+ / self.purity_distribution.distribution(z, mass_proxy=mass_proxy)
)
if average_on is None:
@@ -108,8 +111,8 @@ def function_mapper(
) -> npt.NDArray[np.float64]:
mass = int_args[:, 0]
z = int_args[:, 1]
- mass_proxy = int_args[:,2]
-
+ mass_proxy = int_args[:, 2]
+
sky_area = extra_args[0]
return prediction(mass, z, mass_proxy, sky_area)
@@ -131,7 +134,8 @@ def evaluate_theory_prediction(
"""
self.integrator.integral_bounds = [
(cluster_theory.min_mass, cluster_theory.max_mass),
- this_bin.z_edges, np.log(10) * np.array(this_bin.mass_proxy_edges),
+ this_bin.z_edges,
+ np.log(10) * np.array(this_bin.mass_proxy_edges),
]
self.integrator.extra_args = np.array([sky_area])
diff --git a/setup.py b/setup.py
index 177abd0a..970808c8 100644
--- a/setup.py
+++ b/setup.py
@@ -1,5 +1,4 @@
-"""Setup script for firecrown.
-"""
+"""Setup script for firecrown."""
from setuptools import setup