Request for boost/pk2d output in Baryons module

[nikosarcevic]

Hi, David, all

Is there a chance/is it possible to get the boost for HMCode and BACCO in the same way as we do it for BCM and VanDaalen in Baryon module?
It would make the FISK migration to ccl v3 much easier.
Basically, I need a pk2d object rather than the numpy array.

I need it as pk2d object because this is the format that is required to pass to p_of_a_k argument when constructing the angular cls.

for example, this is exactly what I need:

bcm = ccl.baryons.BaryonsSchneider15()  # baryonic correction model
pk_bcm_boost = bcm.include_baryonic_effects(
    cosmo=cosmology_vanilla,
    pk=cosmology_vanilla.get_nonlin_power()
)

or this

vd = ccl.baryons.BaryonsvanDaalen19()  # Van Daalen 2019
boost_pk_vd = vd.include_baryonic_effects(
   cosmo=cosmology_vanilla,
    pk=cosmology_vanilla.get_nonlin_power()
)

but NOT this

cosmology_hmcode = ccl.CosmologyVanillaLCDM(
   matter_power_spectrum="camb",
   extra_parameters={"camb": {"kmax": 20.0,
                               "halofit_version": "mead2020_feedback",
                               "HMCode_logT_AGN": 7.8}}
)
pk_hmcode = cosmology_hmcode.nonlin_matter_power(k, a)

please tell me this can be somehow done? even if you do not plan to implement it into ccl, can you point me not the direction where I can make mods so within my code I have pk2d for hmcode and bacco

[damonge]

I don't think there's an easy way to do this right now. It would be great to have HMCode implemented natively (and that would solve this problem). Much of the infrastructure is actually there already, but it needs to be put together, which takes time...

[tilmantroester]

Seconding what David said. Having HMCode in the baryon challenge would be nice.

[nikosarcevic]

hmmm, is there a way to pass boosted pk but as a numpy array (NOT as a Pk2D) into ccl.angular_cl()?

I am trying to see if I can somehow construct a Pk2D boost factor myself for baccoemu and hmcode the way it was done for Schneider and VanDaalen but it is not clear to me immediately how to do it straightforwardly.

[nikosarcevic]

@damonge @tilmantroester

so if I call include_baryonic_effects method on the BaccoemuBaryons like this

baccoemu = ccl.BaccoemuBaryons()
boosted_pk_baccoemu = baccoemu.include_baryonic_effects(    
    cosmo=cosmology_vanilla,
    pk=cosmology_vanilla.get_nonlin_power()
)

I get an error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 2
      1 baccoemu = ccl.BaccoemuBaryons()
----> 2 boosted_pk_baccoemu = baccoemu.include_baryonic_effects(    
      3     cosmo=cosmology_vanilla,
      4     pk=cosmology_vanilla.get_nonlin_power()
      5 )

File ~/anaconda3/envs/FISK_LSST/lib/python3.11/site-packages/pyccl/baryons/baryons_base.py:41, in Baryons.include_baryonic_effects(self, cosmo, pk)
     30 def include_baryonic_effects(self, cosmo, pk):
     31     """Apply baryonic effects to a given power spectrum.
     32 
     33     Args:
   (...)
     39         :obj:`~pyccl.pk2d.Pk2D` or :obj:`None`.
     40     """
---> 41     return self._include_baryonic_effects(cosmo, pk)

File ~/anaconda3/envs/FISK_LSST/lib/python3.11/site-packages/pyccl/baryons/baccoemu_baryons.py:190, in BaccoemuBaryons._include_baryonic_effects(self, cosmo, pk)
    188 a_arr, lk_arr, pk_arr = pk.get_spline_arrays()
    189 k_arr = np.exp(lk_arr)
--> 190 fka = self.boost_factor(cosmo, k_arr, a_arr)
    191 pk_arr *= fka
    193 if pk.psp.is_log:

File ~/anaconda3/envs/FISK_LSST/lib/python3.11/site-packages/pyccl/baryons/baccoemu_baryons.py:94, in BaccoemuBaryons.boost_factor(self, cosmo, k, a)
     81 """The baccoemu BCM model boost factor for baryons.
     82 
     83 Args:
   (...)
     90     the power spectrum.
     91 """ # noqa
     93 # Check a ranges
---> 94 self._check_a_range(a)
     96 # First create the dictionary passed to baccoemu
     97 # if a is an array, make sure all the other parameters passed to the
     98 # emulator have the same len
     99 if np.ndim(a) == 1:

File ~/anaconda3/envs/FISK_LSST/lib/python3.11/site-packages/pyccl/baryons/baccoemu_baryons.py:208, in BaccoemuBaryons._check_a_range(self, a)
    206     a_max = max(a)
    207 if a_min < self.a_min or a_max > self.a_max:
