Skip to content

Commit

Permalink
Check a range on BaccoemuBaryons (#1108)
Browse files Browse the repository at this point in the history
* Check a range on BaccoemuBaryons

* Replace isinstance(Iterable) and hasattr(__len__) with np.ndim

* Added test

---------

Co-authored-by: David Alonso <dam.phys@gmail.com>
  • Loading branch information
tilmantroester and damonge committed Jul 21, 2023
1 parent 4ef1507 commit de8127b
Show file tree
Hide file tree
Showing 2 changed files with 33 additions and 6 deletions.
29 changes: 23 additions & 6 deletions pyccl/baryons/baccoemu_baryons.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,19 +31,22 @@ class BaccoemuBaryons(Baryons):
hot gas in haloes
log10_M_inn (:obj:`float`): transition mass of density profiles of
hot gas in haloes (in :math:`M_\\odot`)
verbose (:obj:`bool`): Verbose output from baccoemu.
(default: :obj:`False`)
"""
name = 'BaccoemuBaryons'
__repr_attrs__ = __eq_attrs__ = ("bcm_params",)

def __init__(self, log10_M_c=14.174, log10_eta=-0.3, log10_beta=-0.22,
log10_M1_z0_cen=10.674, log10_theta_out=0.25,
log10_theta_inn=-0.86, log10_M_inn=13.574):
log10_theta_inn=-0.86, log10_M_inn=13.574,
verbose=False):
# avoid tensorflow warnings
import warnings
with warnings.catch_warnings():
warnings.filterwarnings('ignore', category=UserWarning)
import baccoemu
self.mpk = baccoemu.Matter_powerspectrum()
self.mpk = baccoemu.Matter_powerspectrum(verbose=verbose)
self.a_min = self.mpk.emulator['baryon']['bounds'][-1][0]
self.a_max = self.mpk.emulator['baryon']['bounds'][-1][1]
self.k_min = self.mpk.emulator['baryon']['k'][0]
Expand All @@ -62,7 +65,7 @@ def __init__(self, log10_M_c=14.174, log10_eta=-0.3, log10_beta=-0.22,
def _sigma8tot_2_sigma8cold(self, emupars, sigma8tot):
"""Use baccoemu to convert sigma8 total matter to sigma8 cdm+baryons
"""
if hasattr(emupars['omega_cold'], '__len__'):
if np.ndim(emupars['omega_cold']) == 1:
_emupars = {}
for pname in emupars:
_emupars[pname] = emupars[pname][0]
Expand All @@ -87,10 +90,13 @@ def boost_factor(self, cosmo, k, a):
the power spectrum.
""" # noqa

# Check a ranges
self._check_a_range(a)

# First create the dictionary passed to baccoemu
# if a is an array, make sure all the other parameters passed to the
# emulator have the same len
if hasattr(a, '__len__'):
if np.ndim(a) == 1:
emupars = dict(
omega_cold=np.full((len(a)),
cosmo['Omega_c'] + cosmo['Omega_b']),
Expand Down Expand Up @@ -122,12 +128,12 @@ def boost_factor(self, cosmo, k, a):
# power spectrum; so we have to convert from total to cold sigma8
sigma8tot = cosmo['sigma8']
sigma8cold = self._sigma8tot_2_sigma8cold(emupars, sigma8tot)
if hasattr(a, '__len__'):
if np.ndim(a) == 1:
emupars['sigma8_cold'] = np.full((len(a)), sigma8cold)
else:
emupars['sigma8_cold'] = sigma8cold
else:
if hasattr(a, '__len__'):
if np.ndim(a) == 1:
emupars['A_s'] = np.full((len(a)), cosmo['A_s'])
else:
emupars['A_s'] = cosmo['A_s']
Expand Down Expand Up @@ -176,3 +182,14 @@ def _include_baryonic_effects(self, cosmo, pk):
is_logp=pk.psp.is_log,
extrap_order_lok=pk.extrap_order_lok,
extrap_order_hik=pk.extrap_order_hik)

def _check_a_range(self, a):
if np.ndim(a) == 0:
a_min, a_max = a, a
else:
a_min = min(a)
a_max = max(a)
if a_min < self.a_min or a_max > self.a_max:
raise ValueError(f"Requested scale factor outside the bounds of "
f"the emulator: {(a_min, a_max)} outside of "
f"{((self.a_min, self.a_max))}")
10 changes: 10 additions & 0 deletions pyccl/tests/test_baccoemu.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
import numpy as np
import pyccl as ccl

import pytest

BEMULIN_TOLERANCE = 1e-3
BEMUNL_TOLERANCE = 5e-3
BEMBAR_TOLERANCE = 1e-3
Expand Down Expand Up @@ -114,6 +116,14 @@ def test_baccoemu_baryons_changepars():
& (baryons.bcm_params['eta'] == -0.4))


def test_baccoemu_baryons_a_range():
baryons = ccl.BaccoemuBaryons()
cosmo = ccl.CosmologyVanillaLCDM()
k = 1e-1
with pytest.raises(ValueError):
baryons.boost_factor(cosmo, k, baryons.a_min * 0.9)


def test_baccoemu_baryons_As_sigma8():
baryons = ccl.BaccoemuBaryons()
cosmo1 = ccl.Cosmology(
Expand Down

0 comments on commit de8127b

Please sign in to comment.