fixes to firecrown API and tests to align with final CCL v2

docs/conf.py

@@ -133,8 +133,8 @@
 ]
 
 # Copied from github.com/sanderslab/magellanmapper:
-## automate building API .rst files, necessary for ReadTheDocs, as inspired by:
-## https://github.com/readthedocs/readthedocs.org/issues/1139#issuecomment-398083449
+# automate building API .rst files, necessary for ReadTheDocs, as inspired by:
+# https://github.com/readthedocs/readthedocs.org/issues/1139#issuecomment-398083449
 
 
 def run_apidoc(_):
@@ -147,7 +147,6 @@ def run_apidoc(_):
         "-e",
         "-E",
         "-T",
-        #        "--implicit-namespaces",
         "-d",
         "1",
         "-o",
@@ -166,7 +165,3 @@ def run_apidoc(_):
 
         argv.insert(0, apidoc.__file__)
         apidoc.main(argv)
-
-
-# def setup(app):
-#    app.connect('builder-inited', run_apidoc)

environment.yml

@@ -8,6 +8,7 @@ dependencies:
   - cosmosis-build-standard-library
   - numcosmo >=0.16.0
   - coverage
+  - cython < 3.0.0
   - dill
   - fitsio
   - flake8
@@ -24,7 +25,7 @@ dependencies:
     - cobaya
   - portalocker
   - pybobyqa
-  - pyccl
+  - pyccl >= 2.8.0
   - pylint
   - pytest
   - pyyaml

examples/des_y1_3x2pt/des_y1_3x2pt_PT.py

@@ -86,10 +86,9 @@ def build_likelihood(_) -> Tuple[Likelihood, ModelingTools]:
     )
 
     # Create the likelihood from the statistics
-    pt_calculator = pyccl.nl_pt.PTCalculator(
+    pt_calculator = pyccl.nl_pt.EulerianPTCalculator(
         with_NC=True,
         with_IA=True,
-        with_dd=False,
         log10k_min=-4,
         log10k_max=2,
         nk_per_decade=20,
@@ -154,23 +153,28 @@ def run_likelihood() -> None:
     )
 
     # Code that creates a Pk2D object:
-    ptc = pyccl.nl_pt.PTCalculator(
-        with_NC=True, with_IA=True, log10k_min=-4, log10k_max=2, nk_per_decade=20
+    ptc = pyccl.nl_pt.EulerianPTCalculator(
+        with_NC=True,
+        with_IA=True,
+        log10k_min=-4,
+        log10k_max=2,
+        nk_per_decade=20,
+        cosmo=ccl_cosmo,
     )
     ptt_i = pyccl.nl_pt.PTIntrinsicAlignmentTracer(
         c1=(z, c_1), c2=(z, c_2), cdelta=(z, c_d)
     )
     ptt_m = pyccl.nl_pt.PTMatterTracer()
     ptt_g = pyccl.nl_pt.PTNumberCountsTracer(b1=b_1, b2=b_2, bs=b_s)
     # IA
-    pk_im = pyccl.nl_pt.get_pt_pk2d(ccl_cosmo, ptt_i, tracer2=ptt_m, ptc=ptc)
-    pk_ii = pyccl.nl_pt.get_pt_pk2d(ccl_cosmo, ptt_i, ptc=ptc)
-    pk_gi = pyccl.nl_pt.get_pt_pk2d(ccl_cosmo, ptt_g, tracer2=ptt_i, ptc=ptc)
+    pk_im = ptc.get_biased_pk2d(ccl_cosmo, ptt_i, tracer2=ptt_m)
+    pk_ii = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_i, tracer2=ptt_i)
+    pk_gi = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_i)
     # Galaxies
-    pk_gm = pyccl.nl_pt.get_pt_pk2d(ccl_cosmo, ptt_g, tracer2=ptt_m, ptc=ptc)
-    pk_gg = pyccl.nl_pt.get_pt_pk2d(ccl_cosmo, ptt_g, ptc=ptc)
+    pk_gm = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_m)
+    pk_gg = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_g)
     # Magnification: just a matter-matter P(k)
-    pk_mm = pyccl.nl_pt.get_pt_pk2d(ccl_cosmo, ptt_m, tracer2=ptt_m, ptc=ptc)
+    pk_mm = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_m, tracer2=ptt_m)
 
     # Set the parameters for our systematics
     systematics_params = ParamsMap(

examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py

@@ -61,10 +61,10 @@ def build_likelihood(_) -> Tuple[Likelihood, ModelingTools]:
                 )
 
     # Create the likelihood from the statistics
