diff --git a/examples/bacco-values.ini b/examples/bacco-values.ini new file mode 100644 index 00000000..f9e82e5a --- /dev/null +++ b/examples/bacco-values.ini @@ -0,0 +1,24 @@ +[cosmological_parameters] +omega_m = 0.3 +h0 = 0.7 +omega_b = 0.048 +n_s = 0.97 +A_s = 2.19e-9 +w = -1.0 +wa = 0.0 +; New parametrization and names: +mnu = 0.0773 +num_massive_neutrinos = 3 +nnu = 3.046 + +omega_k = 0.0 +tau = 0.0697186 + +[baryon_parameters] +M_c = 9.0 14.0 15.0 +eta = -0.698 -0.3 0.698 +beta = -1.0 -0.22 0.698 +M1_z0_cen= 9.0 10.5 13.0 +theta_out= 0.0 0.25 0.477 +theta_inn = -2.0 -0.86 -0.522 +M_inn = 9.0 13.4 13.5 diff --git a/examples/bacco.ini b/examples/bacco.ini new file mode 100644 index 00000000..f3a74dae --- /dev/null +++ b/examples/bacco.ini @@ -0,0 +1,54 @@ +[runtime] +sampler = test +verbosity = debug + +[apriori] +nsample = 100 + +[output] +filename = output/bacco.txt + +[pipeline] +modules = consistency bbn_consistency camb extrapolate bacco_emulator + +timing=F +debug=F + +values = examples/bacco-values.ini +extra_output = + +[test] +save_dir=output/bacco +fatal_errors=T + +[consistency] +file = utility/consistency/consistency_interface.py + +[bbn_consistency] +file = utility/bbn_consistency/bbn_consistency.py + +[camb] +file = boltzmann/camb/camb_interface.py +mode = power +lmax = 2500 ;max ell to use for cmb calculation +feedback=0 ;amount of output to print +AccuracyBoost=1.1 ;CAMB accuracy boost parameter +do_tensors = T +do_lensing = T +NonLinear = none +zmin_background = 0. +zmax_background = 4. +nz_background = 401 +kmin=1e-4 +kmax = 50.0 +kmax_extrapolate = 500.0 +nk=700 +nz = 150 + +[extrapolate] +file = boltzmann/extrapolate/extrapolate_power.py +kmax = 500. + +[bacco_emulator] +file = structure/baccoemu/baccoemu_interface.py +mode = nonlinear diff --git a/examples/des-campaign.yml b/examples/des-campaign.yml index 79afc0a4..ef4fb03e 100644 --- a/examples/des-campaign.yml +++ b/examples/des-campaign.yml @@ -23,6 +23,35 @@ runs: params: - sampler = test + - name: bacco-nl + parent: fiducial + params: + - bacco.file = structure/baccoemu/baccoemu_interface.py + - bacco.mode = nonlinear + - camb.nonlinear = none + pipeline: + - after camb bacco + + - name: bacco-baryons + parent: fiducial + components: + - bacco_baryon_params + params: + - bacco.file = structure/baccoemu/baccoemu_interface.py + - bacco.mode = baryons + pipeline: + - after camb bacco + + - name: bacco-nl-baryons + parent: fiducial + components: + - bacco_baryon_params + params: + - bacco.file = structure/baccoemu/baccoemu_interface.py + - bacco.mode = nonlinear+baryons + pipeline: + - after camb bacco + # This run is based on the fiducial run above, with # additional changes applied to it on top of the ones above - name: class @@ -120,6 +149,15 @@ runs: # These components can be re-used in multiple runs above. components: + bacco_baryon_params: + values: + - baryon_parameters.M_c = 9.0 14.0 15.0 + - baryon_parameters.eta = -0.698 -0.3 0.698 + - baryon_parameters.beta = -1.0 -0.22 0.698 + - baryon_parameters.M1_z0_cen= 9.0 10.5 13.0 + - baryon_parameters.theta_out= 0.0 0.25 0.477 + - baryon_parameters.theta_inn = -2.0 -0.86 -0.522 + - baryon_parameters.M_inn = 9.0 13.4 13.5 maglim_cuts: params: - 2pt_like.angle_range_xip_1_1 = 2.475 999.0 diff --git a/likelihood/planck2018/.gitignore b/likelihood/planck2018/.gitignore index 2155bb53..e30ec294 100644 --- a/likelihood/planck2018/.gitignore +++ b/likelihood/planck2018/.gitignore @@ -2,3 +2,4 @@ data/ plc-3.0/buildir plc-3.0/include plc-3.0/share +baseline diff --git a/structure/baccoemu/LICENSE b/structure/baccoemu/LICENSE new file mode 100644 index 00000000..1227a037 --- /dev/null +++ b/structure/baccoemu/LICENSE @@ -0,0 +1,20 @@ +Copyright (c) 2020 The Python Packaging Authority + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + diff --git a/structure/baccoemu/baccoemu_interface.py b/structure/baccoemu/baccoemu_interface.py new file mode 100644 index 00000000..a7e86387 --- /dev/null +++ b/structure/baccoemu/baccoemu_interface.py @@ -0,0 +1,207 @@ +from cosmosis.datablock import option_section, names +from scipy.interpolate import RectBivariateSpline +import numpy as np +import baccoemu_vendored as baccoemu +import traceback + + +def setup(options): + mode = options.get_string(option_section, "mode", default="nonlinear") + + # do this only once in the setup phase + emulator = baccoemu.Matter_powerspectrum() + + allowed_modes = ["nonlinear", "baryons", "nonlinear+baryons"] + if mode not in allowed_modes: + raise ValueError(f"unrecognised value of 'mode' parameter in baccoemu: {mode}") + + return mode, emulator + + +def check_params(params): + ranges = { + "omega_cold": [0.15, 0.6], + "omega_baryon": [0.03, 0.07], + "ns": [0.92, 1.01], + "hubble": [0.6, 0.8], + "neutrino_mass": [0.0, 0.4], + } + + for key, (vmin, vmax) in ranges.items(): + if not (params[key] > vmin) and (params[key] < vmax): + raise ValueError(f"BaccoEmu: {key} must be in range [{vmin},{vmax}]") + + +def execute(block, config): + mode, emulator = config + + try: + if mode == "nonlinear": + emulate_nonlinear_power(block, emulator) + elif mode == "nonlinear+baryons": + emulate_boosted_nonlinear_power(block, emulator) + elif mode == "baryons": + emulate_baryonic_boost(block, emulator) + else: + raise RuntimeError("This should not happen") + except ValueError as error: + # print traceback from exception + print(traceback.format_exc()) + return 1 + + return 0 + + +def emulate_nonlinear_power(block, emulator): + z_lin, k_lin, P_lin = block.get_grid("matter_power_lin", "z", "k_h", "p_k") + + # bacco specific stuff + # This is required to avoid a bounds error + zmask = z_lin < 1.5 + a_lin = 1 / (1 + z_lin[zmask]) + + kmask = (k_lin < 4.69) & (k_lin > 0.0001) + cosmo = names.cosmological_parameters + + omega_cold = block[cosmo, "omega_c"] + block[cosmo, "omega_b"] + Omb = block[cosmo, "omega_b"] + + params = { + "omega_cold": omega_cold, + "A_s": block[cosmo, "a_s"], + "omega_baryon": Omb, + "ns": block[cosmo, "n_s"], + "hubble": block[cosmo, "h0"], + "neutrino_mass": block[cosmo, "mnu"], + "w0": block.get_double(cosmo, "w", -1.0), + "wa": 0.0, + "expfactor": a_lin, + } + + # check we're within the allowed bounds for the emulator + check_params(params) + + # evaluate the nonlinear growth factor as a function of k and z + k_bacco, F = emulator.get_nonlinear_boost(k=k_lin[kmask], cold=False, **params) + k_nl = k_lin + z_nl = z_lin + + # interpolate it to the same sampling as the linear matter power spectrum + I = RectBivariateSpline(np.log10(k_bacco), z_lin[zmask], np.log10(F).T) + F_interp = 10 ** I(np.log10(k_nl), z_nl).T + + # apply the factor + P_nl = F_interp * P_lin + + # save everything + block.put_grid("matter_power_nl", "z", z_nl, "k_h", k_nl, "p_k", P_nl) + + return 0 + + +def emulate_boosted_nonlinear_power(block, emulator): + z_lin, k_lin, P_lin = block.get_grid("matter_power_lin", "z", "k_h", "p_k") + + # This is required to avoid a bounds error + zmask = z_lin < 1.5 + a_lin = 1 / (1 + z_lin[zmask]) + + kmask = (k_lin < 4.69) & (k_lin > 0.0001) + + cosmo = names.cosmological_parameters + + # In BACCO omega_cold refers to baryons + CDM + omega_cold = block[cosmo, "omega_c"] + block[cosmo, "omega_b"] + + params = { + "omega_cold": omega_cold, + "A_s": block[cosmo, "a_s"], + "omega_baryon": block[cosmo, "omega_b"], + "ns": block[cosmo, "n_s"], + "hubble": block[cosmo, "h0"], + "neutrino_mass": block[cosmo, "mnu"], + "w0": block.get_double(cosmo, "w", -1.0), + "wa": block.get_double(cosmo, "wa", 0.0), + "expfactor": a_lin, + "M_c": block["baryon_parameters", "m_c"], + "eta": block["baryon_parameters", "eta"], + "beta": block["baryon_parameters", "beta"], + "M1_z0_cen": block["baryon_parameters", "m1_z0_cen"], + "theta_out": block["baryon_parameters", "theta_out"], + "theta_inn": block["baryon_parameters", "theta_inn"], + "M_inn": block["baryon_parameters", "m_inn"], + } + + # check we're within the allowed bounds for the emulator + check_params(params) + + k_nl = k_lin + z_nl = z_lin + + # evaluate the nonlinear growth factor as a function of k and z + k_bacco, F = emulator.get_nonlinear_boost(k=k_lin[kmask], cold=False, **params) + I_nl = RectBivariateSpline(np.log10(k_bacco), z_lin[zmask], np.log10(F).T) + F_interp = 10 ** I_nl(np.log10(k_nl), z_nl).T + + # same thing for baryons + k_bacco, S = emulator.get_baryonic_boost(k=k_lin[kmask], **params) + I_baryon = RectBivariateSpline(np.log10(k_bacco), z_lin[zmask], S.T) + S_interp = I_baryon(np.log10(k_nl), z_nl).T + + # apply the factor + P_nl = F_interp * S_interp * P_lin + + # save everything + block.put_grid("matter_power_nl", "z", z_nl, "k_h", k_nl, "p_k", P_nl) + + return 0 + + +def emulate_baryonic_boost(block, emulator): + z_nl, k_nl, P_nl = block.get_grid("matter_power_nl", "z", "k_h", "p_k") + + # This is required to avoid a bounds error + zmask = z_nl < 1.5 + a_lin = 1 / (1 + z_nl[zmask]) + + kmask = (k_nl < 4.69) & (k_nl > 0.0001) + + cosmo = names.cosmological_parameters + + # In BACCO omega_cold refers to baryons + CDM + omega_cold = block[cosmo, "omega_c"] + block[cosmo, "omega_b"] + + params = { + "omega_cold": omega_cold, + "A_s": block[cosmo, "a_s"], + "omega_baryon": block[cosmo, "omega_b"], + "ns": block[cosmo, "n_s"], + "hubble": block[cosmo, "h0"], + "neutrino_mass": block[cosmo, "mnu"], + "w0": block.get_double(cosmo, "w", -1.0), + "wa": block.get_double(cosmo, "wa", 0.0), + "expfactor": a_lin, + "M_c": block["baryon_parameters", "m_c"], + "eta": block["baryon_parameters", "eta"], + "beta": block["baryon_parameters", "beta"], + "M1_z0_cen": block["baryon_parameters", "m1_z0_cen"], + "theta_out": block["baryon_parameters", "theta_out"], + "theta_inn": block["baryon_parameters", "theta_inn"], + "M_inn": block["baryon_parameters", "m_inn"], + } + + # check we're within the allowed bounds for the emulator + check_params(params) + + # just get the boost factor from bacco + k_bacco, S = emulator.get_baryonic_boost(k=k_nl[kmask], **params) + I_baryon = RectBivariateSpline(np.log10(k_bacco), z_nl[zmask], S.T) + S_interp = I_baryon(np.log10(k_nl), z_nl).T + + # apply the factor + P_nl *= S_interp + + # save everything + block.put_grid("matter_power_nl", "z", z_nl, "k_h", k_nl, "p_k", P_nl) + + return 0 diff --git a/structure/baccoemu/baccoemu_vendored/.gitignore b/structure/baccoemu/baccoemu_vendored/.gitignore new file mode 100644 index 00000000..e4e943ec --- /dev/null +++ b/structure/baccoemu/baccoemu_vendored/.gitignore @@ -0,0 +1,6 @@ +baryon_emu_* +cold_matter_linear_emu_* +no_wiggles_emu_* +sigma8_emu_* +nonlinear_emu_* +total_matter_linear_emu_* diff --git a/structure/baccoemu/baccoemu_vendored/__init__.py b/structure/baccoemu/baccoemu_vendored/__init__.py new file mode 100644 index 00000000..18c66e88 --- /dev/null +++ b/structure/baccoemu/baccoemu_vendored/__init__.py @@ -0,0 +1,16 @@ +import numpy as np +import copy +import pickle +import tqdm +import hashlib +from ._version import __version__ +from .utils import * +from .matter_powerspectrum import * +from .baryonic_boost import * +from .lbias_expansion import * + +import tensorflow +from packaging import version +required_tf_version = '2.6.0' +if version.parse(tensorflow.__version__) < version.parse(required_tf_version): + raise ImportError(f'tensorflow={tensorflow.__version__} is not supported by baccoemu. Please update tensorflow >= 2.6.0') diff --git a/structure/baccoemu/baccoemu_vendored/_version.py b/structure/baccoemu/baccoemu_vendored/_version.py new file mode 100644 index 00000000..9aa3f903 --- /dev/null +++ b/structure/baccoemu/baccoemu_vendored/_version.py @@ -0,0 +1 @@ +__version__ = "2.1.0" diff --git a/structure/baccoemu/baccoemu_vendored/baryonic_boost.py b/structure/baccoemu/baccoemu_vendored/baryonic_boost.py new file mode 100644 index 00000000..3892f0de --- /dev/null +++ b/structure/baccoemu/baccoemu_vendored/baryonic_boost.py @@ -0,0 +1,271 @@ +import numpy as np +import copy +import pickle +from .utils import _transform_space, MyProgressBar + +__all__ = ["get_baryon_fractions"] + +def load_baryonic_emu(fold_name=None, detail_name=None, verbose=True): + """Loads in memory the baryonic boost emulator, described in Aricò et al. 2020c. + + :param verbose: whether to activate the verbose mode, defaults to True. + :type verbose: boolean, optional + + :return: a dictionary containing the emulator object + :rtype: dict + """ + import os + if verbose: + print('Loading Baryonic Emulator...') + basefold = os.path.dirname(os.path.abspath(__file__)) + + from tensorflow.keras.models import load_model + + old_emulator_names = [(basefold + '/' + + "NN_emulator_sg_0.99_15000_PCA6_BNFalse_DO0.0_NL2_data_midmfcorr_10k")] + for old_emulator_name in old_emulator_names: + if os.path.exists(old_emulator_name): + import shutil + shutil.rmtree(old_emulator_name) + + if fold_name is not None: + emulator_name = fold_name + else: + emulator_name = (basefold + '/' + + "baryon_emu_1.0.0") + + if (not os.path.exists(emulator_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading Emulator data (5 Mb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/baryon_emu_1.0.0.tar', + emulator_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_name + '.tar') + + emulator = {} + emulator['emu_type'] = 'nn' + emulator['model'] = load_model(emulator_name, compile=False) + detail_name = 'k_scaler_bounds.pkl' if detail_name is None else detail_name + with open(emulator_name + '/' +detail_name, 'rb') as f: + emulator['scaler'] = pickle.load(f) + emulator['pca'] = pickle.load(f) + emulator['k'] = pickle.load(f) + emulator['values'] = pickle.load(f) + if fold_name == None: + emulator['rotation'] = pickle.load(f) + emulator['bounds'] = pickle.load(f) + + if fold_name == None: + emulator['keys'] = ['omega_cold', 'sigma8_cold', 'omega_baryon', + 'ns', 'hubble', 'neutrino_mass', 'w0', 'wa', 'M_c', 'eta', 'beta', 'M1_z0_cen', + 'theta_inn', 'M_inn', 'theta_out', 'expfactor'] + else: + emulator['keys'] = ['omega_cold', 'sigma8_cold', 'omega_baryon', + 'hubble', 'M_c', 'eta', 'beta', 'M1_z0_cen', + 'theta_inn','M_inn','theta_out', 'expfactor'] + if verbose: + print('Baryonic Emulator loaded in memory.') + return emulator + +def get_baryon_fractions(M_200c, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, M_c=None, eta=None, beta=None, + M1_z0_cen=None, theta_out=None, theta_inn=None, M_inn=None): + """Compute the mass fraction of the different baryonic components, following the + baryonic correction model described in Aricò et al 2020b (see also Aricò et al 2020a). + + :param M_200: Halo mass inside a sphere which encompass a density 200 times larger than the critical density of the Universe. + :type array_like + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + :param M_c: mass fraction of hot gas in haloes + :type M_c: float or array + :param eta: extent of ejected gas + :type eta: float or array + :param beta: mass fraction of hot gas in haloes + :type beta: float or array + :param M1_z0_cen: characteristic halo mass scale for central galaxy + :type M1_z0_cen: float or array + :param theta_out: density profile of hot gas in haloes + :type theta_out: float or array + :param theta_inn: density profile of hot gas in haloes + :type theta_inn: float or array + :param M_inn: density profile of hot gas in haloes + :type M_inn: float or array + + :return: a dictionary containing the baryonic components mass fractions, with the following keywords: + + #. 'ej_gas' -> ejected gas + #. 'cen_galaxy' -> central galaxy + #. 'sat_galaxy' -> satellite galaxy + #. 'bo_gas' -> bound gas + #. 're_gas' -> reaccreted gas + #. 'dark_matter' -> dark matter + #. 'gas' -> total gas + #. 'stellar' -> total stars + #. 'baryon' -> total baryons + + :rtype: dict + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + coord = copy.deepcopy(kwargs) + for par in ['M_c','eta','beta','M1_z0_cen','theta_out','theta_inn','M_inn']: + coord[par] = 10**(coord[par]) + + stellar_pars = _SMHM_relation(coord) + frac_cg = _galaxy_fraction(stellar_pars, M_200c, 'centrals') + frac_sg = _galaxy_fraction(stellar_pars, M_200c, 'satellites') + frac_stars = frac_cg + frac_sg + + + o_m = coord['omega_matter'] if coord['omega_matter'] is not None else coord['omega_cold'] + coord['neutrino_mass'] / 93.14 / coord['hubble']**2 + + baryon_matter_ratio = coord['omega_baryon']/o_m + + if np.any(frac_stars > baryon_matter_ratio): + raise ValueError(""" Your stellar fraction is larger than the baryon fraction! + Please use meaningful stellar parameters""") + + frac_bg = (baryon_matter_ratio-frac_stars)/(1+(coord['M_c']/M_200c)**coord['beta']) + frac_bar = frac_bg + frac_stars + + M_r = 1e16 + beta_r = 2. + frac_re = (baryon_matter_ratio-frac_bar)/(1+(M_r/M_200c)**beta_r) + frac_bg -= frac_re #fraction of reaccreted gas is taken from the bound gas + frac_eg = baryon_matter_ratio - frac_bar + frac_dm = 1 - baryon_matter_ratio + + frac_gas = frac_bg + frac_eg + frac_re + return {'ej_gas': np.array(frac_eg, dtype=np.float32), + 'cen_galaxy': np.array(frac_cg, dtype=np.float32), + 'sat_galaxy': np.array(frac_sg, dtype=np.float32), + 'bo_gas': np.array(frac_bg, dtype=np.float32), + 're_gas': np.array(frac_re, dtype=np.float32), + 'dark_matter': np.array(frac_dm, dtype=np.float32), + 'gas': np.array(frac_gas, dtype=np.float32), + 'stellar': np.array(frac_stars, dtype=np.float32), + 'baryon': np.array(baryon_matter_ratio, dtype=np.float32), + } + +def _SMHM_relation(coordinates): + """ + Internal function which evolve the Stellar Mass to Halo Mass (SMHM) relation parameters + at the correct redshift, following Behroozi et al. 2013. + + :param coordinates: a set of coordinates in parameter space + :type coordinates: dict + :param a: expansion factor + :type a: float + :return: dictionary with the evolved SHAM parameters. + :rtype: dictionary + """ + a = coordinates['expfactor'] + z = 1/a -1 + nu = np.exp(-4*a**2) #Exponential cutoff of evolution of M ∗ (M h ) with scale factor + + pars = {} + pars['M1_z0_cen'] = np.float64(coordinates['M1_z0_cen']) + pars['epsilon_z0_cen'] = np.float32(0.023) + pars['alpha_z0_cen'] = np.float32(-1.779) + pars['gamma_z0_cen'] = np.float32(0.547) + pars['delta_z0_cen'] = np.float32(4.394) + + pars['M1_fsat'] = np.float32(1.59) + pars['epsilon_fsat'] = np.float32(0.2) + pars['alpha_fsat'] = np.float32(0.16) + pars['gamma_fsat'] = np.float32(1.67) + pars['delta_fsat'] = np.float32(0.99) + + ini = {'a':0,'a2':0,'z':0} + stellar_pars = { + 'M1': copy.deepcopy(ini), # Charactheristic halo mass Msolar/h + 'epsilon':copy.deepcopy(ini), # Characteristic stellar mass to halo mass ratio + 'alpha':copy.deepcopy(ini), # Faint-end slope of SMHM relation + 'delta':copy.deepcopy(ini), #Index of subpower law at massive end of SMHM relation + 'gamma':copy.deepcopy(ini) #Strength of subpower law at massive end of SMHM relation + } + + stellar_pars['M1']['a'] = -1.793; stellar_pars['M1']['z'] = -0.251 + stellar_pars['alpha']['a'] = 0.731; + stellar_pars['gamma']['a'] = 1.319; stellar_pars['gamma']['z'] = 0.279 + stellar_pars['delta']['a'] = 2.608; stellar_pars['delta']['a'] = -0.043 + stellar_pars['epsilon']['a'] = -0.006; stellar_pars['epsilon']['a2'] = -0.119; + + #for central galaxies: + for p in stellar_pars.keys(): + p0_cen = np.log10(pars[p+'_z0_cen']) if p in ['M1','epsilon'] else pars[p+'_z0_cen'] + pcen_z = p0_cen + nu*(stellar_pars[p]['a']*(a-1) + stellar_pars[p]['z']*z) + stellar_pars[p]['a2']*(a-1) + pars[p+'_cen'] = 10**(pcen_z) if p in ['M1','epsilon'] else pcen_z + + #for satellites: + for p in stellar_pars.keys(): + pars[p+'_sat'] = pars[p+'_cen'] * pars[p+'_fsat'] + return pars + +def _galaxy_fraction(pars, M_200, type='centrals'): + """ + Function which compute the galaxy mass fractions. + :param pars: set of SHAM parameters + :type pars: dict + :param M_200: Halo mass inside a sphere which encompass a density 200 times larger than the critical density of the Universe. + :type array_like + :param type: 'centrals' to select central galaxies, 'satellites' to select satellites galaxies. + :type string + :return: galaxy mass fractions. + :rtype: array_like + """ + t = '_cen' if type == 'centrals' else '_sat' + return _stellar_fraction(M_200, M1=pars['M1'+t], alpha=pars['alpha'+t], + gamma=pars['gamma'+t],delta=pars['delta'+t],epsilon=pars['epsilon'+t]) + +def _stellar_fraction(M_200, M1=1.526e11, alpha= -1.779, gamma = 0.547, delta = 4.394, epsilon=0.023): + """ + Internal function which compute the galaxy mass fractions. + + :param M_200: Halo mass inside a sphere which encompass a density 200 times larger than the critical density of the Universe. + :type array_like + :return: galaxy mass fractions. + :rtype: array_like + """ + fact = 10**( _g(np.log10(M_200/M1), alpha, gamma, delta) - _g(0, alpha, gamma, delta)) + return epsilon * (M1/M_200) * fact + +def _g(x, alpha= -1.779, gamma = 0.547, delta = 4.394): + return -np.log10( 10**(alpha*x) + 1) + \ + delta * ( (np.log10( 1+np.exp(x) ))**gamma ) / ( 1 + np.exp(10**(-x)) ) diff --git a/structure/baccoemu/baccoemu_vendored/lbias_expansion.py b/structure/baccoemu/baccoemu_vendored/lbias_expansion.py new file mode 100644 index 00000000..5af3f43b --- /dev/null +++ b/structure/baccoemu/baccoemu_vendored/lbias_expansion.py @@ -0,0 +1,748 @@ +import numpy as np +import copy +import pickle +import os +from .utils import _transform_space, MyProgressBar, mean_absolute_exp_percentage_error, accuracy_exp_01, accuracy_exp_005 +from scipy import interpolate +import tensorflow +tensorflow.compat.v1.logging.set_verbosity(tensorflow.compat.v1.logging.ERROR) +from tensorflow.keras.models import load_model +gpus = tensorflow.config.experimental.list_physical_devices('GPU') +if gpus: + for gpu in gpus: + tensorflow.config.experimental.set_memory_growth(gpu, True) +from .matter_powerspectrum import load_smeared_bao_emu + + +__all__ = ["Lbias_expansion"] + +class Lbias_expansion(object): + """ + A class to load and call the baccoemu for the Lagrangian bias expansion terms. + By default, the nonlinear Lagrangian bias expansion terms emulator (described + in Zennaro et al, 2021) and LPT lagrangian bias expansion terms emulator + (described in Aricò et al, 2021) are loaded. + Another emulator loaded by deafult is the IR-resummed linear power spectrum. + + :param lpt: whether to load the LPT emulator, defaults to True + :type lpt: boolean, optional + :param compute_sigma8: whether to load the sigma8 emulator, defaults to True + :type compute_sigma8: boolean, optional + :param smeared_bao: whether to load the smeared bao, defaults to True + :type smeared_bao: boolean, optional + :param nonlinear_boost: whether to load the nonlinear boost emulator, defaults to True + :type nonlinear_boost: boolean, optional + :param compute_sigma8: whether to load the sigma8 emulator, defaults to True + :type compute_sigma8: boolean, optional + :param verbose: whether to activate the verbose mode, defaults to True + :type verbose: boolean, optional + + """ + def __init__(self, lpt=True, smeared_bao=True, nonlinear_boost=True, compute_sigma8=True, verbose=True): + + self.verbose = verbose + + self.compute_lpt = True if lpt else False + self.compute_smeared_bao = True if smeared_bao else False + + self.cosmo_keys = np.array(['omega_cold', 'sigma8_cold', 'omega_baryon', 'ns', + 'hubble', 'neutrino_mass', 'w0', 'wa', 'expfactor']) + + self.lb_term_labels = [r'$1 1$', r'$1 \delta$', r'$1 \delta^2$', r'$1 s^2$', + r'$ 1 \nabla^2\delta$', r'$\delta \delta$', + r'$\delta \delta^2$', r'$\delta s^2$', + r'$\delta \nabla^2\delta$', r'$\delta^2 \delta^2$', + r'$\delta^2 s^2$', r'$\delta^2 \nabla^2\delta$', + r'$s^2 s^2$', r'$s^2 \nabla^2\delta$', + r'$\nabla^2\delta \nabla^2\delta$'] + + self.emulator = {} + if self.compute_lpt: + self.emulator['lpt'] = load_lpt_emu() + + if self.compute_smeared_bao: + self.emulator['smeared_bao'] = load_smeared_bao_emu() + + self.compute_nonlinear_boost = True if nonlinear_boost else False + if self.compute_nonlinear_boost: + self.emulator['nonlinear'] = load_nonlinear_lbias_emu() + + self.compute_sigma8 = True if compute_sigma8 else False + + if self.compute_sigma8: + from .matter_powerspectrum import Matter_powerspectrum + self.matter_powerspectrum_emulator = Matter_powerspectrum(linear=False, smeared_bao=False, + nonlinear_boost=False, baryonic_boost=False, + compute_sigma8=True, verbose=verbose) + + def _get_parameters(self, coordinates, which_emu, grid=None): + """ + Function that returns a dictionary of cosmological parameters, + computing derived cosmological parameters, if not + already present in the given coordinates, and checking the relevant boundaries. + :param coordinates: a set of coordinates in parameter space + :type coordinates: dict + :param which_emu: kind of emulator: options are 'linear', 'nonlinear','baryon','smeared_bao','sigma8' + :type which_emu: str + :param grid: dictionary with parameter and vector of values where to evaluate the emulator, defaults to None + :type grid: array_like, optional + :return: coordinates with derived parameters + :rtype: dict + """ + coordinates = {key: np.atleast_1d(coordinates[key]) for key in set(list(coordinates.keys())) - set(['k', 'k_lin', 'pk_lin'])} + + avail_pars = [coo for coo in coordinates.keys() if coordinates[coo][0] is not None] #parameters currently available + eva_pars = self.emulator[which_emu]['keys'] #parameters strictly needed to evaluate the emulator + req_pars = self.emulator[which_emu]['keys'] #parameters needed for a computation + comp_pars = list(set(req_pars)-set(avail_pars)) #parameters to be computed + deriv_pars = ['omega_cold','sigma8_cold', 'A_s'] #derived parameters that can be computed + miss_pars = list(set(comp_pars)-set(deriv_pars)) #parameters missing from coordinates + extra_pars = list(set(req_pars)-set(eva_pars)) #requested parameters not needed for evaluation + if miss_pars: + print(f"{which_emu} emulator:") + print(f" Please add the parameter(s) {miss_pars} to your coordinates!") + raise KeyError(f"{which_emu} emulator: coordinates need the following parameters: ", miss_pars) + + if ('omega_cold' in avail_pars) & ('omega_matter' in avail_pars): + assert len(coordinates['omega_cold']) == len(coordinates['omega_matter']), 'Both omega_cold and omega_matter were provided, but they have different len' + om_from_oc = coordinates['omega_cold'] + coordinates['neutrino_mass'] / 93.14 /coordinates['hubble']**2 + assert np.all(np.abs(coordinates['omega_matter'] - om_from_oc) < 1e-4), 'Both omega_cold and omega_matter were provided, but they are inconsistent among each other' + + if 'omega_cold' in comp_pars: + if 'omega_matter' not in avail_pars: + raise KeyError('One parameter between omega_matter and omega_cold must be provided!') + + omega_nu = coordinates['neutrino_mass'] / 93.14 /coordinates['hubble']**2 + coordinates['omega_cold'] = coordinates['omega_matter'] - omega_nu + + if ('sigma8_cold' not in avail_pars) & ('A_s' not in avail_pars): + raise KeyError('One parameter between sigma8_cold and A_s must be provided!') + + if ('sigma8_cold' in avail_pars) & ('A_s' in avail_pars): + #commented for the cases where one is computed and same value is repeated + #assert len(np.atleast_1d(coordinates['sigma8_cold'])) == len(atleast_1d(coordinates['A_s'])), 'Both sigma8_cold and A_s were provided, but they have different len' + + ignore_s8_pars = copy.deepcopy(coordinates) + del ignore_s8_pars['sigma8_cold'] + s8_from_A_s = self.matter_powerspectrum_emulator.get_sigma8(**ignore_s8_pars) + assert np.all(np.abs(coordinates['sigma8_cold'] - s8_from_A_s) < 1e-4), 'Both sigma8_cold and A_s were provided, but they are inconsistent among each other' + + if 'sigma8_cold' in comp_pars: + tmp_coords = copy.deepcopy(coordinates) + tmp_coords['cold']=True + coordinates['sigma8_cold'] = np.atleast_1d(self.matter_powerspectrum_emulator.get_sigma8(**tmp_coords)) + + if 'A_s' in comp_pars: + tmp_coords = copy.deepcopy(coordinates) + del tmp_coords['sigma8_cold'] + tmp_coords['A_s'] = 2e-9 + tmp_coords['cold'] = True + coordinates['A_s'] = np.atleast_1d((coordinates['sigma8_cold'] / self.matter_powerspectrum_emulator.get_sigma8(**tmp_coords))**2 * tmp_coords['A_s']) + + + pp = np.squeeze([coordinates[p][0] for p in eva_pars]) + coords_out = copy.deepcopy(coordinates) + + grid = {} + for key in coordinates.keys(): + if len(np.atleast_1d(coordinates[key])) > 1: + grid[key] = np.array(coordinates[key]) + + if len(list(grid.keys()))==0: + grid = None + else: + grid_structure = [] + for key in grid.keys(): + grid_structure.append(len(grid[key])) + grid_structure = np.array(grid_structure) + values, counts = np.unique(grid_structure, return_counts=True) + counts_but_highest = np.delete(counts, np.argmax(counts)) + assert np.all(counts == counts[0]) | np.all(counts_but_highest == 1), 'When passing multiple coordinate sets you should either vary only on parameter, or all parameters should have the same len' + + if grid is not None: + grid_pars = list(grid.keys()) # list of parameters that are varyied in a grid + N = len(grid[grid_pars[0]]) + pp = np.tile(pp, (N, 1)) + for par in grid_pars: + if par in eva_pars: + index = eva_pars.index(par) + pp[:,index] = np.float64(grid[par]) + if par in req_pars: + coords_out[par] = grid[par] + pp = np.float64(pp) + + for i,par in enumerate(eva_pars): + val = pp[i] if grid is None else pp[:,i] + message = 'Param {}={} out of bounds [{}, {}]'.format( + par, val, self.emulator[which_emu]['bounds'][i][0], self.emulator[which_emu]['bounds'][i][1]) + + assert np.all(val >= self.emulator[which_emu]['bounds'][i][0]) & np.all(val <= self.emulator[which_emu]['bounds'][i][1]), message + + if extra_pars: + cc = np.squeeze([coords_out[p] for p in extra_pars]) + if None in cc: + raise ValueError(f'None in parameters: {extra_pars} = {cc}!') + + return coords_out, pp, grid + + def get_galaxy_real_pk(self, bias=None, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, k=None, **kwargs): + """Compute the predicted galaxy auto pk and galaxy-matter cross pk given a set of bias parameters + + :param bias: a list of bias parameters, including b1, b2, bs2, blaplacian + :type bias: array-like + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + :param k: a vector of wavemodes in h/Mpc at which the nonlinear boost will be computed, if None + the default wavemodes of the nonlinear emulator will be used, defaults to None + :type k: array_like, optional + :return: k and P(k), a list of the emulated 15 LPT Lagrangian bias expansion terms + :rtype: tuple + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + import itertools + + assert len(bias) == 4, 'Please, pass a valid bias array, with b1, b2, bs2, blaplacian' + + k, pnn = self.get_nonlinear_pnn(**kwargs) + bias = np.concatenate(([1], bias)) + prod = np.array(list(itertools.combinations_with_replacement(np.arange(len(bias)), r=2))) + + pgal_auto = 0 + for i in range(len(pnn)): + fac = 2 if prod[i, 0] != prod[i,1] else 1 + pgal_auto += bias[prod[i, 0]] * bias[prod[i,1]] * fac * pnn[i] + pgal_cross = np.dot(bias, pnn[:5]) + + return k, pgal_auto, pgal_cross + + + def get_nonlinear_pnn(self, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, k=None, **kwargs): + """Compute the prediction of the nonlinear cold matter power spectrum. + + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + :param k: a vector of wavemodes in h/Mpc at which the nonlinear boost will be computed, if None + the default wavemodes of the nonlinear emulator will be used, defaults to None + :type k: array_like, optional + :return: k and P(k), a list of the emulated 15 LPT Lagrangian bias expansion terms + :rtype: tuple + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + if not self.compute_nonlinear_boost: + raise ValueError("Please enable the l-bias nonlinear boost!") + + coordinates, pp, grid = self._get_parameters(kwargs, 'nonlinear') + emulator = self.emulator['nonlinear'] + + n_log = [0,1,5,9,12,14] + n_switch = [10 , 11] + + _pp = _transform_space(np.array([pp]), space_rotation=False, bounds=emulator['bounds']) + + tmp_coords = copy.deepcopy(kwargs) + tmp_coords['k'] = emulator['k'] + _, pk_lpt = self.get_lpt_pk(**tmp_coords) + _, pk_bao = self.get_smeared_bao_pk(**tmp_coords) + + pk_lpt[0] = pk_bao + pk_lpt[1] = pk_bao + pk_lpt[5] = pk_bao + pk_lpt[4] = -pk_bao * emulator['k']**2 + pk_lpt[8] = -pk_bao * emulator['k']**2 + pk_lpt[14] = pk_bao * emulator['k']**4 + + P_nn = [] + for n in range(15): + prediction = emulator['model'][n](_pp.reshape(-1,9), training=False) + # import pdb; pdb.set_trace() + if n in n_switch: + lpt_term = pk_lpt[n + 2] + else: + lpt_term = pk_lpt[n] + if n in n_log: + _this_P_nn = np.squeeze(np.exp(emulator['scaler'][n].inverse_transform(prediction)) * lpt_term) + else: + _this_P_nn = np.squeeze(emulator['scaler'][n].inverse_transform(prediction) * lpt_term) + if n in n_switch: + if len(pk_lpt.shape) == 2: + _this_P_nn[:25] = _this_P_nn[:25] * pk_lpt[n][:25] / lpt_term[:25] + else: + _this_P_nn[:,:25] = _this_P_nn[:,:25] * pk_lpt[n,:,:25] / lpt_term[:,:25] + P_nn.append(_this_P_nn) + + if k is not None: + if max(k) > 0.75: + raise ValueError(f""" + The maximum k of the l-bias nonlinear emulator must be 0.75 h/Mpc: + the current value is {max(k)} h/Mpc""") + if (min(k) <= 1e-2)&(self.verbose): + print("WARNING: the nonlinear emulator is extrapolating to k < 0.01 h/Mpc!") + + new_P_nn = [] + for n in range(15): + unexpected_negative = np.any(P_nn[n] <= 0.0) # can happen when allowing extrapolation + if (n in n_log) & (unexpected_negative is False): + new_P_nn.append(np.exp(interpolate.interp1d(np.log(emulator['k']), np.log(P_nn[n]), kind='cubic', axis=0 if grid is None else 1, fill_value='extrapolate')(np.log(k)))) + else: + new_P_nn.append(interpolate.interp1d(np.log(emulator['k']), P_nn[n], kind='cubic', axis=0 if grid is None else 1, fill_value='extrapolate')(np.log(k))) + P_nn = np.array(new_P_nn) + else : + k = emulator['k'] + + return k, P_nn + + def get_lpt_pk(self, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, k=None, **kwargs): + """Compute the prediction of the 15 LPT Lagrangian bias expansion terms. + + + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + :param k: a vector of wavemodes in h/Mpc at which the nonlinear boost will be computed, if None + the default wavemodes of the nonlinear emulator will be used, defaults to None + :type k: array_like, optional + :return: k and P(k), a list of the emulated 15 LPT Lagrangian bias expansion terms + :rtype: tuple + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + if not self.compute_lpt: + raise ValueError("Please enable the lpt emulator!") + + emulator = self.emulator['lpt'] + coordinates, pp, grid = self._get_parameters(kwargs, 'lpt') + + sub = emulator['sub'] + scaler = emulator['scaler'] + + P_nn = [] + for n in range(15): + pred = emulator['model'][n](pp.reshape(-1,9), training=False) + prediction = np.squeeze(scaler[n].inverse_transform(pred)) + P_nn.append(prediction) + + if k is not None: + if max(k) > 0.75: + raise ValueError(f""" + The maximum k of the l-bias lpt emulator must be 0.75 h/Mpc: + the current value is {max(k)} h/Mpc""") + if (min(k) <= 1e-2)&(self.verbose): + print("WARNING: the l-bias lpt emulator is extrapolating to k < 0.01 h/Mpc!") + + for n in range(15): + p_interp = interpolate.interp1d(np.log(emulator['k']), P_nn[n], kind='cubic', axis=0 if grid is None else 1, fill_value='extrapolate', + assume_sorted=True) + P_nn[n] = p_interp(np.log(k)) + else : + k = emulator['k'] + + P_nn = np.array([np.exp(P_nn[n])-sub[n] for n in range(15)]) + return k, P_nn + + def get_smeared_bao_pk(self, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, k=None, **kwargs): + """Evaluate the smeared bao emulator at a set of coordinates in parameter space. + + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + :param k: a vector of wavemodes in h/Mpc at which the nonlinear boost will be computed, if None + the default wavemodes of the nonlinear emulator will be used, defaults to None + :type k: array_like, optional + :return: k and P(k), a list of the emulated 15 LPT Lagrangian bias expansion terms + :rtype: tuple + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + if not self.compute_smeared_bao: + raise ValueError("Please enable the smeared bao emulator!") + + emulator = self.emulator['smeared_bao'] + coordinates, pp, grid = self._get_parameters(kwargs, 'smeared_bao') + + ypred = emulator['model'](pp.reshape(-1,9), training=False) + pk_bao = np.squeeze(np.exp(emulator['scaler'].inverse_transform(ypred))) + if k is not None: + if (max(k) > 30.)|(min(k) < 1e-3): + raise ValueError(f""" + A minimum k > 0.001 h/Mpc and a maximum k < 30 h/Mpc + are required for the linear emulator: + the current values are {min(k)} h/Mpc and {max(k)} h/Mpc + """) + + else: + pk_bao = np.exp(interpolate.interp1d(np.log(emulator['k']),np.log(pk_bao), kind='cubic', axis=0 if grid is None else 1, fill_value='extrapolate')(np.log(k))) + else: + k = emulator['k'] + return k, pk_bao + + def _vector_eval_lpt_pk(self, coordinates): + """Get the lpt pk in a fast vectorized way + + Note: no checks are performed to increase speed, be aware of what you do + + :param coordinates: 2D array with a grid of cooordinates + :type coordinates: array + """ + emulator = self.lpt_emulator + sub = emulator['sub'] + scaler = emulator['scaler'] + + P_nn = [] + for n in range(15): + pred = emulator['model'][n](coordinates, training=False) + prediction = np.exp(scaler[n].inverse_transform(pred)) - sub[n] + P_nn.append(prediction) + return np.array(P_nn) + + def _vector_eval_smeared_bao_pk(self, coordinates): + """Get the nonlinear pk in a fast vectorized way + + Note: no checks are performed to increase speed, be aware of what you do + + :param coordinates: 2D array with a grid of cooordinates + :type coordinates: array + """ + emulator = self.smeared_bao_emulator + ypred = emulator['model'](coordinates, training=False) + pk_bao = np.exp(emulator['scaler'].inverse_transform(ypred)) + return pk_bao + + + def _vector_eval_nonlinear_pnn(self, coordinates): + """Get the nonlinear pk in a fast vectorized way + + Note: no checks are performed to increase speed, be aware of what you do + + :param coordinates: 2D array with a grid of cooordinates + :type coordinates: array + """ + from scipy.interpolate import interp1d + emulator = self.nonlinear_emulator + _pp = _transform_space(coordinates, space_rotation=False, bounds=emulator['bounds']) + + pk_lpt = self._vector_eval_lpt_pk(coordinates) + pk_bao = self._vector_eval_smeared_bao_pk(coordinates) + pk_bao_interp = np.exp(interp1d(np.log(self.smeared_bao_emulator['k']), np.log(pk_bao), kind='cubic')(np.log(emulator['k']))) + + pk_lpt[0,:,:] = pk_bao_interp + pk_lpt[1,:,:] = pk_bao_interp + pk_lpt[5,:,:] = pk_bao_interp + pk_lpt[4,:,:] = -pk_bao_interp * emulator['k']**2 + pk_lpt[8,:,:] = -pk_bao_interp * emulator['k']**2 + pk_lpt[14,:,:] = pk_bao_interp * emulator['k']**4 + + P_nn = [] + n_log = [0,1,5,9,12,14] + n_switch = [10 , 11] + for n in range(15): + prediction = emulator['model'][n](_pp, training=False) + if n in n_switch: + lpt_term = pk_lpt[n + 2,:,:] + else: + lpt_term = pk_lpt[n,:,:] + if n in n_log: + _this_P_nn = np.exp(emulator['scaler'][n].inverse_transform(prediction)) * lpt_term + else: + _this_P_nn = emulator['scaler'][n].inverse_transform(prediction) * lpt_term + if n in n_switch: + _this_P_nn[:,:25] = _this_P_nn[:,:25] * pk_lpt[n,:,:25] / lpt_term[:,:25] + P_nn.append(_this_P_nn) + + return np.array(P_nn) + + def _vector_eval_galaxy_real_pk(self, coordinates, bias): + """Get the galaxy pk in a fast vectorized way + + Note: no checks are performed to increase speed, be aware of what you do + + :param coordinates: 2D array with a grid of cooordinates + :type coordinates: array + :param bias: 2D array with a grid of bias + :type bias: array + """ + import itertools + + pnn = self._vector_eval_nonlinear_pnn(coordinates) + bias = np.hstack([np.full((len(coordinates), 1), 1), bias]) + prod = np.array(list(itertools.combinations_with_replacement(np.arange(bias.shape[1]), r=2))) + + pgal_auto = np.zeros((pnn.shape[1], pnn.shape[2])) + pgal_cross = np.zeros((pnn.shape[1], pnn.shape[2])) + for i in range(len(pnn)): + fac = 2 if prod[i, 0] != prod[i, 1] else 1 + pgal_auto += (bias[:, prod[i, 0]] * bias[:, prod[i, 1]] * fac * pnn[i, :, :].T).T + if i < 5: + pgal_cross += (bias[:, i] * pnn[i, :, :].T).T + + return self.nonlinear_emulator['k'], pgal_auto, pgal_cross + + def Hubble(self, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, k=None, **kwargs): + """Hubble function in km / s / Mpc + + Warning: neutrino and dynamical dark energy not yet implemented + Warning: no radiation included + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + O_m = kwargs['omega_matter'] if kwargs['omega_matter'] is not None else kwargs['omega_cold'] + O_Lambda = 1 - O_m + return kwargs['hubble'] * 100 * np.sqrt(O_m / expfactor**3 + O_Lambda) + + def comoving_radial_distance(self, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, k=None, **kwargs): + """Comoving radial distance in Mpc/h + + Warning: neutrino and dynamical dark energy not yet implemented + Warning: no radiation included + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + hpars = copy.deepcopy(kwargs) + del hpars['expfactor'] + + nz = int(1e4) + z_in = 0 + z_fin = 1 / kwargs['expfactor'] - 1 + z_vector = np.linspace(z_in, z_fin, nz) + a_vector = 1 / (1 + z_vector) + H_vector = self.Hubble(expfactor=a_vector, **hpars) + return 3e3 * hpars['hubble']*100 / np.trapz(H_vector, x=z_vector) + + +def load_lpt_emu(verbose=True): + """Loads in memory the lpt emulator described in Aricò et al. 2021. + + :return: a dictionary containing the emulator object + :rtype: dict + """ + + if verbose: + print('Loading l-bias lpt emulator...') + + basefold = os.path.dirname(os.path.abspath(__file__)) + + old_names = [(basefold + '/' + "lpt_emulator")] + for old_name in old_names: + if os.path.exists(old_name): + import shutil + shutil.rmtree(old_name) + + emulator_name = (basefold + '/' + + "lpt_emulator_v2.0.0") + + if (not os.path.exists(emulator_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading emulator data (34 Mb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/lpt_emulator_v2.0.0.tar', + emulator_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_name + '.tar') + + customs={"accuracy_01": accuracy_exp_01, + "accuracy_005": accuracy_exp_005, + "mean_absolute_exp_percentage_error":mean_absolute_exp_percentage_error} + + emulator = {} + emulator['emu_type'] = 'nn' + + emulator['model'] = [] + emulator['sub'] = [] + emulator['scaler'] = [] + for n in range(15): + i_emulator_name = f'{emulator_name}/lpt_emu_field{n}' + + file_to_read = open(f"{i_emulator_name}/details.pickle", "rb") + nn_details = pickle.load(file_to_read) + + emulator['model'].append(load_model(i_emulator_name, custom_objects=customs)) + emulator['scaler'].append(nn_details['scaler']) + emulator['sub'].append(nn_details['subtract']) + + emulator['k'] = nn_details['kk'] + emulator['keys'] = ['omega_cold', 'sigma8_cold', 'omega_baryon', 'ns', + 'hubble', 'neutrino_mass', 'w0', 'wa', 'expfactor'] + emulator['bounds'] = nn_details['bounds'] + + if verbose: + print('L-bias lpt emulator loaded in memory.') + + return emulator + +def load_nonlinear_lbias_emu(emu_type='nn', verbose=True): + """Loads in memory the nonlinear emulator described in Zennaro et al. 2021. + + :param emu_type: type of emulator, can be 'gp' for the gaussian process, ot + 'nn' for the neural network + :type emu_type: str + + :return: a dictionary containing the emulator object + :rtype: dict + """ + if verbose: + print('Loading non-linear l-bias emulator...') + + basefold = os.path.dirname(os.path.abspath(__file__)) + + emulator_name = (basefold + '/' + + "lbias_emulator") + + if (not os.path.exists(emulator_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading emulator data (34Mb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/lbias_emulator.tar', + emulator_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_name + '.tar') + emulator = {} + emulator['emu_type'] = 'nn' + emulator['model'] = [] + for n in range(15): + i_emulator_name = f'{emulator_name}/lbias_emu_field{n}' + emulator['model'].append(load_model(i_emulator_name)) + + with open(emulator_name + '/lbias_emu.pickle', 'rb') as f: + emulator['scaler'] = pickle.load(f) + emulator['npca'] = pickle.load(f) + emulator['k'] = pickle.load(f) + _ = pickle.load(f) # components + emulator['rotation'] = pickle.load(f) + emulator['bounds'] = pickle.load(f) + emulator['keys'] = ['omega_cold', 'sigma8_cold', 'omega_baryon', 'ns', + 'hubble', 'neutrino_mass', 'w0', 'wa', 'expfactor'] + + if verbose: + print('Nonlinear l-bias emulator loaded in memory.') + return emulator diff --git a/structure/baccoemu/baccoemu_vendored/matter_powerspectrum.py b/structure/baccoemu/baccoemu_vendored/matter_powerspectrum.py new file mode 100644 index 00000000..ab1d84b6 --- /dev/null +++ b/structure/baccoemu/baccoemu_vendored/matter_powerspectrum.py @@ -0,0 +1,1595 @@ +import numpy as np +import copy +import pickle +import os + +from numpy.core.shape_base import atleast_1d +from scipy.sparse import coo +from .utils import _transform_space, MyProgressBar, mean_absolute_exp_percentage_error, accuracy_exp_002, accuracy_exp_005 + +import tensorflow +tensorflow.compat.v1.logging.set_verbosity(tensorflow.compat.v1.logging.ERROR) +from tensorflow.keras.models import load_model +from .baryonic_boost import load_baryonic_emu +from scipy import interpolate +gpus = tensorflow.config.experimental.list_physical_devices('GPU') +if gpus: + for gpu in gpus: + tensorflow.config.experimental.set_memory_growth(gpu, True) + +__all__ = ["Matter_powerspectrum"] + +class Matter_powerspectrum(object): + """ + A class to load and call the baccoemu for the matter powerspectrum. + By default, the linear power spectrum (described in Aricò et al. 2021), the nonlinear boost + (described in Angulo et al. 2020), and the baryonic boost (described in Aricò et al. 2020c) are loaded. + + The emulators for A_s, sigma8, sigma12, and the smearing and removing of the BAO are also available. + + At large scales, i.e. k<0.01, the ratios between non-linear/linear and baryonic/non-linear spectra is considered to be 1. + + :param linear: whether to load the linear emulator, defaults to True + :type linear: boolean, optional + :param smeared_bao: whether to load the smeared-BAO emulator, defaults to False + :type smeared_bao: boolean, optional + :param no_wiggles: whether to load the no-wiggles emulator, defaults to True + :type no_wiggles: boolean, optional + :param nonlinear_boost: whether to load the nonlinear boost emulator, defaults to True + :type nonlinear_boost: boolean, optional + :param baryonic_boost: whether to load the baryonic boost emulator, defaults to True + :type baryonic_boost: boolean, optional + :param compute_sigma8: whether to load the sigma8 emulator, defaults to True + :type compute_sigma8: boolean, optional + :param verbose: whether to activate the verbose mode, defaults to True + :type verbose: boolean, optional + + """ + def __init__(self, linear=True, smeared_bao=False, no_wiggles=True, nonlinear_boost = True, + baryonic_boost=True, compute_sigma8=True, + nonlinear_emu_path=None, nonlinear_emu_details=None, + baryonic_emu_path=None, baryonic_emu_details=None, verbose=True): + + self.verbose = verbose + + self.compute_linear = True if linear else False + self.compute_smeared_bao = True if smeared_bao else False + self.compute_no_wiggles = True if no_wiggles else False + self.compute_nonlinear_boost = True if nonlinear_boost else False + self.compute_baryonic_boost = True if baryonic_boost else False + self.compute_sigma8 = True if compute_sigma8 else False + + if self.compute_smeared_bao and self.compute_no_wiggles: + print("""Provide only one between the smeared BAO and the non-wiggles emulators!""") + raise ValueError("""Set either compute_smeared_bao=False (default) or no_wiggles=False!""") + + self.emulator = {} + if self.compute_sigma8: + self.emulator['sigma8'] = load_sigma8_emu(verbose=verbose) + + if self.compute_linear: + self.emulator['linear'] = load_linear_emu(verbose=verbose) + + if self.compute_smeared_bao: + self.emulator['smeared_bao'] = load_smeared_bao_emu(verbose=verbose) + + if self.compute_no_wiggles: + self.emulator['no_wiggles'] = load_no_wiggles_emu(verbose=verbose) + + if self.compute_nonlinear_boost: + self.emulator['nonlinear'] = load_nonlinear_emu(fold_name=nonlinear_emu_path, detail_name=nonlinear_emu_details, verbose=verbose) + + if self.compute_baryonic_boost: + self.emulator['baryon'] = load_baryonic_emu(fold_name=baryonic_emu_path, detail_name=baryonic_emu_details, verbose=verbose) + + def _get_parameters(self, coordinates, which_emu): + """ + Function that returns a dictionary of cosmological parameters, + computing derived cosmological parameters, if not + already present in the given coordinates, and checking the relevant boundaries. + :param coordinates: a set of coordinates in parameter space + :type coordinates: dict + :param which_emu: kind of emulator: options are 'linear', 'nonlinear','baryon','smeared_bao','sigma8' + :type which_emu: str + :return: coordinates with derived parameters + :rtype: dict + """ + coordinates = {key: np.atleast_1d(coordinates[key]) for key in set(list(coordinates.keys())) - set(['k', 'k_lin', 'pk_lin'])} + + avail_pars = [coo for coo in coordinates.keys() if coordinates[coo][0] is not None] #parameters currently available + eva_pars = self.emulator[which_emu]['keys'] #parameters strictly needed to evaluate the emulator + req_pars = self.emulator[which_emu]['keys'] if which_emu != 'linear' else self.emulator[which_emu]['full_keys'] #parameters needed for a computation + comp_pars = list(set(req_pars)-set(avail_pars)) #parameters to be computed + deriv_pars = ['omega_cold','sigma8_cold', 'A_s'] #derived parameters that can be computed + miss_pars = list(set(comp_pars)-set(deriv_pars)) #parameters missing from coordinates + extra_pars = list(set(req_pars)-set(eva_pars)) #requested parameters not needed for evaluation + if miss_pars: + print(f"{which_emu} emulator:") + print(f" Please add the parameter(s) {miss_pars} to your coordinates!") + raise KeyError(f"{which_emu} emulator: coordinates need the following parameters: ", miss_pars) + + if ('omega_cold' in avail_pars) & ('omega_matter' in avail_pars): + assert len(coordinates['omega_cold']) == len(coordinates['omega_matter']), 'Both omega_cold and omega_matter were provided, but they have different len' + om_from_oc = coordinates['omega_cold'] + coordinates['neutrino_mass'] / 93.14 /coordinates['hubble']**2 + assert np.all(np.abs(coordinates['omega_matter'] - om_from_oc) < 1e-4), 'Both omega_cold and omega_matter were provided, but they are inconsistent among each other' + + if 'omega_cold' in comp_pars: + if 'omega_matter' not in avail_pars: + raise KeyError('One parameter between omega_matter and omega_cold must be provided!') + + omega_nu = coordinates['neutrino_mass'] / 93.14 /coordinates['hubble']**2 + coordinates['omega_cold'] = coordinates['omega_matter'] - omega_nu + + if ('sigma8_cold' not in avail_pars) & ('A_s' not in avail_pars): + raise KeyError('One parameter between sigma8_cold and A_s must be provided!') + + if ('sigma8_cold' in avail_pars) & ('A_s' in avail_pars): + #commented for the cases where one is computed and same value is repeated + #assert len(np.atleast_1d(coordinates['sigma8_cold'])) == len(atleast_1d(coordinates['A_s'])), 'Both sigma8_cold and A_s were provided, but they have different len' + + ignore_s8_pars = copy.deepcopy(coordinates) + del ignore_s8_pars['sigma8_cold'] + ignore_s8_pars['cold'] = True + s8_from_A_s = self.get_sigma8(**ignore_s8_pars) + assert np.all(np.abs(coordinates['sigma8_cold'] - s8_from_A_s) < 1e-4), 'Both sigma8_cold and A_s were provided, but they are inconsistent among each other' + + if 'sigma8_cold' in comp_pars: + tmp_coords = copy.deepcopy(coordinates) + tmp_coords['cold']=True + coordinates['sigma8_cold'] = np.atleast_1d(self.get_sigma8(**tmp_coords)) + + if 'A_s' in comp_pars: + tmp_coords = copy.deepcopy(coordinates) + del tmp_coords['sigma8_cold'] + tmp_coords['A_s'] = 2e-9 + tmp_coords['cold'] = True + coordinates['A_s'] = np.atleast_1d((coordinates['sigma8_cold'] / self.get_sigma8(**tmp_coords))**2 * tmp_coords['A_s']) + + pp = np.squeeze([coordinates[p][0] for p in eva_pars]) + coords_out = copy.deepcopy(coordinates) + + grid = {} + for key in coordinates.keys(): + if len(np.atleast_1d(coordinates[key])) > 1: + grid[key] = np.array(coordinates[key]) + + if len(list(grid.keys()))==0: + grid = None + else: + grid_structure = [] + for key in grid.keys(): + grid_structure.append(len(grid[key])) + grid_structure = np.array(grid_structure) + values, counts = np.unique(grid_structure, return_counts=True) + counts_but_highest = np.delete(counts, np.argmax(counts)) + assert np.all(counts == counts[0]) | np.all(counts_but_highest == 1), 'When passing multiple coordinate sets you should either vary only on parameter, or all parameters should have the same len' + + if grid is not None: + grid_pars = list(grid.keys()) # list of parameters that are varyied in a grid + N = len(grid[grid_pars[0]]) + pp = np.tile(pp, (N, 1)) + for par in grid_pars: + if par in eva_pars: + index = eva_pars.index(par) + pp[:,index] = np.float64(grid[par]) + if par in req_pars: + coords_out[par] = grid[par] + pp = np.float64(pp) + + for i,par in enumerate(eva_pars): + val = pp[i] if grid is None else pp[:,i] + message = 'Param {}={} out of bounds [{}, {}]'.format( + par, val, self.emulator[which_emu]['bounds'][i][0], self.emulator[which_emu]['bounds'][i][1]) + + assert np.all(val >= self.emulator[which_emu]['bounds'][i][0]) & np.all(val <= self.emulator[which_emu]['bounds'][i][1]), message + + if extra_pars: + cc = np.squeeze([coords_out[p] for p in extra_pars]) + if None in cc: + raise ValueError(f'None in parameters: {extra_pars} = {cc}!') + + return coords_out, pp, grid + + + def _evaluate_nonlinear(self, **kwargs): + """Evaluate the given emulator at a set of coordinates in parameter space. + + The coordinates must be specified as a dictionary with the following + keywords: + + #. 'omega_cold' + #. 'omega_baryon' + #. 'sigma8_cold' + #. 'hubble' + #. 'ns' + #. 'neutrino_mass' + #. 'w0' + #. 'wa' + #. 'expfactor' + #. 'k' : a vector of wavemodes in h/Mpc at which the nonlinear boost will be computed, if None + the default wavemodes of the nonlinear emulator will be used, defaults to None + #. 'k_lin': a vector of wavemodes in h/Mpc, if None the wavemodes used by + the linear emulator are returned, defaults to None + #. 'pk_lin': a vector of linear matter power spectrum computed at k_lin, either cold or total depending + on the key "cold". If None the linear emulator will be called, defaults to None + #. 'cold': whether to return the cold matter power spectrum or the total one. Default to True + """ + if not self.compute_nonlinear_boost: + raise ValueError("Please enable the nonlinear boost!") + + k = kwargs['k'] if 'k' in kwargs.keys() else None + k_lin = kwargs['k_lin'] if 'k_lin' in kwargs.keys() else None + pk_lin = kwargs['pk_lin'] if 'pk_lin' in kwargs.keys() else None + cold = kwargs['cold'] if 'cold' in kwargs.keys() else True + + emulator = self.emulator['nonlinear'] + + coords, pp, grid = self._get_parameters(kwargs, 'nonlinear') + + _pp = _transform_space(pp, space_rotation=False, bounds=emulator['bounds']) + + yrec = emulator['model'](_pp.reshape(-1,9), training=False) + Q = np.squeeze(np.exp(emulator['scaler'].inverse_transform(yrec))) + + if pk_lin is None: + pk_lin_kwargs = copy.deepcopy(kwargs) + pk_lin_kwargs['k'] = None + k_lin, pk_lin = self.get_linear_pk(**pk_lin_kwargs) + else: + k_lin = np.squeeze(k_lin) + pk_lin = np.squeeze(pk_lin) + if k_lin is None: + raise ValueError("""If the linear power spectrum pk_lin is provided, + also the wavenumbers k_lin at which is computed must be provided """) + elif (np.amin(k_lin)>1e-3) | (np.amax(k_lin) < 10): + raise ValueError(f""" + A minimum k <= 0.001 h/Mpc and a maximum k >= 10 h/Mpc + are required in the linear power spectrum for the calculation + of the non linear boost: + the current values are {np.amin(k_lin)}) h/Mpc and {np.amax(k_lin)} h/Mpc + """) + if cold: + k_lin_cold = k_lin + pk_lin_cold = pk_lin + else: + k_lin_tot = k_lin + pk_lin_tot = pk_lin + total_kwargs = copy.deepcopy(kwargs) + total_kwargs['cold'] = True + k_lin_cold, pk_lin_cold = self.get_linear_pk(**total_kwargs) + + pklin_interp_cold = interpolate.interp1d(np.log(k_lin_cold), np.log(pk_lin_cold), + kind='linear', axis = -1 if grid is None else 1, fill_value='extrapolate') + pk_lin_emu_cold = np.exp(pklin_interp_cold(np.log(emulator['k']))) + + if not cold: + pklin_interp_tot = interpolate.interp1d(np.log(k_lin_tot), np.log(pk_lin_tot), + kind='linear', axis = -1 if grid is None else 1, fill_value='extrapolate') + pk_lin_emu_tot = np.exp(pklin_interp_tot(np.log(emulator['k']))) + pk_lin_emu = pk_lin_emu_tot + else: + pk_lin_emu = pk_lin_emu_cold + + if self.compute_smeared_bao: + smeared_bao_kwargs = copy.deepcopy(kwargs) + smeared_bao_kwargs['k'] = emulator['k'] + _, pk_smeared = self.get_smeared_bao_pk(**smeared_bao_kwargs) + + elif self.compute_no_wiggles: + no_wiggles_kwargs = copy.deepcopy(kwargs) + no_wiggles_kwargs['k'] = emulator['k'] + no_wiggles_kwargs['k_lin'] = k_lin_cold + no_wiggles_kwargs['pk_lin'] = pk_lin_cold + _, pk_no_wiggles = self.get_no_wiggles_pk(**no_wiggles_kwargs) + pk_smeared = _smeared_bao_pk(k_lin=k_lin_cold, pk_lin=pk_lin_cold, k_emu=emulator['k'], pk_lin_emu=pk_lin_emu_cold, pk_nw=pk_no_wiggles, grid=grid) + else: + pk_smeared = _smeared_bao_pk(k_lin=k_lin_cold, pk_lin=pk_lin_cold, k_emu=emulator['k'], pk_lin_emu=pk_lin_emu_cold, grid=grid) + + if cold: + nonlinear_boost = Q * pk_smeared / pk_lin_emu_cold + else: + omega_nu = kwargs['neutrino_mass'] / 93.14 / kwargs['hubble']**2 + omega_matter = kwargs['omega_cold'] + omega_nu if kwargs['omega_cold'] is not None else kwargs['omega_matter'] + f_nu = omega_nu/omega_matter + if grid is None: + add = pk_lin_emu_tot - (1-f_nu)**2 * pk_lin_emu_cold + else: + add = pk_lin_emu_tot - ((1-f_nu)**2 * pk_lin_emu_cold.T).T + + if grid is None: + nonlinear_boost = (Q*pk_smeared*(1-f_nu)**2 + add)/pk_lin_emu_tot + else: + nonlinear_boost = (((Q*pk_smeared).T*(1-f_nu)**2).T + add)/pk_lin_emu_tot + + if k is not None: + k_max = max(k) + kemu_max = self.emulator['nonlinear']['k'].max() + kemu_min = self.emulator['nonlinear']['k'].min() + + if (k_max > kemu_max): + raise ValueError(f"""The nonlinear emulator must be {kemu_max} h/Mpc: the current value is {k_max} h/Mpc""") + + nl_interp = interpolate.interp1d(np.log(emulator['k']), nonlinear_boost, + kind='linear', axis = -1 if grid is None else 1, fill_value='extrapolate', bounds_error=False, + assume_sorted=True) + nonlinear_boost = np.squeeze(nl_interp(np.log(k))) + if grid is None: + nonlinear_boost[k kemu_max) & (self.verbose): + print(f""" WARNING: + The maximum k of the baryonic boost emulator is {kemu_max} h/Mpc: + the baryonic emulator is currently extrapolating to {k_max} h/Mpc; + the extrapolation will likely be not accurate. + """) + + baryonic_interp = interpolate.interp1d(np.log(emulator['k']), baryonic_boost, + kind='linear', axis = 0 if grid is None else 1, fill_value='extrapolate', bounds_error=False, + assume_sorted=True) + baryonic_boost = baryonic_interp(np.log(k)) + if grid is None: + baryonic_boost[k max(emulator['k']))|(min(k) < min(emulator['k'])): + raise ValueError(f""" + A minimum k > {min(emulator['k'])} h/Mpc and a maximum k < {max(emulator['k'])} h/Mpc + are required for the linear emulator: + the current values are {min(k)} h/Mpc and {max(k)} h/Mpc + """) + + else: + if np.any(pk_lin < 0): + pklin_interp = interpolate.interp1d(np.log(emulator['k']), pk_lin, + kind='linear', axis = 0 if grid is None else 1, fill_value='extrapolate', + assume_sorted=True) + pk_lin = pklin_interp(np.log(k)) + else: + pklin_interp = interpolate.interp1d(np.log(emulator['k']), np.log(pk_lin), + kind='linear', axis = 0 if grid is None else 1, fill_value='extrapolate', + assume_sorted=True) + pk_lin = np.exp(pklin_interp(np.log(k))) + else: + k = emulator['k'] + return k, pk_lin + + + def get_smeared_bao_pk(self, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, k=None, **kwargs): + """Evaluate the cold matter power spectrum with smeared bao, calling an emulator + at a set of coordinates in parameter space. + + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + :param k: a vector of wavemodes in h/Mpc at which the nonlinear boost will be computed, if None + the default wavemodes of the linear emulator will be used, defaults to None + :type k: array_like, optional + :param k_lin: a vector of wavemodes in h/Mpc, if None the wavemodes used by + the linear emulator are returned, defaults to None + :type k_lin: array_like, optional + :param pk_lin: a vector of linear power spectrum computed at k_lin, if None + the linear emulator will be called, defaults to None + :type pk_lin: array_like, optional + :param grid: dictionary with parameters and vectors of values where to evaluate the emulator, defaults to None + :type grid: array_like, optional + + :return: k and the linear P(k) + :rtype: tuple + """ + + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + if not self.compute_smeared_bao: + raise ValueError("Please enable the smeared bao emulator!") + + emulator = self.emulator['smeared_bao'] + coordinates, pp, grid = self._get_parameters(kwargs, 'smeared_bao') + + ypred = emulator['model'](pp.reshape(-1,9), training=False) + pk_bao = np.squeeze(np.exp(emulator['scaler'].inverse_transform(ypred))) + + if k is not None: + if (max(k) > max(emulator['k']))|(min(k) < min(emulator['k'])): + raise ValueError(f""" + A minimum k > 0.001 h/Mpc and a maximum k < 30 h/Mpc + are required for the smeared-BAO emulator: + the current values are {min(k)} h/Mpc and {max(k)} h/Mpc + """) + else: + pk_bao_interp = interpolate.interp1d(np.log(emulator['k']), np.log(pk_bao), + kind='linear', axis = 0 if grid is None else 1, fill_value='extrapolate', + assume_sorted=True) + pk_bao = np.exp(pk_bao_interp(np.log(k))) + else: + k = emulator['k'] + return k, pk_bao + + def get_no_wiggles_pk(self, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, k=None, k_lin=None, pk_lin=None, **kwargs): + """Evaluate the cold matter power spectrum with no wiggles (no bao), calling an emulator + at a set of coordinates in parameter space. + + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + :param k: a vector of wavemodes in h/Mpc at which the nonlinear boost will be computed, if None + the default wavemodes of the linear emulator will be used, defaults to None + :type k: array_like, optional + :param k_lin: a vector of wavemodes in h/Mpc, if None the wavemodes used by + the linear emulator are returned, defaults to None + :type k_lin: array_like, optional + :param pk_lin: a vector of linear power spectrum computed at k_lin, if None + the linear emulator will be called, defaults to None + :type pk_lin: array_like, optional + :param grid: dictionary with parameters and vectors of values where to evaluate the emulator, defaults to None + :type grid: array_like, optional + + :return: k and the linear P(k) + :rtype: tuple + """ + + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + if not self.compute_no_wiggles: + raise ValueError("Please enable the no_wiggles emulator!") + + emulator = self.emulator['no_wiggles'] + coordinates, pp, grid = self._get_parameters(kwargs, 'no_wiggles') + + ypred = emulator['model'](pp.reshape(-1,7), training=False) + ypred = emulator['scaler'].inverse_transform(ypred) + ypred[:,emulator['k']>0.8] = 0. + ypred = np.squeeze(ypred) + + if k is not None: + if (min(k) < min(emulator['k'])) and (self.verbose): + print(f"""The no-wiggles emulator is extrapolating from {min(k)} to {min(emulator['k'])} h/Mpc!""") + + pk_no_wiggles_interp = interpolate.interp1d(np.log(emulator['k']), ypred, + kind='linear', axis = 0 if grid is None else 1, fill_value=0., bounds_error=False, + assume_sorted=True) + pk_no_wiggles_plin_ratio = np.exp(pk_no_wiggles_interp(np.log(k))) + else: + k = emulator['k'] + pk_no_wiggles_plin_ratio = np.squeeze(np.exp(ypred)) + + if (k_lin is not None) & (pk_lin is not None): + do_interp = False + if len(k_lin) != len(k): + do_interp = True + else: + if np.any(k_lin != k): + do_interp = True + if do_interp: + pk_lin_interp = interpolate.interp1d(np.log(k_lin), np.log(pk_lin), + kind='linear', axis = 0 if grid is None else 1, fill_value='extrapolate', + assume_sorted=True) + pk_lin = np.exp(pk_lin_interp(np.log(k))) + + else: + linear_kwargs = copy.deepcopy(kwargs) + linear_kwargs['k'] = k + k_lin, pk_lin = self.get_linear_pk(**linear_kwargs) + + pk_no_wiggles = pk_no_wiggles_plin_ratio * pk_lin + + return k, pk_no_wiggles + + def get_sigma8(self, cold=True, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, **kwargs): + """Return sigma8, calling an emulator + at a set of coordinates in parameter space. + + :param cold: whether to return sigma8_cold (cdm+baryons) or total matter + :type cold: bool + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + + :return: sigma8 + :rtype: float or array + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + if self.compute_sigma8: + emulator = self.emulator['sigma8'] + else: + raise ValueError("Please enable the sigma8 emulator!") + + coordinates, pp, grid = self._get_parameters(kwargs, 'sigma8') + ypred = np.squeeze(emulator['model'](pp.reshape(-1,7), training=False)) + if hasattr(coordinates['A_s'], '__len__'): + ypred = (ypred.T * np.sqrt(coordinates['A_s']/1.e-9)).T + else: + ypred = ypred * np.sqrt(coordinates['A_s']/1.e-9) + s8_index = 0 if cold else 1 + res = ypred[...,s8_index] + res = float(res) if res.ndim==0 else res + return res + + + def get_sigma12(self, cold=True, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, **kwargs): + """Return sigma12, calling an emulator + at a set of coordinates in parameter space. + + :param cold: whether to return sigma8_cold (cdm+baryons) or total matter + :type cold: bool + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + + :return: sigma12 + :rtype: float or array + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + if self.compute_sigma8: + emulator = self.emulator['sigma8'] + else: + raise ValueError("Please enable the sigma8 emulator!") + + coordinates, pp, grid = self._get_parameters(kwargs, 'sigma8') + ypred = np.squeeze(emulator['model'](pp.reshape(-1,7), training=False)) + ypred = ypred * np.sqrt(coordinates['A_s']/1.e-9) + s8_index = 2 if cold else 3 + res = ypred[...,s8_index] + res = float(res) if res.ndim==0 else res + return res + + + def get_A_s(self, omega_cold=None, omega_matter=None, omega_baryon=None, + sigma8_cold=None, A_s=None, hubble=None, ns=None, neutrino_mass=None, + w0=None, wa=None, expfactor=None, **kwargs): + """Return A_s, corresponding to a given sigma8_cold + at a set of coordinates in parameter space. + + :param omega_cold: omega cold matter (cdm + baryons), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_cold: float or array + :param omega_matter: omega total matter (cdm + baryons + neutrinos), either omega_cold + or omega_matter should be specified, if both are specified + they should be consistent + :type omega_matter: float or array + :param sigma8_cold: rms of cold (cdm + baryons) linear perturbations, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type sigma8_cold: float or array + :param A_s: primordial scalar amplitude at k=0.05 1/Mpc, either sigma8_cold + or A_s should be specified, if both are specified they should be + consistent + :type A_s: float or array + :param hubble: adimensional Hubble parameters, h=H0/(100 km/s/Mpc) + :type hubble: float or array + :param ns: scalar spectral index + :type ns: float or array + :param neutrino_mass: total neutrino mass + :type neutrino_mass: float or array + :param w0: dark energy equation of state redshift 0 parameter + :type w0: float or array + :param wa: dark energy equation of state redshift dependent parameter + :type wa: float or array + :param expfactor: expansion factor a = 1 / (1 + z) + :type expfactor: float or array + + :return: A_s + :rtype: float or array + """ + _kwargs = locals() + kwargs = {key: _kwargs[key] for key in set(list(_kwargs.keys())) - set(['self'])} + + coordinates, pp, grid = self._get_parameters(kwargs, 'linear') + + res = np.squeeze(coordinates['A_s']) + + return float(res) if res.ndim == 0 else res + + +def load_linear_emu(verbose=True): + """Loads in memory the linear emulator described in Aricò et al. 2021. + + :return: a dictionary containing the emulator object + :rtype: dict + """ + + if verbose: + print('Loading linear emulator...') + + basefold = os.path.dirname(os.path.abspath(__file__)) + emulator_cold_name = (basefold + '/' + + "cold_matter_linear_emu_1.0.0") + emulator_tot_name = (basefold + '/' + + "total_matter_linear_emu_1.0.0") + + old_names = [(basefold + '/' + "linear_emulator")] + for old_name in old_names: + if os.path.exists(old_name): + import shutil + shutil.rmtree(old_name) + + if (not os.path.exists(emulator_cold_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading emulator data (2 Mb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/cold_matter_linear_emu_1.0.0.tar', + emulator_cold_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_cold_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_cold_name + '.tar') + + if (not os.path.exists(emulator_tot_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading emulator data (2 Mb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/total_matter_linear_emu_1.0.0.tar', + emulator_tot_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_tot_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_tot_name + '.tar') + + customs = { + "accuracy_exp_002": accuracy_exp_002, + "accuracy_exp_005": accuracy_exp_005, + "mean_absolute_exp_percentage_error":mean_absolute_exp_percentage_error + } + metrics_list = ["accuracy",accuracy_exp_002, accuracy_exp_005] + + emulator = {} + emulator['emu_type'] = 'nn' + emulator['model_cold'] = load_model(emulator_cold_name, custom_objects=customs, compile=False) + emulator['model_cold'].compile(optimizer='adam', loss=mean_absolute_exp_percentage_error, metrics=metrics_list) + + file_to_read = open(f"{emulator_cold_name}/details.pickle", "rb") + nn_cold_details = pickle.load(file_to_read) + + emulator['k'] = nn_cold_details['kk'] + emulator['scaler_cold'] = nn_cold_details['scaler'] + + emulator['model_tot'] = load_model(emulator_tot_name, custom_objects=customs, compile=False) + emulator['model_tot'].compile(optimizer='adam', loss=mean_absolute_exp_percentage_error, metrics=metrics_list) + + file_to_read = open(f"{emulator_tot_name}/details.pickle", "rb") + nn_tot_details = pickle.load(file_to_read) + + emulator['scaler_tot'] = nn_tot_details['scaler'] + emulator['keys'] = ['omega_cold', 'omega_baryon', + 'hubble', 'neutrino_mass', 'w0', 'wa', 'expfactor'] + emulator['full_keys'] = ['omega_cold', 'omega_baryon', 'A_s', 'ns', + 'hubble', 'neutrino_mass', 'w0', 'wa', 'expfactor'] + + emulator['bounds'] = nn_tot_details['bounds']#{key: nn_cold_details['bounds'][i] for i, key in enumerate(emulator['keys'])} + + if verbose: + print('Linear emulator loaded in memory.') + + return emulator + +def load_sigma8_emu(verbose=True): + """Loads in memory an emulator to quickly pass from A_s to sigma8. + + :return: a dictionary containing the emulator object + :rtype: dict + """ + + if verbose: + print('Loading sigma8 emulator...') + + basefold = os.path.dirname(os.path.abspath(__file__)) + + old_names = [(basefold + '/' + "sigma8_emu_1.0.0"),] + for old_name in old_names: + if os.path.exists(old_name): + import shutil + shutil.rmtree(old_name) + + emulator_name = (basefold + '/' + + "sigma8_emu_2.0.0") + + if (not os.path.exists(emulator_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading emulator data (141 Kb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/sigma8_emu_2.0.0.tar', + emulator_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_name + '.tar') + + + customs = { + "accuracy_exp_002": accuracy_exp_002, + "accuracy_exp_005": accuracy_exp_005, + "mean_absolute_exp_percentage_error":mean_absolute_exp_percentage_error + } + metrics_list = ["accuracy",accuracy_exp_002, accuracy_exp_005] + + emulator = {} + emulator['emu_type'] = 'nn' + emulator['model'] = load_model(emulator_name, custom_objects=customs, compile=False) + emulator['model'].compile(optimizer='adam', loss=mean_absolute_exp_percentage_error, metrics=metrics_list) + + file_to_read = open(f"{emulator_name}/details.pickle", "rb") + nn_details = pickle.load(file_to_read) + emulator['bounds'] = nn_details['bounds'] + emulator['keys'] = ['omega_cold', 'omega_baryon', 'ns', 'hubble', 'neutrino_mass', 'w0', 'wa'] + + if verbose: + print('Sigma8 emulator loaded in memory.') + + return emulator + +def load_smeared_bao_emu(verbose=True): + """Loads in memory the smeared bao emulator. + + :return: a dictionary containing the emulator object + :rtype: dict + """ + + if verbose: + print('Loading smeared bao emulator...') + + basefold = os.path.dirname(os.path.abspath(__file__)) + + old_names = [(basefold + '/' + "smeared_bao_emu"), + (basefold + '/' + "smeared_bao_emu_1")] + for old_name in old_names: + if os.path.exists(old_name): + import shutil + shutil.rmtree(old_name) + + emulator_name = (basefold + '/' + + "smeared_bao_emu_1.0.0") + + if (not os.path.exists(emulator_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading emulator data (1.2 Mb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/smeared_bao_emu_1.0.0.tar', + emulator_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_name + '.tar') + + customs = { + "accuracy_exp_002": accuracy_exp_002, + "accuracy_exp_005": accuracy_exp_005, + "mean_absolute_exp_percentage_error":mean_absolute_exp_percentage_error + } + metrics_list = ["accuracy",accuracy_exp_002, accuracy_exp_005] + + emulator = {} + emulator['emu_type'] = 'nn' + emulator['model'] = load_model(emulator_name, custom_objects=customs, compile=False) + emulator['model'].compile(optimizer='adam', loss=mean_absolute_exp_percentage_error, metrics=metrics_list) + + file_to_read = open(f"{emulator_name}/details.pickle", "rb") + nn_details = pickle.load(file_to_read) + emulator['k'] = nn_details['kk'] + emulator['scaler'] = nn_details['scaler'] + emulator['bounds'] = nn_details['bounds'] + emulator['keys'] = ['omega_cold', 'sigma8_cold', 'omega_baryon', 'ns', + 'hubble', 'neutrino_mass', 'w0', 'wa', 'expfactor'] + + if verbose: + print('Smeared bao emulator loaded in memory.') + + return emulator + +def load_no_wiggles_emu(verbose=True): + """Loads in memory the no-wiggles emulator. + + :return: a dictionary containing the emulator object + :rtype: dict + """ + + if verbose: + print('Loading no-wiggles emulator...') + import deepdish as dd + + basefold = os.path.dirname(os.path.abspath(__file__)) + + emulator_name = (basefold + '/' + + "no_wiggles_emu_1.0.0") + + if (not os.path.exists(emulator_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading emulator data (6.7 Mb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/no_wiggles_emu_1.0.0.tar', + emulator_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_name + '.tar') + + customs = { + "accuracy_exp_002": accuracy_exp_002, + "accuracy_exp_005": accuracy_exp_005, + "mean_absolute_exp_percentage_error":mean_absolute_exp_percentage_error + } + metrics_list = ["accuracy",accuracy_exp_002, accuracy_exp_005] + + emulator = {} + emulator['emu_type'] = 'nn' + emulator['model'] = load_model(emulator_name, custom_objects=customs, compile=False) + emulator['model'].compile(optimizer='adam', loss=mean_absolute_exp_percentage_error, metrics=metrics_list) + + file_to_read = open(f"{emulator_name}/details.pickle", "rb") + nn_details = pickle.load(file_to_read) + + emulator['k'] = nn_details['kk'] + emulator['scaler'] = nn_details['scaler'] + emulator['bounds'] = nn_details['bounds'] + emulator['keys'] = nn_details['keys'] + + if verbose: + print('No-wiggles emulator loaded in memory.') + + return emulator + +def load_nonlinear_emu(verbose=True, fold_name=None, detail_name=None): + """Loads in memory the nonlinear emulator described in Angulo et al. 2020. + + :return: a dictionary containing the emulator object + :rtype: dict + """ + if verbose: + print('Loading non-linear emulator...') + + if fold_name is None: + basefold = os.path.dirname(os.path.abspath(__file__)) + + old_emulator_names = [(basefold + '/' + + "NN_emulator_data_iter4_big_160.pickle_sg_0.95_2000_rot_bao"), + (basefold + '/' + + "NN_emulator_data_iter4_big_160.pickle_sg_0.99_2000_PCA5_BNFalse_DO0rot_bao"), + (basefold + '/' + + "NN_emulator_data_iter4_big_160.pickle_sg_0.99_2000_PCA6_BNFalse_DO0rot_bao")] + for old_emulator_name in old_emulator_names: + if os.path.exists(old_emulator_name): + import shutil + shutil.rmtree(old_emulator_name) + + emulator_name = (basefold + '/' + + "nonlinear_emu_1.0.1") + detail_name = 'details.pickle' + + if (not os.path.exists(emulator_name)): + import urllib.request + import tarfile + import ssl + ssl._create_default_https_context = ssl._create_unverified_context + print('Downloading emulator data (3Mb)...') + urllib.request.urlretrieve( + 'https://bacco.dipc.org/nonlinear_emu_1.0.1.tar', + emulator_name + '.tar', + MyProgressBar()) + tf = tarfile.open(emulator_name+'.tar', 'r') + tf.extractall(path=basefold) + tf.close() + os.remove(emulator_name + '.tar') + else: + emulator_name = fold_name + emulator = {} + emulator['emu_type'] = 'nn' + emulator['model'] = load_model(emulator_name, compile=False) + with open(os.path.join(emulator_name, detail_name), 'rb') as f: + emulator['scaler'] = pickle.load(f) + _pca = pickle.load(f) + emulator['k'] = pickle.load(f) + _trcoords = pickle.load(f) + emulator['bounds'] = pickle.load(f) + + emulator['keys'] = ['omega_cold', 'sigma8_cold', 'omega_baryon', 'ns', + 'hubble', 'neutrino_mass', 'w0', 'wa', 'expfactor'] + if verbose: + print('Nonlinear emulator loaded in memory.') + return emulator + +def _compute_camb_spectrum(params, kmax=50, k_per_logint=0, cold=True): + """ + Calls camb with the current cosmological parameters and returns a + dictionary with the following keys: + kk, pk + :param cold: whether to return the cold matter power spectrum or the total one. Default to cold. + :type cold: bool, optional + + Through the species keyword the following power spectra can be obtained: + matter, cdm, baryons, neutrinos, cold matter (cdm+baryons), photons, + divergence of the cdm velocity field, divergence of the baryon velocity + field, divergence of the cdm-baryon relative velocity field + """ + import camb + + if 'tau' not in params.keys(): + params['tau'] = 0.0952 + if 'num_massive_neutrinos' not in params.keys(): + params['num_massive_neutrinos'] = 3 if params['neutrino_mass'] != 0.0 else 0 + if 'Neffective' not in params.keys(): + params['Neffective'] = 3.046 + if 'omega_k' not in params.keys(): + params['omega_k'] = 0 + if 'omega_cdm' not in params.keys(): + if 'omega_cold' in params.keys(): + params['omega_cdm'] = params['omega_cold'] - params['omega_baryon'] + elif 'omega_matter' in params.keys(): + params['omega_cdm'] = params['omega_matter'] - params['omega_baryon'] - params['neutrino_mass'] / 93.14 / params['hubble']**2 + else: + raise ValueError('At least one among omega_matter and omega_cold should be specified') + + assert params['omega_k'] == 0, 'Non flat geometries are not supported' + + expfactor = params['expfactor'] + + # Set up a new set of parameters for CAMB + pars = camb.CAMBparams() + + # This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency + # Set neutrino-related parameters + # camb.nonlinear.Halofit('takahashi') + pars.set_cosmology( + H0=100 * params['hubble'], + ombh2=(params['omega_baryon'] * params['hubble']**2), + omch2=(params['omega_cdm'] * params['hubble']**2), + omk=params['omega_k'], + neutrino_hierarchy='degenerate', + num_massive_neutrinos=params['num_massive_neutrinos'], + mnu=params['neutrino_mass'], + standard_neutrino_neff=params['Neffective'], + tau=params['tau']) + + + if 'A_s' in params.keys(): + if params['A_s'] is not None: + A_s = params['A_s'] + ReNormalizeInputSpectrum = False + else: + A_s = 2.e-9 + ReNormalizeInputSpectrum = True + else: + A_s = 2.e-9 + ReNormalizeInputSpectrum = True + + pars.set_dark_energy( + w=params['w0'], + wa=params['wa']) + + redshifts = [(1 / expfactor - 1)] + if expfactor < 1: + redshifts.append(0) + + pars.InitPower.set_params(ns=params['ns'], As=A_s) + pars.YHe = 0.24 + pars.set_matter_power( + redshifts=redshifts, + kmax=kmax, + k_per_logint=k_per_logint) + + pars.WantCls = False + pars.WantScalars = False + pars.Want_CMB = False + pars.DoLensing = False + + # calculate results for these parameters + results = camb.get_results(pars) + + if cold: + index = 7 # cdm + baryons + else: + index = 6 + kh, z, pk = results.get_linear_matter_power_spectrum(var1=(1 + index), + var2=(1 + index), + hubble_units=True, + have_power_spectra=False, + params=None) + pk = pk[-1, :] + + if ReNormalizeInputSpectrum: + sigma8 = results.get_sigmaR(8, z_indices=-1, var1=(1 + index), var2=(1 + index)) + if cold: + Normalization = (params['sigma8_cold'] / sigma8)**2 + else: + Normalization = (params['sigma8'] / sigma8)**2 + pk *= Normalization + + mask = (kh > 1e-4) + + return {'k': kh[mask], 'pk': pk[mask]} + +def compute_camb_pk(coordinates, k=None, cold=True): + """Compute the linear prediction of the matter power spectrum using camb + + The coordinates must be specified as a dictionary with the following + keywords: + + #. 'omega_cold' + #. 'omega_baryon' + #. 'sigma8' + #. 'hubble' + #. 'ns' + #. 'neutrino_mass' + #. 'w0' + #. 'wa' + #. 'expfactor' + + :param coordinates: a set of coordinates in parameter space + :type coordinates: dict + :param k: a vector of wavemodes in h/Mpc, if None the wavemodes used by + camb are returned, defaults to None + :type k: array_like, optional + :param cold: whether to return the cold matter power spectrum or the total one. Default to cold. + :type cold: bool, optional + :return: k and linear pk + :rtype: tuple + """ + _pk_dict = _compute_camb_spectrum(coordinates, cold=cold) + if k is not None: + _k = k + _interp = interpolate.interp1d(np.log(_pk_dict['k']), np.log(_pk_dict['pk']), kind='cubic') + _pk = np.exp(_interp(np.log(_k))) + else: + _k = _pk_dict['k'] + _pk = _pk_dict['pk'] + return _k, _pk + +def _nowiggles_pk(k_lin=None, pk_lin=None, k_emu=None): + """De-wiggled linear prediction of the cold matter power spectrum + + The BAO feature is removed by identifying and removing its corresponding + bump in real space, by means of a DST, and consequently transforming + back to Fourier space. + See: + - Baumann et al 2018 (https://arxiv.org/pdf/1712.08067.pdf) + - Giblin et al 2019 (https://arxiv.org/pdf/1906.02742.pdf) + + :param k_lin: a vector of wavemodes in h/Mpc, if None the wavemodes used by + camb are returned, defaults to None + :type k_lin: array_like, optional + :param pk_lin: a vector of linear power spectrum computed at k_lin, if None + camb will be called, defaults to None + :type pk_lin: array_like, optional + + :param k_emu: a vector of wavemodes in h/Mpc, if None the wavemodes used by + the emulator are returned, defaults to None + :type k_emu: array_like, optional + + :return: dewiggled pk computed at k_emu + :rtype: array_like + """ + + from scipy.fftpack import dst, idst + + nk = int(2**15) + kmin = k_lin.min() + kmax = 10 + klin = np.linspace(kmin, kmax, nk) + + pkcamb_cs = interpolate.splrep(np.log(k_lin), np.log(pk_lin), s=0) + pklin = np.exp(interpolate.splev(np.log(klin), pkcamb_cs, der=0, ext=0)) + + f = np.log10(klin * pklin) + + dstpk = dst(f, type=2) + + even = dstpk[0::2] + odd = dstpk[1::2] + + i_even = np.arange(len(even)).astype(int) + i_odd = np.arange(len(odd)).astype(int) + + even_cs = interpolate.splrep(i_even, even, s=0) + odd_cs = interpolate.splrep(i_odd, odd, s=0) + + even_2nd_der = interpolate.splev(i_even, even_cs, der=2, ext=0) + odd_2nd_der = interpolate.splev(i_odd, odd_cs, der=2, ext=0) + + # these indexes have been fudged for the k-range considered + # [~1e-4, 10], any other choice would require visual inspection + imin_even = i_even[100:300][np.argmax(even_2nd_der[100:300])] - 20 + imax_even = i_even[100:300][np.argmin(even_2nd_der[100:300])] + 70 + imin_odd = i_odd[100:300][np.argmax(odd_2nd_der[100:300])] - 20 + imax_odd = i_odd[100:300][np.argmin(odd_2nd_der[100:300])] + 75 + + i_even_holed = np.concatenate((i_even[:imin_even], i_even[imax_even:])) + i_odd_holed = np.concatenate((i_odd[:imin_odd], i_odd[imax_odd:])) + + even_holed = np.concatenate((even[:imin_even], even[imax_even:])) + odd_holed = np.concatenate((odd[:imin_odd], odd[imax_odd:])) + + even_holed_cs = interpolate.splrep(i_even_holed, even_holed * (i_even_holed+1)**2, s=0) + odd_holed_cs = interpolate.splrep(i_odd_holed, odd_holed * (i_odd_holed+1)**2, s=0) + + even_smooth = interpolate.splev(i_even, even_holed_cs, der=0, ext=0) / (i_even + 1)**2 + odd_smooth = interpolate.splev(i_odd, odd_holed_cs, der=0, ext=0) / (i_odd + 1)**2 + + dstpk_smooth = [] + for ii in range(len(i_even)): + dstpk_smooth.append(even_smooth[ii]) + dstpk_smooth.append(odd_smooth[ii]) + dstpk_smooth = np.array(dstpk_smooth) + + pksmooth = idst(dstpk_smooth, type=2) / (2 * len(dstpk_smooth)) + pksmooth = 10**(pksmooth) / klin + + k_highk = k_lin[k_lin > 5] + p_highk = pk_lin[k_lin > 5] + + k_extended = np.concatenate((klin[klin < 5], k_highk)) + p_extended = np.concatenate((pksmooth[klin < 5], p_highk)) + + pksmooth_cs = interpolate.splrep(np.log(k_extended), np.log(p_extended), s=0) + pksmooth_interp = np.exp(interpolate.splev(np.log(k_emu), pksmooth_cs, der=0, ext=0)) + + return pksmooth_interp + +def _smeared_bao_pk(k_lin=None, pk_lin=None, k_emu=None, pk_lin_emu=None, pk_nw=None, grid=None): + """Prediction of the cold matter power spectrum using a Boltzmann solver with smeared BAO feature + + :param k_lin: a vector of wavemodes in h/Mpc, if None the wavemodes used by + camb are returned, defaults to None + :type k_lin: array_like, optional + :param pk_lin: a vector of linear power spectrum computed at k_lin, if None + camb will be called, defaults to None + :type pk_lin: array_like, optional + + :param k_emu: a vector of wavemodes in h/Mpc, if None the wavemodes used by + the emulator are returned, defaults to None + :type k_emu: array_like, optional + :param pk_emu: a vector of linear power spectrum computed at k_emu, defaults to None + :type pk_emu: array_like, optional + :param pk_nw: a vector of no-wiggles power spectrum computed at k_emu, defaults to None + :type pk_nw: array_like, optional + :param grid: dictionary with parameter and vector of values where to evaluate the emulator, defaults to None + :type grid: array_like, optional + + :return: smeared BAO pk computed at k_emu + :rtype: array_like + """ + from scipy.integrate import trapz + + if grid is None: + sigma_star_2 = trapz(k_lin * pk_lin, x=np.log(k_lin)) / (3 * np.pi**2) + k_star_2 = 1 / sigma_star_2 + G = np.exp(-0.5 * (k_emu**2 / k_star_2)) + if pk_nw is None: + pk_nw = _nowiggles_pk(k_lin=k_lin, pk_lin=pk_lin, k_emu=k_emu) + else: + sigma_star_2 = trapz(k_lin[None,:] * pk_lin, x=np.log(k_lin[None:,]), axis=1) / (3 * np.pi**2) + k_star_2 = 1 / sigma_star_2 + G = np.exp(-0.5 * (k_emu**2 / k_star_2[:,None])) + if pk_nw is None: + pk_nw = np.array([_nowiggles_pk(k_lin=k_lin, pk_lin=pkl, k_emu=k_emu) for pkl in pk_lin]) + return pk_lin_emu * G + pk_nw * (1 - G) diff --git a/structure/baccoemu/baccoemu_vendored/tqdm.patch b/structure/baccoemu/baccoemu_vendored/tqdm.patch new file mode 100644 index 00000000..f612b6e1 --- /dev/null +++ b/structure/baccoemu/baccoemu_vendored/tqdm.patch @@ -0,0 +1,41 @@ +diff --git a/structure/baccoemu/baccoemu_vendored/__init__.py b/structure/baccoemu/baccoemu_vendored/__init__.py +index 85b3d19..18c66e8 100644 +--- a/structure/baccoemu/baccoemu_vendored/__init__.py ++++ b/structure/baccoemu/baccoemu_vendored/__init__.py +@@ -1,7 +1,7 @@ + import numpy as np + import copy + import pickle +-import progressbar ++import tqdm + import hashlib + from ._version import __version__ + from .utils import * +diff --git a/structure/baccoemu/baccoemu_vendored/utils.py b/structure/baccoemu/baccoemu_vendored/utils.py +index 9b448a1..4e6cd65 100644 +--- a/structure/baccoemu/baccoemu_vendored/utils.py ++++ b/structure/baccoemu/baccoemu_vendored/utils.py +@@ -1,7 +1,7 @@ + import numpy as np + import copy + import pickle +-import progressbar ++import tqdm + import hashlib + from ._version import __version__ + +@@ -68,11 +68,10 @@ class MyProgressBar(): + + def __call__(self, block_num, block_size, total_size): + if not self.pbar: +- self.pbar=progressbar.ProgressBar(maxval=total_size) +- self.pbar.start() +- ++ self.pbar=tqdm.trange(total_size) ++ self.pbar.set_description("Downloading") + downloaded = block_num * block_size + if downloaded < total_size: + self.pbar.update(downloaded) + else: +- self.pbar.finish() ++ self.pbar.close() diff --git a/structure/baccoemu/baccoemu_vendored/utils.py b/structure/baccoemu/baccoemu_vendored/utils.py new file mode 100644 index 00000000..4e6cd65e --- /dev/null +++ b/structure/baccoemu/baccoemu_vendored/utils.py @@ -0,0 +1,77 @@ +import numpy as np +import copy +import pickle +import tqdm +import hashlib +from ._version import __version__ + +def _md5(fname): + hash_md5 = hashlib.md5() + with open(fname, "rb") as f: + for chunk in iter(lambda: f.read(4096), b""): + hash_md5.update(chunk) + return hash_md5.hexdigest() + +def _transform_space(x, space_rotation=False, rotation=None, bounds=None): + """Normalize coordinates to [0,1] intervals and if necessary apply a rotation + + :param x: coordinates in parameter space + :type x: ndarray + :param space_rotation: whether to apply the rotation matrix defined through + the rotation keyword, defaults to False + :type space_rotation: bool, optional + :param rotation: rotation matrix, defaults to None + :type rotation: ndarray, optional + :param bounds: ranges within which the emulator hypervolume is defined, + defaults to None + :type bounds: ndarray, optional + :return: normalized and (if required) rotated coordinates + :rtype: ndarray + """ + if space_rotation: + #Get x into the eigenbasis + R = rotation['rotation_matrix'].T + xR = copy.deepcopy(np.array([np.dot(R, xi) + for xi in x])) + xR = xR - rotation['rot_points_means'] + xR = xR/rotation['rot_points_stddevs'] + return xR + else: + return (x - bounds[:, 0])/(bounds[:, 1] - bounds[:, 0]) + +def accuracy_exp_002(y_true, y_pred): + dataset = K.abs(K.exp(y_pred)/K.exp(y_true)-1) + tot = dataset >= 0 + sel = dataset <= 0.002 + return K.shape(dataset[sel])[0] /K.shape(dataset[tot])[0] + +def accuracy_exp_005(y_true, y_pred): + dataset = K.abs(K.exp(y_pred)/K.exp(y_true)-1) + tot = dataset >= 0 + sel = dataset <= 0.005 + return K.shape(dataset[sel])[0] /K.shape(dataset[tot])[0] + +def accuracy_exp_01(y_true, y_pred): + dataset = K.abs(K.exp(y_pred)/K.exp(y_true)-1) + tot = dataset >= 0 + sel = dataset <= 0.01 + return K.shape(dataset[sel])[0] /K.shape(dataset[tot])[0] + +def mean_absolute_exp_percentage_error(y_true, y_pred): + diff = K.abs((K.exp(y_true) - K.exp(y_pred)) / K.clip(K.exp(y_true), + K.epsilon(),None)) + return K.mean(diff, axis=-1) + +class MyProgressBar(): + def __init__(self): + self.pbar = None + + def __call__(self, block_num, block_size, total_size): + if not self.pbar: + self.pbar=tqdm.trange(total_size) + self.pbar.set_description("Downloading") + downloaded = block_num * block_size + if downloaded < total_size: + self.pbar.update(downloaded) + else: + self.pbar.close() diff --git a/structure/baccoemu/module.yaml b/structure/baccoemu/module.yaml new file mode 100644 index 00000000..2b147db6 --- /dev/null +++ b/structure/baccoemu/module.yaml @@ -0,0 +1,139 @@ +#This is a template for module description files +name: "bacco_emulator" +version: "7e0ca8b556da6ad8e11168026078b1d29920adcf" +purpose: "Emulate the non-linear, baryonified, matter power spectrum" +url: "https://bitbucket.org/rangulo/baccoemu" +interface: "baccoemu_interface.py" +attribution: [Giovanni Aricò, Raul E. Angulo, Sergio Contreras, Lurdes Ondaro-Mallea, Marcos Pellejero-Ibañez, Matteo Zennaro, Jens Stücker, Simon Samuroff] +rules: + Please cite the papers below if you use this code in your research" +cite: + - "https://doi.org/10.1093/mnras/stab1911" + - "https://doi.org/10.1093/mnras/stab2018" + - "https://doi.org/10.12688/openreseurope.14310.2" + +assumptions: + - "Neural network emulation of NL baryonic power effcts" + - "w0waCDM cosmology" + +explanation: | + Baccoemu is a collection of cosmological neural-network emulators for large-scale structure statistics + in a wide cosmological parameter space, including dynamical dark energy and massive neutrinos. + These emulators were developed as part of the Bacco project. + + We imported the Bacco Emulator code into this directory because the PyPI version of the code is not recent, + and it has some dependencies that we would like to install from conda-forge rather than PyPI (tensorflow). + We have to make one change to it to make it work on python 3.8 and 3.9 - replacing the progressbar2 + library with TQDM. A patch file with this change is included. + + The license for this package says copyright belongs to the python packaging authority, but that seems unlikely. + +# List of parameters that can go in the params.ini file in the section for this module +params: + mode: + meaning: | + What to emulate; the NL DM-only power, the baryonic boost, or both. + Choose from 'nonlinear', 'baryons', 'nonlinear+baryons'. If 'baryons' is chosen then + the NL power is read from the block (e.g. from camb). Otherwise the NL power is + emulated. + type: str + default: pk_nl_only + +#Inputs for a given choice of a parameter, from the values.ini or from other modules +#If no such choices, just do one of these omitting mode=something part: +inputs: + matter_power_lin: + z: + meaning: Redshifts of samples. + type: real 1d + k_h: + meaning: Wavenumbers k of samples in Mpc/h. + type: real 1d + p_k: + meaning: Linear power spectrum at samples in (Mpc/h)^-3. + type: real 2d + matter_power_nl: + z: + meaning: Redshifts of samples. Only if mode = = "baryons". + type: real 1d + k_h: + meaning: Wavenumbers k of samples in Mpc/h. Only if mode = = "baryons". + type: real 1d + p_k: + meaning: Linear power spectrum at samples in (Mpc/h)^-3. Only if mode = = "baryons". + type: real 2d + cosmological_parameters: + A_s: + meaning: Amplitude of the primordial power spectrum. + type: real + default: + omega_c: + meaning: Cold dark matter density parameter + type: real + default: + omega_b: + meaning: Baryon density parameter + type: real + default: + n_s: + meaning: Primordial power spectrum spectral index. + type: real + default: + h0: + meaning: Hubble parameter. + type: real + default: + mnu: + meaning: Sum of neutrino masses in eV. + type: real + default: + w0: + meaning: Dark energy equation of state parameter at z=0 + type: real + default: -1.0 + wa: + meaning: Dark energy equation of state rate of change with scale factor + type: real + default: 0.0 + baryon_parameters: + M_c: + meaning: Gas mass parameter; log scale in Msun/h. Range 9.0 to 15.0 + type: real + default: + eta: + meaning: Ejected gas density exponential cut-off scale, in range -0.698 to 0.698 + type: real + default: + beta: + meaning: Gas mass index parameter; log scale, in range -1.0 to 0.698 + type: real + default: + M1_z0_cen: + meaning: characteristic halo mass scale for central galaxies, in range 9.0 to 13.0; log scale in Msun/h + type: real + default: + theta_out: + meaning: Scaling of r200 to give outer power-law profile scale of hot gas radius. Range 0.0 to 0.477 + type: real + default: + theta_inn: + meaning: Scaling of r200 to give inner power-law profile scale of hot gas radius. Range -2.0 to -0.522 + type: real + default: + M_inn: + meaning: Characteristic mass of inner hot gas; log scale in Msun/h. Range 9.0 to 13.5. + type: real + default: + + +outputs: + matter_power_nl: + z: + meaning: Redshifts of samples. + type: real 1d + k_h: + meaning: Wavenumbers k of samples in Mpc/h. + type: real 1d + p_k: + meaning: Non-linear power spectrum at samples in (Mpc/h)^-3, potentially including baryon corrections. + type: real 2d diff --git a/tests/test_cosmosis_standard_library.py b/tests/test_cosmosis_standard_library.py index b215a188..0a4c0e87 100644 --- a/tests/test_cosmosis_standard_library.py +++ b/tests/test_cosmosis_standard_library.py @@ -133,3 +133,21 @@ def test_kids(capsys): run_cosmosis("examples/kids-1000.ini") check_likelihood(capsys, "-47.6") check_no_camb_warnings(capsys) + +def test_bacco(): + # The baseline version just does non-linear matter power + run_cosmosis("examples/bacco.ini") + + # This variant emulates NL power with baryonic effects + run_cosmosis("examples/bacco.ini", + override={ + ("bacco_emulator", "mode"): "nonlinear+baryons", + }) + + # This variant uses camb to get the NL power and only emulates the baryonic effects + run_cosmosis("examples/bacco.ini", + override={ + ("bacco_emulator", "mode"): "baryons", + ("camb", "nonlinear"): "pk", + ("camb", "halofit_version"): "takahashi", + })