Baryons objects as convenience input for Cosmology (#1113)

* implemented, documented, and tested

* one more test, faster tests

* lint

* removed bcm option

pyccl/cosmology.py

@@ -17,7 +17,7 @@
 from . import (
     CCLError, CCLObject, CCLParameters, CosmologyParams,
     DEFAULT_POWER_SPECTRUM, DefaultParams, Pk2D, check, lib,
-    unlock_instance, emulators)
+    unlock_instance, emulators, baryons)
 from . import physical_constants as const
 
 
@@ -174,6 +174,9 @@ class Cosmology(CCLObject):
             The transfer function to use. Defaults to 'boltzmann_camb'.
         matter_power_spectrum (:obj:`str` or :class:`~pyccl.emulators.emu_base.EmulatorPk`):
             The matter power spectrum to use. Defaults to 'halofit'.
+        baryonic_effects (:class:`~pyccl.baryons.baryons_base.Baryons` or `None`):
+            The baryonic effects model to use. Options are `None` (no baryonic effects), or
+            a :class:`~pyccl.baryons.baryons_base.Baryons` object.
         extra_parameters (:obj:`dict`): Dictionary holding extra
             parameters. Currently supports extra parameters for CAMB.
             Details described below. Defaults to None.
@@ -195,10 +198,15 @@ class Cosmology(CCLObject):
 
         extra_parameters = {"camb": {"halofit_version": "mead2020_feedback",
                                      "HMCode_logT_AGN": 7.8}}
+
+    .. note :: If using camb to compute the non-linear power spectrum with HMCode
+               to include baryonic effects, you should not include any extra
+               baryonic effects (i.e. set `baryonic_effects=None`).
     """ # noqa
     from ._core.repr_ import build_string_Cosmology as __repr__
     __eq_attrs__ = ("_params_init_kwargs", "_config_init_kwargs",
-                    "_accuracy_params",)
+                    "_accuracy_params", "lin_pk_emu", 'nl_pk_emu',
+                    "baryons",)
 
     def __init__(
             self, *, Omega_c=None, Omega_b=None, h=None, n_s=None,
@@ -208,6 +216,7 @@ def __init__(
             mu_0=0, sigma_0=0, c1_mg=1, c2_mg=1, lambda_mg=0,
             transfer_function='boltzmann_camb',
             matter_power_spectrum='halofit',
+            baryonic_effects=None,
             extra_parameters=None,
             T_ncdm=DefaultParams.T_ncdm):
 
@@ -228,6 +237,12 @@ def __init__(
             self.nl_pk_emu = matter_power_spectrum
             matter_power_spectrum = 'emulator'
 
+        self.baryons = baryonic_effects
+        if not isinstance(self.baryons, baryons.Baryons):
+            if self.baryons is not None:
+                raise ValueError("`baryonic_effects` must be `None` "
+                                 "or a `Baryons` instance.")
+
         # going to save these for later
         self._params_init_kwargs = dict(
             Omega_c=Omega_c, Omega_b=Omega_b, h=h, n_s=n_s, sigma8=sigma8,
@@ -562,6 +577,9 @@ def _compute_nonlin_power(self):
         elif mps == 'emulator':
             pk = self.nl_pk_emu.get_pk2d(self)
 
+        # Include baryonic effects
+        if self.baryons is not None:
+            pk = self.baryons.include_baryonic_effects(self, pk)
         return pk
 
     @unlock_instance(mutate=False)

pyccl/tests/test_baryons.py

@@ -1,9 +1,9 @@
-import pytest
 import numpy as np
 import pyccl as ccl
+import pytest
 
 
-COSMO = ccl.CosmologyVanillaLCDM()
+COSMO = ccl.CosmologyVanillaLCDM(transfer_function='bbks')
 bar = ccl.BaryonsSchneider15()
 
 
@@ -41,3 +41,27 @@ def test_baryons_from_name():
     bar2 = ccl.Baryons.from_name('Schneider15')
     assert bar.name == bar2.name
     assert bar2.name == 'Schneider15'
+
+
+def test_baryons_in_cosmology():
+    # Test that applying baryons during cosmology creation works.
+    # 1. Outside of cosmo
+    cosmo_nb = ccl.CosmologyVanillaLCDM(
+        transfer_function='bbks', baryonic_effects=None)
+    pk_nb = cosmo_nb.get_nonlin_power()
+    pk_wb = bar.include_baryonic_effects(cosmo_nb, pk_nb)
+    # 2. In cosmo - from object.
+    cosmo_wb2 = ccl.CosmologyVanillaLCDM(
+        transfer_function='bbks', baryonic_effects=bar)
+    pk_wb2 = cosmo_wb2.get_nonlin_power()
+
+    ks = np.geomspace(1E-2, 10, 128)
+    pk_wb = pk_wb(ks, 1.0)
+    pk_wb2 = pk_wb2(ks, 1.0)
+
+    assert np.allclose(pk_wb, pk_wb2, atol=0, rtol=1E-6)
+
+
+def test_baryons_in_cosmology_error():
+    with pytest.raises(ValueError):
+        ccl.CosmologyVanillaLCDM(baryonic_effects=3.1416)