-    pt_calculator = pyccl.nl_pt.PTCalculator(
+    pt_calculator = pyccl.nl_pt.EulerianPTCalculator(
         with_NC=False,
         with_IA=True,
-        with_dd=True,
+        # with_dd=True,
         log10k_min=-4,
         log10k_max=2,
         nk_per_decade=20,
@@ -115,16 +115,21 @@ def run_likelihood() -> None:
     )
 
     # Code that creates a Pk2D object:
-    ptc = pyccl.nl_pt.PTCalculator(
-        with_NC=True, with_IA=True, log10k_min=-4, log10k_max=2, nk_per_decade=20
+    ptc = pyccl.nl_pt.EulerianPTCalculator(
+        with_NC=True,
+        with_IA=True,
+        log10k_min=-4,
+        log10k_max=2,
+        nk_per_decade=20,
+        cosmo=ccl_cosmo,
     )
     ptt_i = pyccl.nl_pt.PTIntrinsicAlignmentTracer(
         c1=(z, c_1), c2=(z, c_2), cdelta=(z, c_d)
     )
     ptt_m = pyccl.nl_pt.PTMatterTracer()
     # IAs x matter
-    pk_im = pyccl.nl_pt.get_pt_pk2d(ccl_cosmo, ptt_i, tracer2=ptt_m, ptc=ptc)
-    pk_ii = pyccl.nl_pt.get_pt_pk2d(ccl_cosmo, ptt_i, ptc=ptc)
+    pk_im = ptc.get_biased_pk2d(tracer1=ptt_i, tracer2=ptt_m)
+    pk_ii = ptc.get_biased_pk2d(tracer1=ptt_i, tracer2=ptt_i)
 
     # Set the parameters for our systematics
     systematics_params = ParamsMap(

examples/des_y1_3x2pt/numcosmo_firecrown_model.yml

@@ -20,15 +20,6 @@ scalar_params:
     absolute_tolerance: 0.0
     default_value: 0.0
     fit_type: !NcmParamType 'free'
-- !ScalarParameter
-    symbol: \alpha_g
-    name: alphag
-    lower_bound: -5.0
-    upper_bound: 5.0
-    scale: 0.1
-    absolute_tolerance: 0.0
-    default_value: -1.0
-    fit_type: !NcmParamType 'fixed'
 - !ScalarParameter
     symbol: z_\mathrm{piv}
     name: z_piv

firecrown/connector/mapping.py

@@ -192,7 +192,7 @@ def asdict(self) -> Dict[str, Union[Optional[float], List[float]]]:
             "Omega_g": self.Omega_g,
             "Neff": self.Neff,
             "m_nu": self.m_nu,
-            "m_nu_type": self.m_nu_type,
+            "mass_split": self.m_nu_type,
             "w0": self.w0,
             "wa": self.wa,
             "T_CMB": self.T_CMB,

firecrown/likelihood/gauss_family/statistic/two_point.py

@@ -444,13 +444,9 @@ def calculate_pk(
             # Compute perturbation power spectrum
 
             pt_calculator = tools.get_pt_calculator()
-            pk = pyccl.nl_pt.get_pt_pk2d(
-                tools.get_ccl_cosmology(),
-                tracer0.pt_tracer,
+            pk = pt_calculator.get_biased_pk2d(
+                tracer1=tracer0.pt_tracer,
                 tracer2=tracer1.pt_tracer,
-                nonlin_pk_type="nonlinear",
-                ptc=pt_calculator,
-                update_ptc=False,
             )
         elif tracer0.has_hm or tracer1.has_hm:
             # Compute halo model power spectrum

firecrown/modeling_tools.py

@@ -16,10 +16,10 @@ class ModelingTools:
     def __init__(
         self,
         *,
-        pt_calculator: Optional[pyccl.nl_pt.PTCalculator] = None,
+        pt_calculator: Optional[pyccl.nl_pt.EulerianPTCalculator] = None,
     ):
         self.ccl_cosmo: Optional[pyccl.Cosmology] = None
-        self.pt_calculator: Optional[pyccl.nl_pt.PTCalculator] = pt_calculator
+        self.pt_calculator: Optional[pyccl.nl_pt.EulerianPTCalculator] = pt_calculator
         self.powerspectra: Dict[str, pyccl.Pk2D] = {}
 
     def add_pk(self, name: str, powerspectrum: pyccl.Pk2D):
@@ -68,8 +68,7 @@ def prepare(self, ccl_cosmo: pyccl.Cosmology) -> None:
         self.ccl_cosmo = ccl_cosmo
 
         if self.pt_calculator is not None:
-            pk_lin_z0 = pyccl.linear_matter_power(ccl_cosmo, self.pt_calculator.ks, 1.0)
-            self.pt_calculator.update_pk(pk_lin_z0)
+            self.pt_calculator.update_ingredients(ccl_cosmo)
 
     @final
     def reset(self) -> None:
@@ -84,7 +83,7 @@ def get_ccl_cosmology(self) -> pyccl.Cosmology:
             raise RuntimeError("Cosmology has not been set")
         return self.ccl_cosmo
 
-    def get_pt_calculator(self) -> pyccl.nl_pt.PTCalculator:
+    def get_pt_calculator(self) -> pyccl.nl_pt.EulerianPTCalculator:
         """Return the perturbation theory calculator object."""
 
         if self.pt_calculator is None:

firecrown/models/cluster_abundance.py

@@ -117,8 +117,8 @@ def mf_d2N_dV_dlnM(self, ccl_cosmo: ccl.Cosmology, logM: float, z: float) -> flo
         if self.halo_mass_function is None:
             self.halo_mass_function = ccl.halos.MassFunc.from_name(
                 self.halo_mass_function_name
-            )(ccl_cosmo, **self.halo_mass_function_args)
-        nm = self.halo_mass_function.get_mass_function(ccl_cosmo, mass, a)
+            )(**self.halo_mass_function_args)
+        nm = self.halo_mass_function(ccl_cosmo, mass, a)
         return nm
 
     def mf_d2N_dz_dlnM(self, ccl_cosmo: ccl.Cosmology, logM: float, z: float) -> float:

tests/test_cluster.py

@@ -22,7 +22,7 @@ def fixture_ccl_cosmo():
     """Fixture for a CCL cosmology object."""
 
     return ccl.Cosmology(
-        Omega_c=0.22, Omega_b=0.0448, h=0.71, sigma8=0.8, n_s=0.963, Neff=3.04
+        Omega_c=0.22, Omega_b=0.0448, h=0.71, sigma8=0.8, n_s=0.963, Neff=3.44
     )
 
 
@@ -90,7 +90,7 @@ def fixture_cluster_mass_rich_args(parameters):
 def fixture_cluster_abundance(parameters):
     """Fixture for cluster objects."""
 
-    hmd_200 = ccl.halos.MassDef200c()
+    hmd_200 = ccl.halos.MassDef200c
     hmf_args: Dict[str, Any] = {}
     hmf_name = "Bocquet16"
     sky_area = 489

tests/test_modeling_tools.py

@@ -0,0 +1,35 @@
+"""
+Tests for the module firecrown.modeling_tools
+"""
+import pytest
+import pyccl
+from firecrown.modeling_tools import ModelingTools
+
+
+@pytest.fixture(name="dummy_powerspectrum")
+def make_dummy_powerspectrum() -> pyccl.Pk2D:
+    """Create an empty power spectrum. This is the only type we can create
+    without supplying a cosmology."""
+    return pyccl.Pk2D.__new__(pyccl.Pk2D)
+
+
+def test_default_constructed_state():
+    tools = ModelingTools()
+    # Default constructed state is pretty barren...
+    assert tools.ccl_cosmo is None
+    assert tools.pt_calculator is None
+    assert len(tools.powerspectra) == 0
+
+
+def test_default_constructed_no_tools():
+    tools = ModelingTools()
+    with pytest.raises(RuntimeError):
+        _ = tools.get_pk("nonesuch")
+
+
+def test_no_adding_pk_twice(dummy_powerspectrum: pyccl.Pk2D):
+    tools = ModelingTools()
+    tools.add_pk("silly", dummy_powerspectrum)
+    assert tools.powerspectra["silly"] == dummy_powerspectrum
+    with pytest.raises(KeyError):
+        tools.add_pk("silly", dummy_powerspectrum)

tests/test_pt_systematics.py

@@ -61,15 +61,6 @@ def test_pt_systematics(weak_lensing_source, number_counts_source, sacc_data):
         TwoPoint("galaxy_density_xi", number_counts_source, number_counts_source),
     ]
 
-    pt_calculator = pt.PTCalculator(
-        with_NC=True,
-        with_IA=True,
-        with_dd=False,
-        log10k_min=-4,
-        log10k_max=2,
-        nk_per_decade=4,
-    )
-
     likelihood = ConstGaussian(statistics=stats)
     likelihood.read(sacc_data)
     src0_tracer = sacc_data.get_tracer("src0")
@@ -81,6 +72,14 @@ def test_pt_systematics(weak_lensing_source, number_counts_source, sacc_data):
     ccl_cosmo = ccl.CosmologyVanillaLCDM()
     ccl_cosmo.compute_nonlin_power()
 
+    pt_calculator = pt.EulerianPTCalculator(
+        with_NC=True,
+        with_IA=True,
+        log10k_min=-4,
+        log10k_max=2,
+        nk_per_decade=4,
+        cosmo=ccl_cosmo,
+    )
     modeling_tools = ModelingTools(pt_calculator=pt_calculator)
     modeling_tools.prepare(ccl_cosmo)
 
@@ -99,20 +98,17 @@ def test_pt_systematics(weak_lensing_source, number_counts_source, sacc_data):
         ccl_cosmo, z=z, a1=a_1, a1delta=a_d, a2=a_2, Om_m2_for_c2=False
     )
 
-    # Code that creates a Pk2D object:
-    ptc = pt.PTCalculator(
-        with_NC=True, with_IA=True, log10k_min=-4, log10k_max=2, nk_per_decade=4
-    )
+    # Code that creates Pk2D objects:
     ptt_i = pt.PTIntrinsicAlignmentTracer(c1=(z, c_1), c2=(z, c_2), cdelta=(z, c_d))
     ptt_m = pt.PTMatterTracer()
     ptt_g = pt.PTNumberCountsTracer(b1=b_1, b2=b_2, bs=b_s)
     # IA
-    pk_im = pt.get_pt_pk2d(ccl_cosmo, ptt_i, tracer2=ptt_m, ptc=ptc)
-    pk_ii = pt.get_pt_pk2d(ccl_cosmo, ptt_i, ptc=ptc)
-    pk_gi = pt.get_pt_pk2d(ccl_cosmo, ptt_g, tracer2=ptt_i, ptc=ptc)
+    pk_im = pt_calculator.get_biased_pk2d(tracer1=ptt_i, tracer2=ptt_m)
+    pk_ii = pt_calculator.get_biased_pk2d(tracer1=ptt_i, tracer2=ptt_i)
+    pk_gi = pt_calculator.get_biased_pk2d(tracer1=ptt_g, tracer2=ptt_i)
     # Galaxies
-    pk_gm = pt.get_pt_pk2d(ccl_cosmo, ptt_g, tracer2=ptt_m, ptc=ptc)
-    pk_gg = pt.get_pt_pk2d(ccl_cosmo, ptt_g, ptc=ptc)
+    pk_gm = pt_calculator.get_biased_pk2d(tracer1=ptt_g, tracer2=ptt_m)
+    pk_gg = pt_calculator.get_biased_pk2d(tracer1=ptt_g, tracer2=ptt_g)
 
     # Set the parameters for our systematics
     systematics_params = ParamsMap(
@@ -234,15 +230,6 @@ def test_pt_mixed_systematics(sacc_data):
     )
 
     # Create the likelihood from the statistics
-    pt_calculator = pt.PTCalculator(
-        with_NC=True,
-        with_IA=True,
-        with_dd=False,
-        log10k_min=-4,
-        log10k_max=2,
-        nk_per_decade=4,
-    )
-
     likelihood = ConstGaussian(statistics=[stat])
     likelihood.read(sacc_data)
 
@@ -255,6 +242,14 @@ def test_pt_mixed_systematics(sacc_data):
     ccl_cosmo = ccl.CosmologyVanillaLCDM()
     ccl_cosmo.compute_nonlin_power()
 
+    pt_calculator = pt.EulerianPTCalculator(
+        with_NC=True,
+        with_IA=True,
+        log10k_min=-4,
+        log10k_max=2,
+        nk_per_decade=4,
+        cosmo=ccl_cosmo,
+    )
     modeling_tools = ModelingTools(pt_calculator=pt_calculator)
     modeling_tools.prepare(ccl_cosmo)
 
@@ -270,14 +265,11 @@ def test_pt_mixed_systematics(sacc_data):
         ccl_cosmo, z=z, a1=a_1, a1delta=a_d, a2=a_2, Om_m2_for_c2=False
     )
 
-    # Code that creates a Pk2D object:
-    ptc = pt.PTCalculator(
-        with_NC=True, with_IA=True, log10k_min=-4, log10k_max=2, nk_per_decade=4
-    )
+    # Code that creates Pk2D objects:
     ptt_i = pt.PTIntrinsicAlignmentTracer(c1=(z, c_1), c2=(z, c_2), cdelta=(z, c_d))
     ptt_m = pt.PTMatterTracer()
     # IA
-    pk_mi = pt.get_pt_pk2d(ccl_cosmo, ptt_m, tracer2=ptt_i, ptc=ptc)
+    pk_mi = pt_calculator.get_biased_pk2d(tracer1=ptt_m, tracer2=ptt_i)
 
     # Set the parameters for our systematics
     systematics_params = ParamsMap(