Allow sigma8 as an input when using CAMB for nonlinear power spectra (#…

…1106)

* Allow sigma8 as input with nonlinear power spectrum from CAMB.

* pacify flake8

* try fix gha

* try fix gha

* try fix gha

* try fix gha

* try fix gha

* Update environment.yml

* Update environment.yml

* testing abacus solution

* Update ci.yml

* testing abacus solution

* pip -v

* Install classy from source

* Rescale sigma8 when using MG

* Loosen test accuracies to be consistent with internal precision settings.

* More benchmark fixing.

* Use CCL sigma8 instead of CAMB sigma8

---------

Co-authored-by: David Alonso <dam.phys@gmail.com>
Co-authored-by: Hsin Fan <57552401+hsinfan1996@users.noreply.github.com>

benchmarks/test_ptpk_lpt.py

@@ -52,4 +52,6 @@ def test_pt_pk(comb):
         k = kin[ind]
         dpk = data[iz][i_d+1][ind]
         tpk = pk(k, a, COSMO)
-        assert np.all(np.fabs(tpk / dpk - 1) < LPTPK_TOLERANCE)
+        assert np.all(
+            np.fabs(tpk / dpk - 1) < LPTPK_TOLERANCE
+        )

pyccl/boltzmann.py

@@ -7,7 +7,7 @@
 except ModuleNotFoundError:
     pass  # prevent nans from isitgr
 
-from . import CCLError, Pk2D, check, lib
+from . import CCLError, Pk2D, check, lib, sigma8
 
 
 def get_camb_pk_lin(cosmo, *, nonlin=False):
@@ -48,9 +48,8 @@ def get_camb_pk_lin(cosmo, *, nonlin=False):
     if np.isfinite(cosmo["A_s"]):
         A_s_fid = cosmo["A_s"]
     elif np.isfinite(cosmo["sigma8"]):
-        # in this case, CCL will internally normalize for us when we init
-        # the linear power spectrum - so we just get close
-        A_s_fid = 2.43e-9 * (cosmo["sigma8"] / 0.87659)**2
+        A_s_fid = 2.1e-9
+        sigma8_target = cosmo["sigma8"]
     else:
         raise CCLError(
             "Could not normalize the linear power spectrum. "
@@ -140,61 +139,20 @@ def get_camb_pk_lin(cosmo, *, nonlin=False):
         wa=cosmo['wa']
     )
 
-    if nonlin:
-        cp.NonLinearModel = camb.nonlinear.Halofit()
-        halofit_version = extra_camb_params.get("halofit_version", "mead")
-        options = {k: extra_camb_params[k] for k in
-                   ["HMCode_A_baryon",
-                    "HMCode_eta_baryon",
-                    "HMCode_logT_AGN"] if k in extra_camb_params}
-        cp.NonLinearModel.set_params(halofit_version=halofit_version,
-                                     **options)
-
-    cp.set_matter_power(
-        redshifts=[_z for _z in zs],
-        kmax=extra_camb_params.get("kmax", 10.0),
-        nonlinear=nonlin)
-    if not nonlin:
-        assert cp.NonLinear == camb.model.NonLinear_none
-
     cp.set_for_lmax(extra_camb_params.get("lmax", 5000))
     cp.InitPower.set_params(
         As=A_s_fid,
         ns=cosmo['n_s'])
 
-    # run CAMB and get results
-    camb_res = camb.get_results(cp)
-    k, z, pk = camb_res.get_linear_matter_power_spectrum(
-        hubble_units=True, nonlinear=False)
-
-    # convert to non-h inverse units
-    k *= cosmo['h']
-    pk /= (h2 * cosmo['h'])
-
-    # now build interpolant
-    nk = k.shape[0]
-    lk_arr = np.log(k)
-    a_arr = 1.0 / (1.0 + z)
-    na = a_arr.shape[0]
-    sinds = np.argsort(a_arr)
-    a_arr = a_arr[sinds]
-    ln_p_k_and_z = np.zeros((na, nk), dtype=np.float64)
-    for i, sind in enumerate(sinds):
-        ln_p_k_and_z[i, :] = np.log(pk[sind, :])
+    cp.set_matter_power(
+        redshifts=[_z for _z in zs],
+        kmax=extra_camb_params.get("kmax", 10.0))
 