--> 208     raise ValueError(f"Requested scale factor outside the bounds of "
    209                      f"the emulator: {(a_min, a_max)} outside of "
    210                      f"{((self.a_min, self.a_max))}")

ValueError: Requested scale factor outside the bounds of the emulator: (0.01, 1.0) outside of (0.4, 1.0)

any ideas how to get around this?

[tilmantroester]

That's the issue I had mentioned, the Bacco range for a doesn't cover the default ranges in CCL. Options I see are setting the CCL spline ranges for the power spectrum to the BaccoEmu ranges, creating the Pk2D object manually from a grid, like in the notebook, or add the baryonified version of BaccoEmu as a CCL emulator for the non-linear power spectrum.

[nikosarcevic]

That's the issue I had mentioned, the Bacco range for a doesn't cover the default ranges in CCL. Options I see are setting the CCL spline ranges for the power spectrum to the BaccoEmu ranges, creating the Pk2D object manually from a grid, like in the notebook, or add the baryonified version of BaccoEmu as a CCL emulator for the non-linear power spectrum.

so option 1 is to use something like I made in the CCLX notebook:

cosmology_vanilla = ccl.CosmologyVanillaLCDM()  # cosmology object
k = np.logspace(-3, 1, 100)  # wavenumber
a = 1.  # scale factor a z=0
baccoemu = ccl.BaccoemuBaryons()
cosmology_baccoemu = ccl.CosmologyVanillaLCDM()
bacco_boost = baccoemu.boost_factor(cosmology_baccoemu, k, a)
pk_bacco = cosmology_baccoemu.nonlin_matter_power(k, a) * bacco_boost

so I would need to convert bacco_boost to a Pk2D. The thing is im not sure how to do it

another option is to use pk_bacco and pass it into ccl.Cosmology object? Not sure I understood that correctly?

[tilmantroester]

That's option 2.
In that case you'd turn pk_bacco (but for a range of a, not a single a as in your example) into a Pk2D object and pass that to angular_cl

[nikosarcevic]

That's option 2. In that case you'd turn pk_bacco (but for a range of a, not a single a as in your example) into a Pk2D object and pass that to angular_cl

ok, so just to confirm -- I just need to make a list of scale factors (all equal to 1) and make a 2D array with k (nk x na) grid and then use that to construct pk2d?

[tilmantroester]

No, the scale factors need to span the range required for the Limber integrals. Have a look at the _include_baryonic_effects method in the BaccoemuBaryons class.

[nikosarcevic]

ok I coded this, I get the Pk2D within the k and a bounds

def create_baccoemu_pk2d(cosmo, baccoemu, num_points=100):
    """
    Create a Pk2D object using BaccoemuBaryons for baryonic effects in CCL,
    with automatically generated k and a arrays.

    Args:
        cosmo (ccl.Cosmology): A CCL Cosmology object.
        baccoemu (BaccoemuBaryons): A BaccoemuBaryons object for baryonic effects.
        num_points (int): Number of points for both the k and a arrays.

    Returns:
        Pk2D: A Pk2D object representing the power spectrum with baryonic effects.
    """
    # Generate k array (logarithmically spaced within valid range)
    k_arr = np.logspace(np.log10(baccoemu.k_min), np.log10(baccoemu.k_max), num_points)

    # Generate a array (linearly spaced within valid range)
    a_arr = np.linspace(baccoemu.a_min, baccoemu.a_max, num_points)

    # Initialize the 2D array for the power spectrum
    pk_arr = np.zeros((num_points, num_points))

    # Compute the power spectrum for each scale factor
    for i, a in enumerate(a_arr):
        # Check if the scale factor is within the bounds
        baccoemu._check_a_range(a)

        # Compute the boost factor for baryonic effects
        bacco_boost = baccoemu.boost_factor(cosmo, k_arr, a)

        # Compute the matter power spectrum at this scale factor
        pk_lin = ccl.nonlin_matter_power(cosmo, k_arr, a)

        # Apply the baryonic boost factor
        pk_arr[i, :] = pk_lin * bacco_boost

    # Convert k_arr to logarithmic values since Pk2D expects log(k)
    lk_arr = np.log(k_arr)

    # Create the Pk2D object
    pk2d = ccl.Pk2D(a_arr=a_arr, lk_arr=lk_arr, pk_arr=pk_arr, is_logp=False)
    return pk2d

baccoemu = ccl.BaccoemuBaryons()
pk2d_baccoemu = create_baccoemu_pk2d(cosmology_vanilla, baccoemu)
pk2d_baccoemu

this works, but I would appreciate a sanity check from you two @tilmantroester @damonge