Add Bacco emulator

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

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

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

likelihood/planck2018/.gitignore

@@ -2,3 +2,4 @@ data/
 plc-3.0/buildir
 plc-3.0/include
 plc-3.0/share
+baseline

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. 
+

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

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_*

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')

structure/baccoemu/baccoemu_vendored/_version.py

@@ -0,0 +1 @@
+__version__ = "2.1.0"