-    pk_lin = Pk2D(
-        a_arr=a_arr,
-        lk_arr=lk_arr,
-        pk_arr=ln_p_k_and_z,
-        is_logp=True,
-        extrap_order_lok=1,
-        extrap_order_hik=2)
+    camb_res = camb.get_transfer_functions(cp)
 
-    if not nonlin:
-        return pk_lin
-    else:
+    def construct_Pk2D(camb_res, nonlin=False):
         k, z, pk = camb_res.get_linear_matter_power_spectrum(
-            hubble_units=True, nonlinear=True)
+            hubble_units=True, nonlinear=nonlin)
 
         # convert to non-h inverse units
         k *= cosmo['h']
@@ -211,14 +169,44 @@ def get_camb_pk_lin(cosmo, *, nonlin=False):
         for i, sind in enumerate(sinds):
             ln_p_k_and_z[i, :] = np.log(pk[sind, :])
 
-        pk_nonlin = Pk2D(
+        pk = Pk2D(
             a_arr=a_arr,
             lk_arr=lk_arr,
             pk_arr=ln_p_k_and_z,
             is_logp=True,
             extrap_order_lok=1,
             extrap_order_hik=2)
 
+        return pk
+
+    if np.isfinite(cosmo["sigma8"]):
+        camb_res.calc_power_spectra()
+        pk = construct_Pk2D(camb_res, nonlin=False)
+        sigma8_tmp = sigma8(cosmo, p_of_k_a=pk)
+        camb_res.Params.InitPower.As *= sigma8_target**2 / sigma8_tmp**2
+
+    if nonlin:
+        camb_res.Params.NonLinear = camb.model.NonLinear_pk
+        camb_res.Params.NonLinearModel = camb.nonlinear.Halofit()
+        halofit_version = extra_camb_params.get("halofit_version", "mead")
+        options = {k: extra_camb_params[k] for k in
+                   ["HMCode_A_baryon",
+                    "HMCode_eta_baryon",
+                    "HMCode_logT_AGN"] if k in extra_camb_params}
+        camb_res.Params.NonLinearModel.set_params(
+            halofit_version=halofit_version,
+            **options)
+    else:
+        assert camb_res.Params.NonLinear == camb.model.NonLinear_none
+
+    camb_res.calc_power_spectra()
+    pk_lin = construct_Pk2D(camb_res, nonlin=False)
+
+    if not nonlin:
+        return pk_lin
+    else:
+        pk_nonlin = construct_Pk2D(camb_res, nonlin=True)
+
         return pk_lin, pk_nonlin
 
 

pyccl/cosmology.py

@@ -348,7 +348,7 @@ def _build_parameters(
 
         # Check to make sure specified amplitude parameter is consistent.
         if [A_s, sigma8].count(np.nan) != 1:
-            raise ValueError("Set either A_s or sigma8 and not both.")
+            raise ValueError("Set either A_s or sigma8 but not both.")
 
         # Check if any compulsory parameters are not set.
         compul = {"Omega_c": Omega_c, "Omega_b": Omega_b, "h": h, "n_s": n_s}
@@ -485,7 +485,14 @@ def _compute_linear_power(self):
         pk = None
         rescale_s8 = True
         rescale_mg = True
-        if trf == 'boltzmann_class':
+        if trf == "boltzmann_camb":
+            rescale_s8 = False
+            # For MG, the input sigma8 includes the effects of MG, while the
+            # sigma8 that CAMB uses is the GR definition. So we need to rescale
+            # sigma8 afterwards.
+            if self["mu_0"] != 0:
+                rescale_s8 = True
+        elif trf == 'boltzmann_class':
             pk = self.get_class_pk_lin()
         elif trf == 'boltzmann_isitgr':
             rescale_mg = False
@@ -505,11 +512,10 @@ def _compute_linear_power(self):
         # status variable to use it later if the transfer function is CAMB too.
         pkl = None
         if self._config_init_kwargs["matter_power_spectrum"] == "camb":
-            if not np.isfinite(self["A_s"]):
-                raise CCLError("CAMB doesn't rescale non-linear power spectra "
-                               "consistently without A_s.")
-            # no rescaling because A_s is necessarily provided
-            rescale_mg = rescale_s8 = False
+            rescale_mg = False
+            if self["mu_0"] != 0:
+                raise ValueError("Can't rescale non-linear power spectrum "
+                                 "from CAMB for mu-Sigma MG.")
             name = "delta_matter:delta_matter"
             pkl, self._pk_nl[name] = self.get_camb_pk_lin(nonlin=True)
 

pyccl/tests/test_cosmology.py

@@ -295,3 +295,12 @@ def test_ccl_physical_constants_smoke():
 def test_ccl_global_parameters_repr():
     ccl.spline_params.reload()
     assert eval(repr(ccl.spline_params)) == ccl.spline_params._bak
+
+
+def test_camb_sigma8_input():
+    sigma8 = 0.85
+    cosmo = ccl.Cosmology(
+        Omega_c=0.25, Omega_b=0.05, h=0.67, n_s=0.96, sigma8=sigma8,
+        transfer_function="boltzmann_camb"
+    )
+    assert np.isclose(cosmo.sigma8(), sigma8)

pyccl/tests/test_nonlin_camb_power.py

@@ -53,24 +53,6 @@ def test_nonlin_camb_power():
         assert np.allclose(pk_camb, pk_nonlin_ccl, rtol=3e-5)
 
 
-def test_nonlin_camb_power_with_sigma8():
-    Omega_c = 0.25
-    Omega_b = 0.05
-    n_s = 0.97
-    h = 0.7
-
-    ccl_cosmo = ccl.Cosmology(Omega_c=Omega_c, Omega_b=Omega_b, h=h, m_nu=0.0,
-                              sigma8=0.8, n_s=n_s,
-                              transfer_function="boltzmann_camb",
-                              matter_power_spectrum="camb")
-
-    k = np.logspace(-3, 1, 10)
-
-    # Check that non-linear power spectrum isn't being used with sigma8
-    with pytest.raises(ccl.errors.CCLError):
-        ccl.nonlin_matter_power(ccl_cosmo, k, 1.0)
-
-
 def test_nonlin_camb_power_raises():
     # Test that it raises when (trf, mps) == (no camb, camb).
     with pytest.raises(ccl.CCLError):

pyccl/tests/test_power.py

@@ -99,7 +99,7 @@ def test_power_sigma8norm_norms_consistent(tf):
         transfer_function=tf)
 
     # make sure they come out the same-ish
-    assert np.allclose(ccl.sigma8(cosmo), ccl.sigma8(cosmo_s8))
+    assert np.allclose(sigma8, ccl.sigma8(cosmo_s8))
 
     # and that the power spectra look right
     a = 0.8

pyccl/tests/test_power_mu_sigma_sigma8norm.py

@@ -59,6 +59,7 @@ def test_power_mu_sigma_sigma8norm_norms_consistent(tf):
 
     # make sure they come out the same-ish
     assert np.allclose(ccl.sigma8(cosmo), ccl.sigma8(cosmo_s8))
+
     if tf != 'boltzmann_isitgr':
         # and that the power spectra look right
         a = 0.8
@@ -68,3 +69,22 @@ def test_power_mu_sigma_sigma8norm_norms_consistent(tf):
             ccl.linear_matter_power(cosmo, 1e-4, a) /
             ccl.linear_matter_power(cosmo_s8, 1e-4, a))
         assert np.allclose(pk_rat, gfac)
+
+
+def test_nonlin_camb_MG_error():
+    Omega_c = 0.25
+    Omega_b = 0.05
+    n_s = 0.97
+    h = 0.7
+
+    ccl_cosmo = ccl.Cosmology(Omega_c=Omega_c, Omega_b=Omega_b, h=h, m_nu=0.0,
+                              A_s=2.1e-9, n_s=n_s,
+                              transfer_function="boltzmann_camb",
+                              matter_power_spectrum="camb",
+                              mu_0=0.1, sigma_0=0.2)
+
+    k = np.logspace(-3, 1, 10)
+
+    # Check that non-linear power spectrum isn't being used with sigma8
+    with pytest.raises(ValueError):
+        ccl.nonlin_matter_power(ccl_cosmo, k, 1.0)