From 2874450452d77649340d8b2e47a09ff6b291695b Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 19 Oct 2021 17:09:03 +0100 Subject: [PATCH 01/29] initial commit of bias theory class --- soliket/bias.py | 56 ++++++++++++++++++++++++++++++++++++++ soliket/tests/test_bias.py | 29 ++++++++++++++++++++ 2 files changed, 85 insertions(+) create mode 100644 soliket/bias.py create mode 100644 soliket/tests/test_bias.py diff --git a/soliket/bias.py b/soliket/bias.py new file mode 100644 index 00000000..2edd9fd5 --- /dev/null +++ b/soliket/bias.py @@ -0,0 +1,56 @@ +r"""Class to calculate bias models for haloes and galaxies as cobaya Theory classes + +""" + +import numpy as np +from typing import Sequence, Union +from cobaya.theory import Theory + +class Bias(Theory): + + def initialize(self): + kmax: float = 0 # Maximum k (1/Mpc units) for Pk, or zero if not needed + nonlinear: bool = False # whether to get non-linear Pk from CAMB/Class + z: Union[Sequence, np.ndarray] = [] # redshift sampling + + # def get_requirements(self): + + def must_provide(self, **requirements): + + self.kmax = max(self.kmax, options.get('kmax', self.kmax)) + self.z = np.unique(np.concatenate( + (np.atleast_1d(options.get("z", self._default_z_sampling)), + np.atleast_1d(self.z)))) + + # Dictionary of the things needed from CAMB/CLASS + needs = {} + + needs['Pk_grid'] = { + 'vars_pairs': self._var_pairs or [('delta_tot', 'delta_tot')], + 'nonlinear': (True, False) if self.nonlinear else False, + 'z': self.z, + 'k_max': self.kmax + } + + return needs + + # def calculate(self, state, want_derived=True, **params_values_dict): + + # def get_Pkgg(self): + + # def get_Pkgm(self): + +class Linear_bias(Bias): + + params = {'b_lin' : None} + + def calculate(self, state, want_derived=True, **params_values_dict): + + state['Pk_gg'] = self.b_lin * self.b_lin * Pk_grid + state['Pk_gm'] = self.b_lin * Pk_grid + + def get_Pkgg(self): + return self._current_state['Pk_gg'] + + def get_Pkgm(self): + return self._current_state['Pk_gm'] diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py new file mode 100644 index 00000000..e83e7baa --- /dev/null +++ b/soliket/tests/test_bias.py @@ -0,0 +1,29 @@ +# pytest -k bias -v . + +import pytest +import numpy as np + +from cobaya.model import get_model + +def test_bias_import(): + from soliket.bias import Bias + +def test_linear_bias_import(): + from soliket.bias import Linear_bias + +def test_linear_bias(): + + from soliket.bias import Linear_bias + + info = {"params": {"b_lin": 1., + "ombh2": 0.0245, + "omch2": 0.1225, + "tau": 0.05}, + "likelihood": {"one" : None}, + "theory": {#"camb" : None, + "linear_bias": {"external": Linear_bias} + }, + "sampler": {"evaluate": None} + } + + model = get_model(info) \ No newline at end of file From 4d2bd56155a8bbf80b3a6c07290e9d7c6ff6fbaf Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Fri, 29 Oct 2021 10:53:40 +0100 Subject: [PATCH 02/29] linear bias model and testing --- soliket/bias.py | 60 +++++++++++++++++++++++++++--------- soliket/tests/test_bias.py | 49 ++++++++++++++++++++++++----- soliket/tests/test_bias.yaml | 10 ++++++ 3 files changed, 98 insertions(+), 21 deletions(-) create mode 100644 soliket/tests/test_bias.yaml diff --git a/soliket/bias.py b/soliket/bias.py index 2edd9fd5..16521945 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -6,17 +6,30 @@ from typing import Sequence, Union from cobaya.theory import Theory + class Bias(Theory): + kmax: float = 10. # Maximum k (1/Mpc units) for Pk, or zero if not needed + nonlinear: bool = False # whether to get non-linear Pk from CAMB/Class + z: Union[Sequence, np.ndarray] = [] # redshift sampling + + extra_args: dict = {} # extra (non-parameter) arguments passed to ccl.Cosmology() + + _logz = np.linspace(-3, np.log10(1100), 150) + _default_z_sampling = 10**_logz + _default_z_sampling[0] = 0 + def initialize(self): - kmax: float = 0 # Maximum k (1/Mpc units) for Pk, or zero if not needed - nonlinear: bool = False # whether to get non-linear Pk from CAMB/Class - z: Union[Sequence, np.ndarray] = [] # redshift sampling - # def get_requirements(self): + self._var_pairs = set() + + def get_requirements(self): + return {} def must_provide(self, **requirements): + options = requirements.get('linear_bias') or {} + self.kmax = max(self.kmax, options.get('kmax', self.kmax)) self.z = np.unique(np.concatenate( (np.atleast_1d(options.get("z", self._default_z_sampling)), @@ -25,6 +38,11 @@ def must_provide(self, **requirements): # Dictionary of the things needed from CAMB/CLASS needs = {} + self.nonlinear = self.nonlinear or options.get('nonlinear', False) + self._var_pairs.update( + set((x, y) for x, y in + options.get('vars_pairs', [('delta_tot', 'delta_tot')]))) + needs['Pk_grid'] = { 'vars_pairs': self._var_pairs or [('delta_tot', 'delta_tot')], 'nonlinear': (True, False) if self.nonlinear else False, @@ -34,23 +52,37 @@ def must_provide(self, **requirements): return needs - # def calculate(self, state, want_derived=True, **params_values_dict): + def _get_Pk_mm(self): - # def get_Pkgg(self): + for pair in self._var_pairs: - # def get_Pkgm(self): + if self.nonlinear: + k, z, Pk_nonlin = self.provider.get_Pk_grid(var_pair=pair, + nonlinear=True) + Pk_mm = np.flip(Pk_nonlin, axis=0) + else: + k, z, Pk_lin = self.provider.get_Pk_grid(var_pair=pair, + nonlinear=False) + Pk_mm = np.flip(Pk_lin, axis=0) -class Linear_bias(Bias): + return Pk_mm - params = {'b_lin' : None} - - def calculate(self, state, want_derived=True, **params_values_dict): - - state['Pk_gg'] = self.b_lin * self.b_lin * Pk_grid - state['Pk_gm'] = self.b_lin * Pk_grid + # def calculate(self, state, want_derived=True, **params_values_dict): def get_Pkgg(self): return self._current_state['Pk_gg'] def get_Pkgm(self): return self._current_state['Pk_gm'] + + +class Linear_bias(Bias): + + params = {'b_lin': None} + + def calculate(self, state, want_derived=True, **params_values_dict): + + Pk_mm = self._get_Pk_mm() + + state['Pk_gg'] = params_values_dict['b_lin']**2. * Pk_mm + state['Pk_gm'] = params_values_dict['b_lin'] * Pk_mm diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index e83e7baa..718c8eda 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -4,26 +4,61 @@ import numpy as np from cobaya.model import get_model +from cobaya.run import run + def test_bias_import(): from soliket.bias import Bias + def test_linear_bias_import(): from soliket.bias import Linear_bias -def test_linear_bias(): + +def test_linear_bias_model(): + + from soliket.bias import Linear_bias + + info = {"params": { + "b_lin": 1., + "H0": 70., + "ombh2": 0.0245, + "omch2": 0.1225, + "ns": 0.96, + "As": 2.2e-9, + "tau": 0.05 + }, + "likelihood": {"one": None}, + "theory": {"camb": None, + "linear_bias": {"external": Linear_bias} + }, + "sampler": {"evaluate": None}, + "debug": True + } + + model = get_model(info) # noqa F841 + + +def test_linear_bias_run(): from soliket.bias import Linear_bias - info = {"params": {"b_lin": 1., + info = {"params": { + "b_lin": 1., + "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, - "tau": 0.05}, - "likelihood": {"one" : None}, - "theory": {#"camb" : None, + "ns": 0.96, + "As": 2.2e-9, + "tau": 0.05 + }, + "likelihood": {"one": None}, + "theory": {"camb": None, "linear_bias": {"external": Linear_bias} }, - "sampler": {"evaluate": None} + "sampler": {"evaluate": None}, + "debug": True } - model = get_model(info) \ No newline at end of file + model = get_model(info) # noqa F841 + updated_info, sampler = run(info) diff --git a/soliket/tests/test_bias.yaml b/soliket/tests/test_bias.yaml new file mode 100644 index 00000000..796bdfae --- /dev/null +++ b/soliket/tests/test_bias.yaml @@ -0,0 +1,10 @@ +params: + mnu: 0.0 + sigma8: +theory: + camb: + linear_bias +debug: false +stop_at_error: true +sampler: + evaluate: From 04d83c6ed462f0743138ce7e7b27905d09d91855 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Mon, 1 Nov 2021 17:26:36 +0000 Subject: [PATCH 03/29] test can compute model successfully --- soliket/bias.py | 15 +++++++++------ soliket/tests/test_bias.py | 18 +++++++++++++++++- soliket/tests/test_bias.yaml | 1 - 3 files changed, 26 insertions(+), 8 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 16521945..a9bf1743 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -2,6 +2,7 @@ """ +import pdb import numpy as np from typing import Sequence, Union from cobaya.theory import Theory @@ -22,6 +23,7 @@ class Bias(Theory): def initialize(self): self._var_pairs = set() + # self._var_pairs = [('delta_tot', 'delta_tot')] def get_requirements(self): return {} @@ -50,6 +52,7 @@ def must_provide(self, **requirements): 'k_max': self.kmax } + assert len(self._var_pairs) < 2, "CCL doesn't support other Pk yet" return needs def _get_Pk_mm(self): @@ -57,22 +60,22 @@ def _get_Pk_mm(self): for pair in self._var_pairs: if self.nonlinear: - k, z, Pk_nonlin = self.provider.get_Pk_grid(var_pair=pair, + k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, nonlinear=True) - Pk_mm = np.flip(Pk_nonlin, axis=0) + # Pk_mm = np.flip(Pk_nonlin, axis=0) else: - k, z, Pk_lin = self.provider.get_Pk_grid(var_pair=pair, + k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, nonlinear=False) - Pk_mm = np.flip(Pk_lin, axis=0) + # Pk_mm = np.flip(Pk_lin, axis=0) return Pk_mm # def calculate(self, state, want_derived=True, **params_values_dict): - def get_Pkgg(self): + def get_Pk_gg(self): return self._current_state['Pk_gg'] - def get_Pkgm(self): + def get_Pk_gm(self): return self._current_state['Pk_gm'] diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 718c8eda..24732cdc 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -61,4 +61,20 @@ def test_linear_bias_run(): } model = get_model(info) # noqa F841 - updated_info, sampler = run(info) + model.add_requirements({"Pk_grid": {"z": 0., "k_max": 10., + "nonlinear": False, + "vars_pairs": ('delta_tot', 'delta_tot') + }, + "Pk_gg": None + }) + + model.logposterior(info['params']) # force computation of model + + lhood = model.likelihood['one'] + + k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + nonlinear=False) + + Pk_gg = lhood.provider.get_Pk_gg() + + assert np.allclose(Pk_mm_lin, Pk_gg) diff --git a/soliket/tests/test_bias.yaml b/soliket/tests/test_bias.yaml index 796bdfae..e5e00a6f 100644 --- a/soliket/tests/test_bias.yaml +++ b/soliket/tests/test_bias.yaml @@ -3,7 +3,6 @@ params: sigma8: theory: camb: - linear_bias debug: false stop_at_error: true sampler: From 5f867a5c60a394508cece196d6e97c0ff2559397 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Mon, 1 Nov 2021 17:28:56 +0000 Subject: [PATCH 04/29] test Pk_gm too --- soliket/tests/test_bias.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 24732cdc..e51c3aad 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -65,7 +65,8 @@ def test_linear_bias_run(): "nonlinear": False, "vars_pairs": ('delta_tot', 'delta_tot') }, - "Pk_gg": None + "Pk_gg": None, + "Pk_gm": None }) model.logposterior(info['params']) # force computation of model @@ -76,5 +77,7 @@ def test_linear_bias_run(): nonlinear=False) Pk_gg = lhood.provider.get_Pk_gg() + Pk_gm = lhood.provider.get_Pk_gm() - assert np.allclose(Pk_mm_lin, Pk_gg) + assert np.allclose(Pk_mm_lin * info["params"]["b_lin"]**2., Pk_gg) + assert np.allclose(Pk_mm_lin * info["params"]["b_lin"], Pk_gm) From 4bc42ed102d806177069f698138e9942ce799338 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Mon, 1 Nov 2021 17:30:59 +0000 Subject: [PATCH 05/29] change test name --- soliket/tests/test_bias.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index e51c3aad..ef8de0ab 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -39,7 +39,7 @@ def test_linear_bias_model(): model = get_model(info) # noqa F841 -def test_linear_bias_run(): +def test_linear_bias_compute(): from soliket.bias import Linear_bias From 27643f2f6d07b8d2a7bbd394506e91f809411867 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Mon, 1 Nov 2021 17:38:08 +0000 Subject: [PATCH 06/29] rename to grid --- soliket/bias.py | 12 ++++++------ soliket/tests/test_bias.py | 10 +++++----- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index a9bf1743..9eab0289 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -72,11 +72,11 @@ def _get_Pk_mm(self): # def calculate(self, state, want_derived=True, **params_values_dict): - def get_Pk_gg(self): - return self._current_state['Pk_gg'] + def get_Pk_gg_grid(self): + return self._current_state['Pk_gg_grid'] - def get_Pk_gm(self): - return self._current_state['Pk_gm'] + def get_Pk_gm_grid(self): + return self._current_state['Pk_gm_grid'] class Linear_bias(Bias): @@ -87,5 +87,5 @@ def calculate(self, state, want_derived=True, **params_values_dict): Pk_mm = self._get_Pk_mm() - state['Pk_gg'] = params_values_dict['b_lin']**2. * Pk_mm - state['Pk_gm'] = params_values_dict['b_lin'] * Pk_mm + state['Pk_gg_grid'] = params_values_dict['b_lin']**2. * Pk_mm + state['Pk_gm_grid'] = params_values_dict['b_lin'] * Pk_mm diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index ef8de0ab..b6e9f52c 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -39,7 +39,7 @@ def test_linear_bias_model(): model = get_model(info) # noqa F841 -def test_linear_bias_compute(): +def test_linear_bias_compute_grid(): from soliket.bias import Linear_bias @@ -65,8 +65,8 @@ def test_linear_bias_compute(): "nonlinear": False, "vars_pairs": ('delta_tot', 'delta_tot') }, - "Pk_gg": None, - "Pk_gm": None + "Pk_gg_grid": None, + "Pk_gm_grid": None }) model.logposterior(info['params']) # force computation of model @@ -76,8 +76,8 @@ def test_linear_bias_compute(): k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), nonlinear=False) - Pk_gg = lhood.provider.get_Pk_gg() - Pk_gm = lhood.provider.get_Pk_gm() + Pk_gg = lhood.provider.get_Pk_gg_grid() + Pk_gm = lhood.provider.get_Pk_gm_grid() assert np.allclose(Pk_mm_lin * info["params"]["b_lin"]**2., Pk_gg) assert np.allclose(Pk_mm_lin * info["params"]["b_lin"], Pk_gm) From 8896b5230c4d24eb78afd9633aa8d21155a933e6 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Mon, 1 Nov 2021 18:12:10 +0000 Subject: [PATCH 07/29] added interpolators (need tidying) --- soliket/bias.py | 55 +++++++++++++++++++++++++++++++++++--- soliket/tests/test_bias.py | 51 +++++++++++++++++++++++++++++++++-- 2 files changed, 101 insertions(+), 5 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 9eab0289..aa2d5b73 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -6,6 +6,7 @@ import numpy as np from typing import Sequence, Union from cobaya.theory import Theory +from cobaya.theories.cosmo.boltzmannbase import PowerSpectrumInterpolator class Bias(Theory): @@ -52,7 +53,7 @@ def must_provide(self, **requirements): 'k_max': self.kmax } - assert len(self._var_pairs) < 2, "CCL doesn't support other Pk yet" + assert len(self._var_pairs) < 2, "Bias doesn't support other Pk yet" return needs def _get_Pk_mm(self): @@ -60,11 +61,11 @@ def _get_Pk_mm(self): for pair in self._var_pairs: if self.nonlinear: - k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, + self.k, self.z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, nonlinear=True) # Pk_mm = np.flip(Pk_nonlin, axis=0) else: - k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, + self.k, self.z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, nonlinear=False) # Pk_mm = np.flip(Pk_lin, axis=0) @@ -72,6 +73,54 @@ def _get_Pk_mm(self): # def calculate(self, state, want_derived=True, **params_values_dict): + def get_Pk_gg_interpolator(self, extrap_kmax=None): + + key = ("Pk_gg_interpolator", + self.nonlinear, extrap_kmax) + tuple(sorted(self._var_pairs)) + + pk = self.get_Pk_gg_grid() + log_p = True + sign = 1 + if np.any(pk < 0): + if np.all(pk < 0): + sign = -1 + else: + log_p = False + if log_p: + pk = np.log(sign * pk) + # elif extrap_kmax > k[-1]: + # raise LoggedError(self.log, + # 'Cannot do log extrapolation with zero-crossing pk ' + # 'for %s, %s' % var_pair) + result = PowerSpectrumInterpolator(self.z, self.k, pk, logP=log_p, logsign=sign, + extrap_kmax=extrap_kmax) + self.current_state[key] = result + return result + + def get_Pk_gm_interpolator(self, extrap_kmax=None): + + key = ("Pk_gm_interpolator", + self.nonlinear, extrap_kmax) + tuple(sorted(self._var_pairs)) + + pk = self.get_Pk_gm_grid() + log_p = True + sign = 1 + if np.any(pk < 0): + if np.all(pk < 0): + sign = -1 + else: + log_p = False + if log_p: + pk = np.log(sign * pk) + # elif extrap_kmax > k[-1]: + # raise LoggedError(self.log, + # 'Cannot do log extrapolation with zero-crossing pk ' + # 'for %s, %s' % var_pair) + result = PowerSpectrumInterpolator(self.z, self.k, pk, logP=log_p, logsign=sign, + extrap_kmax=extrap_kmax) + self.current_state[key] = result + return result + def get_Pk_gg_grid(self): return self._current_state['Pk_gg_grid'] diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index b6e9f52c..37695bf7 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -20,7 +20,7 @@ def test_linear_bias_model(): from soliket.bias import Linear_bias info = {"params": { - "b_lin": 1., + "b_lin": 1.1, "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, @@ -44,7 +44,7 @@ def test_linear_bias_compute_grid(): from soliket.bias import Linear_bias info = {"params": { - "b_lin": 1., + "b_lin": 1.1, "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, @@ -81,3 +81,50 @@ def test_linear_bias_compute_grid(): assert np.allclose(Pk_mm_lin * info["params"]["b_lin"]**2., Pk_gg) assert np.allclose(Pk_mm_lin * info["params"]["b_lin"], Pk_gm) + + +def test_linear_bias_compute_interpolator(): + + from soliket.bias import Linear_bias + + info = {"params": { + "b_lin": 1.1, + "H0": 70., + "ombh2": 0.0245, + "omch2": 0.1225, + "ns": 0.96, + "As": 2.2e-9, + "tau": 0.05 + }, + "likelihood": {"one": None}, + "theory": {"camb": None, + "linear_bias": {"external": Linear_bias} + }, + "sampler": {"evaluate": None}, + "debug": True + } + + model = get_model(info) # noqa F841 + model.add_requirements({"Pk_grid": {"z": 0., "k_max": 10., + "nonlinear": False, + "vars_pairs": ('delta_tot', 'delta_tot') + }, + "Pk_gg_interpolator": None, + "Pk_gm_interpolator": None + }) + + model.logposterior(info['params']) # force computation of model + + lhood = model.likelihood['one'] + + k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + nonlinear=False) + + Pk_gg = lhood.provider.get_Pk_gg_interpolator().P + Pk_gg = Pk_gg(z, k) + + Pk_gm = lhood.provider.get_Pk_gm_interpolator().P + Pk_gm = Pk_gm(z, k) + + assert np.allclose(Pk_mm_lin * info["params"]["b_lin"]**2., Pk_gg) + assert np.allclose(Pk_mm_lin * info["params"]["b_lin"], Pk_gm) From d7ca6f1a99f4ab854092b0ff88cc37c6460963c1 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Wed, 3 Nov 2021 15:16:38 +0000 Subject: [PATCH 08/29] tidied tests, extrap error check on interpolator --- soliket/bias.py | 29 ++++++++------- soliket/tests/test_bias.py | 74 ++++++++++++++------------------------ 2 files changed, 43 insertions(+), 60 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index aa2d5b73..7ebbd186 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -7,6 +7,7 @@ from typing import Sequence, Union from cobaya.theory import Theory from cobaya.theories.cosmo.boltzmannbase import PowerSpectrumInterpolator +from cobaya.log import LoggedError class Bias(Theory): @@ -73,10 +74,12 @@ def _get_Pk_mm(self): # def calculate(self, state, want_derived=True, **params_values_dict): - def get_Pk_gg_interpolator(self, extrap_kmax=None): + def get_Pk_gg_interpolator(self, + var_pair=("delta_tot", "delta_tot"), + extrap_kmax=None): key = ("Pk_gg_interpolator", - self.nonlinear, extrap_kmax) + tuple(sorted(self._var_pairs)) + self.nonlinear, extrap_kmax) + tuple(sorted(var_pair)) pk = self.get_Pk_gg_grid() log_p = True @@ -88,19 +91,21 @@ def get_Pk_gg_interpolator(self, extrap_kmax=None): log_p = False if log_p: pk = np.log(sign * pk) - # elif extrap_kmax > k[-1]: - # raise LoggedError(self.log, - # 'Cannot do log extrapolation with zero-crossing pk ' - # 'for %s, %s' % var_pair) + elif extrap_kmax > self.k[-1]: + raise LoggedError(self.log, + 'Cannot do log extrapolation with zero-crossing pk ' + 'for %s, %s' % var_pair) result = PowerSpectrumInterpolator(self.z, self.k, pk, logP=log_p, logsign=sign, extrap_kmax=extrap_kmax) self.current_state[key] = result return result - def get_Pk_gm_interpolator(self, extrap_kmax=None): + def get_Pk_gm_interpolator(self, + var_pair=("delta_tot", "delta_tot"), + extrap_kmax=None): key = ("Pk_gm_interpolator", - self.nonlinear, extrap_kmax) + tuple(sorted(self._var_pairs)) + self.nonlinear, extrap_kmax) + tuple(sorted(var_pair)) pk = self.get_Pk_gm_grid() log_p = True @@ -112,10 +117,10 @@ def get_Pk_gm_interpolator(self, extrap_kmax=None): log_p = False if log_p: pk = np.log(sign * pk) - # elif extrap_kmax > k[-1]: - # raise LoggedError(self.log, - # 'Cannot do log extrapolation with zero-crossing pk ' - # 'for %s, %s' % var_pair) + elif extrap_kmax > self.k[-1]: + raise LoggedError(self.log, + 'Cannot do log extrapolation with zero-crossing pk ' + 'for %s, %s' % var_pair) result = PowerSpectrumInterpolator(self.z, self.k, pk, logP=log_p, logsign=sign, extrap_kmax=extrap_kmax) self.current_state[key] = result diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 37695bf7..1dbd94d1 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -6,6 +6,20 @@ from cobaya.model import get_model from cobaya.run import run +info = {"params": { + "b_lin": 1.1, + "H0": 70., + "ombh2": 0.0245, + "omch2": 0.1225, + "ns": 0.96, + "As": 2.2e-9, + "tau": 0.05 + }, + "likelihood": {"one": None}, + "sampler": {"evaluate": None}, + "debug": True + } + def test_bias_import(): from soliket.bias import Bias @@ -19,22 +33,10 @@ def test_linear_bias_model(): from soliket.bias import Linear_bias - info = {"params": { - "b_lin": 1.1, - "H0": 70., - "ombh2": 0.0245, - "omch2": 0.1225, - "ns": 0.96, - "As": 2.2e-9, - "tau": 0.05 - }, - "likelihood": {"one": None}, - "theory": {"camb": None, - "linear_bias": {"external": Linear_bias} - }, - "sampler": {"evaluate": None}, - "debug": True - } + info["theory"] = { + "camb": None, + "linear_bias": {"external": Linear_bias} + } model = get_model(info) # noqa F841 @@ -43,22 +45,10 @@ def test_linear_bias_compute_grid(): from soliket.bias import Linear_bias - info = {"params": { - "b_lin": 1.1, - "H0": 70., - "ombh2": 0.0245, - "omch2": 0.1225, - "ns": 0.96, - "As": 2.2e-9, - "tau": 0.05 - }, - "likelihood": {"one": None}, - "theory": {"camb": None, - "linear_bias": {"external": Linear_bias} - }, - "sampler": {"evaluate": None}, - "debug": True - } + info["theory"] = { + "camb": None, + "linear_bias": {"external": Linear_bias} + } model = get_model(info) # noqa F841 model.add_requirements({"Pk_grid": {"z": 0., "k_max": 10., @@ -87,22 +77,10 @@ def test_linear_bias_compute_interpolator(): from soliket.bias import Linear_bias - info = {"params": { - "b_lin": 1.1, - "H0": 70., - "ombh2": 0.0245, - "omch2": 0.1225, - "ns": 0.96, - "As": 2.2e-9, - "tau": 0.05 - }, - "likelihood": {"one": None}, - "theory": {"camb": None, - "linear_bias": {"external": Linear_bias} - }, - "sampler": {"evaluate": None}, - "debug": True - } + info["theory"] = { + "camb": None, + "linear_bias": {"external": Linear_bias} + } model = get_model(info) # noqa F841 model.add_requirements({"Pk_grid": {"z": 0., "k_max": 10., From a774e0a43e08f842a5346d1c10d52ee6b52a8176 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Mon, 29 Nov 2021 16:50:02 +0000 Subject: [PATCH 09/29] added LPT class --- soliket/bias.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/soliket/bias.py b/soliket/bias.py index 7ebbd186..15003c54 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -143,3 +143,15 @@ def calculate(self, state, want_derived=True, **params_values_dict): state['Pk_gg_grid'] = params_values_dict['b_lin']**2. * Pk_mm state['Pk_gm_grid'] = params_values_dict['b_lin'] * Pk_mm + +class LPT_bias(Bias): + + params = {b1, b2, b3 etc} + + def calculate(self, state, want_derived=True, **params_values_dict): + + Pk_mm = self._get_Pk_mm() + + state['Pk_gg_grid'] = some function of (Pk_mm) + + From 18b42c9ee51cb2c4749615db76a0f0ecc63abfda Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 22 Feb 2022 12:04:35 +0000 Subject: [PATCH 10/29] added velocileptors LPT bias model Pk_gg --- soliket/bias.py | 68 +++++++++++++++++++++++++++++++++--- soliket/tests/test_bias.py | 18 ++++++++++ soliket/tests/test_bias.yaml | 6 ++++ 3 files changed, 88 insertions(+), 4 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 15003c54..28fc2f8c 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -8,6 +8,7 @@ from cobaya.theory import Theory from cobaya.theories.cosmo.boltzmannbase import PowerSpectrumInterpolator from cobaya.log import LoggedError +from velocileptors.EPT.cleft_kexpanded_resummed_fftw import RKECLEFT class Bias(Theory): @@ -57,6 +58,16 @@ def must_provide(self, **requirements): assert len(self._var_pairs) < 2, "Bias doesn't support other Pk yet" return needs + def _get_growth(self): + for pair in self._var_pairs: + + k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, + nonlinear=False) + + assert(z[0] == 0) + + return np.mean(Pk_mm/Pk_mm[:1], axis = -1)**0.5 + def _get_Pk_mm(self): for pair in self._var_pairs: @@ -146,12 +157,61 @@ def calculate(self, state, want_derived=True, **params_values_dict): class LPT_bias(Bias): - params = {b1, b2, b3 etc} + params = {'b11' : None, 'b21' : None, 'bs1' : None, 'b12' : None, 'b22' : None, 'bs2' : None} - def calculate(self, state, want_derived=True, **params_values_dict): + def init_cleft(self): - Pk_mm = self._get_Pk_mm() + self.cleft_obj = RKECLEFT(self.k, self._get_Pk_mm()) + + self.lpt_table = [] + + for D in self._get_growth(): + self.cleft_obj.make_ptable(D=D, kmin=self.k[0], kmax=self.k[-1], nk=self.k.size) + self.lpt_table.append(self.cleft_obj.pktable) - state['Pk_gg_grid'] = some function of (Pk_mm) + self.lpt_table = np.array(self.lpt_table) + def _get_Pk_gg(self, **params_values_dict): + + b11 = params_values_dict['b11'] + b21 = params_values_dict['b21'] + bs1 = params_values_dict['bs1'] + b12 = params_values_dict['b12'] + b22 = params_values_dict['b22'] + bs2 = params_values_dict['bs2'] + + bL11 = b11 - 1 + bL12 = b12 - 1 + + if self.nonlinear: + pgg = (b11 * b12)[:, None] * Pk_mm + else: + Pdmdm = self.lpt_table[:, :, 1] + Pdmd1 = 0.5 * self.lpt_table[:, :, 2] + Pd1d1 = self.lpt_table[:, :, 3] + pgg = (Pdmdm + (bL11 + bL12)[:, None] * Pdmd1 + + (bL11 * bL12)[:, None] * Pd1d1) + + Pdmd2 = 0.5 * self.lpt_table[:, :, 4] + Pd1d2 = 0.5 * self.lpt_table[:, :, 5] + Pd2d2 = self.lpt_table[:, :, 6] * self.wk_low[None, :] + Pdms2 = 0.25 * self.lpt_table[:, :, 7] + Pd1s2 = 0.25 * self.lpt_table[:, :, 8] + Pd2s2 = 0.25 * self.lpt_table[:, :, 9] * self.wk_low[None, :] + Ps2s2 = 0.25 * self.lpt_table[:, :, 10] * self.wk_low[None, :] + + pgg += ((b21 + b22)[:, None] * Pdmd2 + + (bs1 + bs2)[:, None] * Pdms2 + + (bL11 * b22 + bL12 * b21)[:, None] * Pd1d2 + + (bL11 * bs2 + bL12 * bs1)[:, None] * Pd1s2 + + (b21 * b22)[:, None] * Pd2d2 + + (b21 * bs2 + b22 * bs1)[:, None] * Pd2s2 + + (bs1 * bs2)[:, None] * Ps2s2) + + return pgg + + def calculate(self, state, want_derived=True, **params_values_dict): + + Pk_mm = self._get_Pk_mm() + state['Pk_gg_grid'] = self._get_Pk_gg(Pk_mm, params_values_dict) diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 1dbd94d1..9b7e2eac 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -8,6 +8,12 @@ info = {"params": { "b_lin": 1.1, + "b11": 1., + "b21": 1., + "bs1": 1., + "b12": 1., + "b22": 1., + "bs2": 1., "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, @@ -106,3 +112,15 @@ def test_linear_bias_compute_interpolator(): assert np.allclose(Pk_mm_lin * info["params"]["b_lin"]**2., Pk_gg) assert np.allclose(Pk_mm_lin * info["params"]["b_lin"], Pk_gm) + + +def test_LPT_bias_model(): + + from soliket.bias import LPT_bias + + info["theory"] = { + "camb": None, + "lpt_bias": {"external": LPT_bias} + } + + model = get_model(info) # noqa F841 diff --git a/soliket/tests/test_bias.yaml b/soliket/tests/test_bias.yaml index e5e00a6f..12eb3d4a 100644 --- a/soliket/tests/test_bias.yaml +++ b/soliket/tests/test_bias.yaml @@ -1,4 +1,10 @@ params: + b11: 1. + b21: 1. + bs1: 1. + b12: 1. + b22: 1. + bs2: 1. mnu: 0.0 sigma8: theory: From 475e812e31ca9fe4d7cbcdfcb6a8315b1ca4ca34 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 22 Feb 2022 12:29:48 +0000 Subject: [PATCH 11/29] codestyle --- soliket/bias.py | 17 ++++++++++++----- soliket/tests/test_bias.py | 1 + 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 2c8e6349..f0c4d273 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -65,7 +65,7 @@ def _get_growth(self): assert(z[0] == 0) - return np.mean(Pk_mm/Pk_mm[:1], axis = -1)**0.5 + return np.mean(Pk_mm / Pk_mm[:1], axis=-1)**0.5 def _get_Pk_mm(self): @@ -100,9 +100,15 @@ def calculate(self, state, want_derived=True, **params_values_dict): state['Pk_gg_grid'] = params_values_dict['b_lin']**2. * Pk_mm state['Pk_gm_grid'] = params_values_dict['b_lin'] * Pk_mm + class LPT_bias(Bias): - params = {'b11' : None, 'b21' : None, 'bs1' : None, 'b12' : None, 'b22' : None, 'bs2' : None} + params = {'b11': None, + 'b21': None, + 'bs1': None, + 'b12': None, + 'b22': None, + 'bs2': None} def init_cleft(self): @@ -111,12 +117,13 @@ def init_cleft(self): self.lpt_table = [] for D in self._get_growth(): - self.cleft_obj.make_ptable(D=D, kmin=self.k[0], kmax=self.k[-1], nk=self.k.size) + self.cleft_obj.make_ptable(D=D, + kmin=self.k[0], kmax=self.k[-1], nk=self.k.size) self.lpt_table.append(self.cleft_obj.pktable) self.lpt_table = np.array(self.lpt_table) - def _get_Pk_gg(self, **params_values_dict): + def _get_Pk_gg(self, Pk_mm, **params_values_dict): b11 = params_values_dict['b11'] b21 = params_values_dict['b21'] @@ -152,7 +159,7 @@ def _get_Pk_gg(self, **params_values_dict): (b21 * b22)[:, None] * Pd2d2 + (b21 * bs2 + b22 * bs1)[:, None] * Pd2s2 + (bs1 * bs2)[:, None] * Ps2s2) - + return pgg def calculate(self, state, want_derived=True, **params_values_dict): diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 18245b33..c55167b9 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -78,6 +78,7 @@ def test_linear_bias_compute_grid(): assert np.allclose(Pk_mm_lin * info["params"]["b_lin"]**2., Pk_gg) assert np.allclose(Pk_mm_lin * info["params"]["b_lin"], Pk_gm) + def test_LPT_bias_model(): from soliket.bias import LPT_bias From 29fa1d4bd2a29e0b27a8b32e0def96c39bded4ee Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Thu, 24 Feb 2022 10:35:26 +0000 Subject: [PATCH 12/29] added LPT Pgg grid --- soliket/bias.py | 34 +++++++++++++++++++++------------- soliket/tests/test_bias.py | 27 +++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 13 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index f0c4d273..850bb5da 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -112,25 +112,28 @@ class LPT_bias(Bias): def init_cleft(self): - self.cleft_obj = RKECLEFT(self.k, self._get_Pk_mm()) + Pk_mm = self._get_Pk_mm() + + self.cleft_obj = RKECLEFT(self.k, Pk_mm[0], + extrap_min=np.floor(np.log10(self.k[0])), + extrap_max=np.ceil(np.log10(self.k[-1]))) self.lpt_table = [] for D in self._get_growth(): - self.cleft_obj.make_ptable(D=D, - kmin=self.k[0], kmax=self.k[-1], nk=self.k.size) + self.cleft_obj.make_ptable(D=D, kmin=self.k[0], kmax=self.k[-1], nk=self.k.size) self.lpt_table.append(self.cleft_obj.pktable) self.lpt_table = np.array(self.lpt_table) def _get_Pk_gg(self, Pk_mm, **params_values_dict): - b11 = params_values_dict['b11'] - b21 = params_values_dict['b21'] - bs1 = params_values_dict['bs1'] - b12 = params_values_dict['b12'] - b22 = params_values_dict['b22'] - bs2 = params_values_dict['bs2'] + b11 = params_values_dict['b11'] * np.ones_like(Pk_mm) + b21 = params_values_dict['b21'] * np.ones_like(Pk_mm) + bs1 = params_values_dict['bs1'] * np.ones_like(Pk_mm) + b12 = params_values_dict['b12'] * np.ones_like(Pk_mm) + b22 = params_values_dict['b22'] * np.ones_like(Pk_mm) + bs2 = params_values_dict['bs2'] * np.ones_like(Pk_mm) bL11 = b11 - 1 bL12 = b12 - 1 @@ -146,11 +149,11 @@ def _get_Pk_gg(self, Pk_mm, **params_values_dict): Pdmd2 = 0.5 * self.lpt_table[:, :, 4] Pd1d2 = 0.5 * self.lpt_table[:, :, 5] - Pd2d2 = self.lpt_table[:, :, 6] * self.wk_low[None, :] + Pd2d2 = self.lpt_table[:, :, 6]# * self.wk_low[None, :] Pdms2 = 0.25 * self.lpt_table[:, :, 7] Pd1s2 = 0.25 * self.lpt_table[:, :, 8] - Pd2s2 = 0.25 * self.lpt_table[:, :, 9] * self.wk_low[None, :] - Ps2s2 = 0.25 * self.lpt_table[:, :, 10] * self.wk_low[None, :] + Pd2s2 = 0.25 * self.lpt_table[:, :, 9]# * self.wk_low[None, :] + Ps2s2 = 0.25 * self.lpt_table[:, :, 10]# * self.wk_low[None, :] pgg += ((b21 + b22)[:, None] * Pdmd2 + (bs1 + bs2)[:, None] * Pdms2 + @@ -164,6 +167,11 @@ def _get_Pk_gg(self, Pk_mm, **params_values_dict): def calculate(self, state, want_derived=True, **params_values_dict): + try: + self.cleft_obj + except: + self.init_cleft() + Pk_mm = self._get_Pk_mm() - state['Pk_gg_grid'] = self._get_Pk_gg(Pk_mm, params_values_dict) + state['Pk_gg_grid'] = self._get_Pk_gg(Pk_mm, **params_values_dict) diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index c55167b9..76e29fea 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -89,3 +89,30 @@ def test_LPT_bias_model(): } model = get_model(info) # noqa F841 + +def test_LPT_bias_compute_grid(): + + from soliket.bias import LPT_bias + + info["theory"] = { + "camb": None, + "LPT_bias": {"external": LPT_bias} + } + + model = get_model(info) # noqa F841 + model.add_requirements({"Pk_grid": {"z": 0., "k_max": 10., + "nonlinear": False, + "vars_pairs": ('delta_tot', 'delta_tot'), + "extrap_kmin": 1.e-3 + }, + "Pk_gg_grid": None + }) + + model.logposterior(info['params']) # force computation of model + + lhood = model.likelihood['one'] + + k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + nonlinear=False) + + Pk_gg = lhood.provider.get_Pk_gg_grid() From 9c294974c8ed10c9c77a74ae19e080de836d9ba2 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Mon, 28 Feb 2022 16:52:55 +0000 Subject: [PATCH 13/29] interpolation of spectra from ptable --- soliket/bias.py | 58 +++++++++++++++++++++++--------------- soliket/tests/test_bias.py | 8 ++++++ 2 files changed, 44 insertions(+), 22 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 850bb5da..8530dc4d 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -5,6 +5,8 @@ import pdb import numpy as np from typing import Sequence, Union +from scipy.interpolate import interp1d + from cobaya.theory import Theory from cobaya.log import LoggedError from velocileptors.EPT.cleft_kexpanded_resummed_fftw import RKECLEFT @@ -110,10 +112,15 @@ class LPT_bias(Bias): 'b22': None, 'bs2': None} - def init_cleft(self): + def init_cleft(self, k_filter=None): Pk_mm = self._get_Pk_mm() + if k_filter is not None: + self.wk_low = 1-np.exp(-(self.k/k_filter)**2) + else: + self.wk_low = np.ones_like(self.k) + self.cleft_obj = RKECLEFT(self.k, Pk_mm[0], extrap_min=np.floor(np.log10(self.k[0])), extrap_max=np.ceil(np.log10(self.k[-1]))) @@ -128,49 +135,56 @@ def init_cleft(self): def _get_Pk_gg(self, Pk_mm, **params_values_dict): - b11 = params_values_dict['b11'] * np.ones_like(Pk_mm) - b21 = params_values_dict['b21'] * np.ones_like(Pk_mm) - bs1 = params_values_dict['bs1'] * np.ones_like(Pk_mm) - b12 = params_values_dict['b12'] * np.ones_like(Pk_mm) - b22 = params_values_dict['b22'] * np.ones_like(Pk_mm) - bs2 = params_values_dict['bs2'] * np.ones_like(Pk_mm) + b11 = params_values_dict['b11']# * np.ones_like(Pk_mm) + b21 = params_values_dict['b21']# * np.ones_like(Pk_mm) + bs1 = params_values_dict['bs1']# * np.ones_like(Pk_mm) + b12 = params_values_dict['b12']# * np.ones_like(Pk_mm) + b22 = params_values_dict['b22']# * np.ones_like(Pk_mm) + bs2 = params_values_dict['bs2']# * np.ones_like(Pk_mm) bL11 = b11 - 1 bL12 = b12 - 1 if self.nonlinear: - pgg = (b11 * b12)[:, None] * Pk_mm + pgg = (b11 * b12) * Pk_mm else: Pdmdm = self.lpt_table[:, :, 1] Pdmd1 = 0.5 * self.lpt_table[:, :, 2] Pd1d1 = self.lpt_table[:, :, 3] - pgg = (Pdmdm + (bL11 + bL12)[:, None] * Pdmd1 + - (bL11 * bL12)[:, None] * Pd1d1) + pgg = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) Pdmd2 = 0.5 * self.lpt_table[:, :, 4] Pd1d2 = 0.5 * self.lpt_table[:, :, 5] - Pd2d2 = self.lpt_table[:, :, 6]# * self.wk_low[None, :] + Pd2d2 = self.lpt_table[:, :, 6] * self.wk_low[None, :] Pdms2 = 0.25 * self.lpt_table[:, :, 7] Pd1s2 = 0.25 * self.lpt_table[:, :, 8] - Pd2s2 = 0.25 * self.lpt_table[:, :, 9]# * self.wk_low[None, :] - Ps2s2 = 0.25 * self.lpt_table[:, :, 10]# * self.wk_low[None, :] + Pd2s2 = 0.25 * self.lpt_table[:, :, 9] * self.wk_low[None, :] + Ps2s2 = 0.25 * self.lpt_table[:, :, 10] * self.wk_low[None, :] + + pgg += ((b21 + b22) * Pdmd2 + + (bs1 + bs2) * Pdms2 + + (bL11 * b22 + bL12 * b21) * Pd1d2 + + (bL11 * bs2 + bL12 * bs1) * Pd1s2 + + (b21 * b22) * Pd2d2 + + (b21 * bs2 + b22 * bs1) * Pd2s2 + + (bs1 * bs2) * Ps2s2) + + cleft_k = self.lpt_table[:, :, 0] + + # pdb.set_trace() + + # spectra = {r'$(1,1)$':self.lpt_table[:,:,1],r'$(1,b_1)$':0.5*self.lpt_table[:,:,2], r'$(b_1,b_1)$': self.lpt_table[:,:,3],r'$(1,b_2)$':0.5*self.lpt_table[:,:,4], r'$(b_1,b_2)$': 0.5*self.lpt_table[:,:,5], r'$(b_2,b_2)$': self.lpt_table[:,:,6],r'$(1,b_s)$':0.5*self.lpt_table[:,:,7], r'$(b_1,b_s)$': 0.5*self.lpt_table[:,:,8], r'$(b_2,b_s)$':0.5*self.lpt_table[:,:,9], r'$(b_s,b_s)$':self.lpt_table[:,:,10],r'$(1,b_3)$':0.5*self.lpt_table[:,:,11],r'$(b_1,b_3)$': 0.5*self.lpt_table[:,:,12]} - pgg += ((b21 + b22)[:, None] * Pdmd2 + - (bs1 + bs2)[:, None] * Pdms2 + - (bL11 * b22 + bL12 * b21)[:, None] * Pd1d2 + - (bL11 * bs2 + bL12 * bs1)[:, None] * Pd1s2 + - (b21 * b22)[:, None] * Pd2d2 + - (b21 * bs2 + b22 * bs1)[:, None] * Pd2s2 + - (bs1 * bs2)[:, None] * Ps2s2) + pgg_interpolated = interp1d(cleft_k[0], pgg, kind='cubic', fill_value='extrapolate')(self.k) - return pgg + return pgg_interpolated def calculate(self, state, want_derived=True, **params_values_dict): try: self.cleft_obj except: - self.init_cleft() + self.init_cleft(k_filter=1.e-1) Pk_mm = self._get_Pk_mm() diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 76e29fea..556f0138 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -1,5 +1,6 @@ # pytest -k bias -v . +import pdb import pytest import numpy as np @@ -116,3 +117,10 @@ def test_LPT_bias_compute_grid(): nonlinear=False) Pk_gg = lhood.provider.get_Pk_gg_grid() + + from matplotlib import pyplot as plt + plt.ion() + pdb.set_trace() + + plt.loglog(k, Pk_mm_lin[0]) + plt.loglog(k, Pk_gg[0]) From 3814833536a750a8a84c527930256b202a95f43f Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 1 Mar 2022 11:52:50 +0000 Subject: [PATCH 14/29] added pgm --- soliket/bias.py | 37 +++++++++++++++++++++++++++++++------ soliket/tests/test_bias.py | 13 +++++++------ 2 files changed, 38 insertions(+), 12 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 8530dc4d..9d384aea 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -135,12 +135,12 @@ def init_cleft(self, k_filter=None): def _get_Pk_gg(self, Pk_mm, **params_values_dict): - b11 = params_values_dict['b11']# * np.ones_like(Pk_mm) - b21 = params_values_dict['b21']# * np.ones_like(Pk_mm) - bs1 = params_values_dict['bs1']# * np.ones_like(Pk_mm) - b12 = params_values_dict['b12']# * np.ones_like(Pk_mm) - b22 = params_values_dict['b22']# * np.ones_like(Pk_mm) - bs2 = params_values_dict['bs2']# * np.ones_like(Pk_mm) + b11 = params_values_dict['b1g1']# * np.ones_like(Pk_mm) + b21 = params_values_dict['b2g1']# * np.ones_like(Pk_mm) + bs1 = params_values_dict['bsg1']# * np.ones_like(Pk_mm) + b12 = params_values_dict['b1g2']# * np.ones_like(Pk_mm) + b22 = params_values_dict['b2g2']# * np.ones_like(Pk_mm) + bs2 = params_values_dict['bsg2']# * np.ones_like(Pk_mm) bL11 = b11 - 1 bL12 = b12 - 1 @@ -179,6 +179,31 @@ def _get_Pk_gg(self, Pk_mm, **params_values_dict): return pgg_interpolated + def get_pgm(self, Pk_mm, **params_values_dict): + + b1 = params_values_dict['b1g1']# * np.ones_like(Pk_mm) + b2 = params_values_dict['b2g1']# * np.ones_like(Pk_mm) + bs = params_values_dict['bsg1']# * np.ones_like(Pk_mm) + + bL1 = b1-1 + if Pnl is None: + Pdmdm = self.lpt_table[:, :, 1] + Pdmd1 = 0.5*self.lpt_table[:, :, 2] + pgm = Pdmdm + bL1[:, None] * Pdmd1 + else: + pgm = b1[:, None]*Pnl + Pdmd2 = 0.5*self.lpt_table[:, :, 4] + Pdms2 = 0.25*self.lpt_table[:, :, 7] + + pgm += (b2[:, None] * Pdmd2 + + bs[:, None] * Pdms2) + + cleft_k = self.lpt_table[:, :, 0] + + pgm_interpolated = interp1d(cleft_k[0], pgm, kind='cubic', fill_value='extrapolate')(self.k) + + return pgm_interpolated + def calculate(self, state, want_derived=True, **params_values_dict): try: diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 556f0138..1717f804 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -9,12 +9,12 @@ info = {"params": { "b_lin": 1.1, - "b11": 1., - "b21": 1., - "bs1": 1., - "b12": 1., - "b22": 1., - "bs2": 1., + "b1g1": 1., + "b2g1": 1., + "bsg1": 1., + "b1g2": 1., + "b2g2": 1., + "bsg2": 1., "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, @@ -117,6 +117,7 @@ def test_LPT_bias_compute_grid(): nonlinear=False) Pk_gg = lhood.provider.get_Pk_gg_grid() + Pk_gm = lhood.provider.get_Pk_gm_grid() from matplotlib import pyplot as plt plt.ion() From 42b75f567853fbdfb05fff9891b58d192152a413 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Fri, 11 Mar 2022 14:00:53 +0000 Subject: [PATCH 15/29] working (untested) pk_gm, sensible bias defaults --- soliket/bias.py | 30 +++++++++++++++++------------- soliket/tests/test_bias.py | 29 ++++++++++++++++------------- 2 files changed, 33 insertions(+), 26 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 9d384aea..75346942 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -105,12 +105,12 @@ def calculate(self, state, want_derived=True, **params_values_dict): class LPT_bias(Bias): - params = {'b11': None, - 'b21': None, - 'bs1': None, - 'b12': None, - 'b22': None, - 'bs2': None} + params = {'b1g1': None, + 'b2g1': None, + 'bsg1': None, + 'b1g2': None, + 'b2g2': None, + 'bsg2': None} def init_cleft(self, k_filter=None): @@ -179,24 +179,25 @@ def _get_Pk_gg(self, Pk_mm, **params_values_dict): return pgg_interpolated - def get_pgm(self, Pk_mm, **params_values_dict): + def _get_Pk_gm(self, Pk_mm, **params_values_dict): b1 = params_values_dict['b1g1']# * np.ones_like(Pk_mm) b2 = params_values_dict['b2g1']# * np.ones_like(Pk_mm) bs = params_values_dict['bsg1']# * np.ones_like(Pk_mm) bL1 = b1-1 - if Pnl is None: + if self.nonlinear: + pgm = b1 * Pk_mm + else: Pdmdm = self.lpt_table[:, :, 1] Pdmd1 = 0.5*self.lpt_table[:, :, 2] - pgm = Pdmdm + bL1[:, None] * Pdmd1 - else: - pgm = b1[:, None]*Pnl + pgm = Pdmdm + bL1 * Pdmd1 + Pdmd2 = 0.5*self.lpt_table[:, :, 4] Pdms2 = 0.25*self.lpt_table[:, :, 7] - pgm += (b2[:, None] * Pdmd2 + - bs[:, None] * Pdms2) + pgm += (b2 * Pdmd2 + + bs * Pdms2) cleft_k = self.lpt_table[:, :, 0] @@ -211,6 +212,9 @@ def calculate(self, state, want_derived=True, **params_values_dict): except: self.init_cleft(k_filter=1.e-1) + # do init at init, then update_power_spectrum in calculate + Pk_mm = self._get_Pk_mm() state['Pk_gg_grid'] = self._get_Pk_gg(Pk_mm, **params_values_dict) + state['Pk_gm_grid'] = self._get_Pk_gm(Pk_mm, **params_values_dict) diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 1717f804..5b22240e 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -8,13 +8,13 @@ from cobaya.run import run info = {"params": { - "b_lin": 1.1, - "b1g1": 1., - "b2g1": 1., - "bsg1": 1., - "b1g2": 1., - "b2g2": 1., - "bsg2": 1., + "b_lin": 1.0, + "b1g1": 1.0, + "b2g1": 0.0, + "bsg1": 0.0, + "b1g2": 1.0, + "b2g2": 0.0, + "bsg2": 0.0, "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, @@ -117,11 +117,14 @@ def test_LPT_bias_compute_grid(): nonlinear=False) Pk_gg = lhood.provider.get_Pk_gg_grid() - Pk_gm = lhood.provider.get_Pk_gm_grid() + # Pk_gm = lhood.provider.get_Pk_gm_grid() from matplotlib import pyplot as plt - plt.ion() - pdb.set_trace() - - plt.loglog(k, Pk_mm_lin[0]) - plt.loglog(k, Pk_gg[0]) + # plt.ion() + # pdb.set_trace() + + # plt.loglog(k, Pk_mm_lin[0]) + # plt.loglog(k, Pk_gg[0]) + # plt.xlim([1.e-3, 1.e0]) + # plt.ylim([1.e3, 1.e5]) + # plt.savefig('plots/lpt-bias-pk.png', dpi=300, bbox_inches='tight') From 0d3f9378c4e407d3fb3a8a467dae27aa17654fe8 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Fri, 11 Mar 2022 14:46:09 +0000 Subject: [PATCH 16/29] split cleft object to init and update --- soliket/bias.py | 49 +++++++++++++++++++++----------------- soliket/tests/test_bias.py | 2 +- 2 files changed, 28 insertions(+), 23 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 75346942..24046cb1 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -27,6 +27,7 @@ class Bias(Theory): def initialize(self): self._var_pairs = set() + self.k = None # self._var_pairs = [('delta_tot', 'delta_tot')] def get_requirements(self): @@ -59,28 +60,25 @@ def must_provide(self, **requirements): assert len(self._var_pairs) < 2, "Bias doesn't support other Pk yet" return needs - def _get_growth(self): - for pair in self._var_pairs: - - k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, - nonlinear=False) - - assert(z[0] == 0) - - return np.mean(Pk_mm / Pk_mm[:1], axis=-1)**0.5 - def _get_Pk_mm(self): + def _get_Pk_mm(self, update_growth=True): for pair in self._var_pairs: if self.nonlinear: - self.k, self.z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, - nonlinear=True) - # Pk_mm = np.flip(Pk_nonlin, axis=0) + k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, + nonlinear=True) else: - self.k, self.z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, - nonlinear=False) - # Pk_mm = np.flip(Pk_lin, axis=0) + k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=pair, + nonlinear=False) + if self.k is None: + self.k = k + if self.z is None: + self.z = z + + if update_growth: + assert(z[0] == 0) + self.Dz = np.mean(Pk_mm / Pk_mm[:1], axis=-1)**0.5 return Pk_mm @@ -97,7 +95,7 @@ class Linear_bias(Bias): def calculate(self, state, want_derived=True, **params_values_dict): - Pk_mm = self._get_Pk_mm() + Pk_mm = self._get_Pk_mm(update_growth=False) # growth not needed for linear bias state['Pk_gg_grid'] = params_values_dict['b_lin']**2. * Pk_mm state['Pk_gm_grid'] = params_values_dict['b_lin'] * Pk_mm @@ -114,25 +112,32 @@ class LPT_bias(Bias): def init_cleft(self, k_filter=None): - Pk_mm = self._get_Pk_mm() + # Pk_mm = self._get_Pk_mm() if k_filter is not None: self.wk_low = 1-np.exp(-(self.k/k_filter)**2) else: self.wk_low = np.ones_like(self.k) + + # self.update_lpt_table(Pk_mm) + + def update_lpt_table(self, Pk_mm): + self.cleft_obj = RKECLEFT(self.k, Pk_mm[0], extrap_min=np.floor(np.log10(self.k[0])), extrap_max=np.ceil(np.log10(self.k[-1]))) self.lpt_table = [] - for D in self._get_growth(): + for D in self.Dz: self.cleft_obj.make_ptable(D=D, kmin=self.k[0], kmax=self.k[-1], nk=self.k.size) self.lpt_table.append(self.cleft_obj.pktable) self.lpt_table = np.array(self.lpt_table) + + def _get_Pk_gg(self, Pk_mm, **params_values_dict): b11 = params_values_dict['b1g1']# * np.ones_like(Pk_mm) @@ -207,14 +212,14 @@ def _get_Pk_gm(self, Pk_mm, **params_values_dict): def calculate(self, state, want_derived=True, **params_values_dict): + Pk_mm = self._get_Pk_mm() + try: self.cleft_obj except: self.init_cleft(k_filter=1.e-1) - # do init at init, then update_power_spectrum in calculate - - Pk_mm = self._get_Pk_mm() + self.update_lpt_table(Pk_mm) state['Pk_gg_grid'] = self._get_Pk_gg(Pk_mm, **params_values_dict) state['Pk_gm_grid'] = self._get_Pk_gm(Pk_mm, **params_values_dict) diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 5b22240e..46e6c437 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -101,7 +101,7 @@ def test_LPT_bias_compute_grid(): } model = get_model(info) # noqa F841 - model.add_requirements({"Pk_grid": {"z": 0., "k_max": 10., + model.add_requirements({"Pk_grid": {"z": 0., "k_max": 1., "nonlinear": False, "vars_pairs": ('delta_tot', 'delta_tot'), "extrap_kmin": 1.e-3 From ad41b15a82981e2f2903880c45c71c412ac95ca8 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 26 Apr 2022 11:24:51 +0100 Subject: [PATCH 17/29] model and tests complete --- soliket/bias.py | 109 +++++++++++++++++++++---------------- soliket/tests/test_bias.py | 32 ++++------- 2 files changed, 72 insertions(+), 69 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 24046cb1..b0e56fcd 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -105,24 +105,36 @@ class LPT_bias(Bias): params = {'b1g1': None, 'b2g1': None, - 'bsg1': None, + # 'bsg1': None, 'b1g2': None, 'b2g2': None, - 'bsg2': None} + # 'bsg2': None + } def init_cleft(self, k_filter=None): - # Pk_mm = self._get_Pk_mm() + self.update_lpt_table() if k_filter is not None: - self.wk_low = 1-np.exp(-(self.k/k_filter)**2) + self.wk_low = 1 - np.exp(-(self.k / k_filter)**2) else: self.wk_low = np.ones_like(self.k) - # self.update_lpt_table(Pk_mm) - - def update_lpt_table(self, Pk_mm): + def update_lpt_table(self): + + k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + nonlinear=False) + + if self.k is None: + self.k = k + else: + assert np.allclose(k, self.k) # check we are consistent with ks + + if self.z is None: + self.z = z + else: + assert np.allclose(z, self.z) # check we are consistent with zs self.cleft_obj = RKECLEFT(self.k, Pk_mm[0], extrap_min=np.floor(np.log10(self.k[0])), @@ -130,33 +142,37 @@ def update_lpt_table(self, Pk_mm): self.lpt_table = [] + assert(z[0] == 0) + self.Dz = np.mean(Pk_mm / Pk_mm[:1], axis=-1)**0.5 + for D in self.Dz: - self.cleft_obj.make_ptable(D=D, kmin=self.k[0], kmax=self.k[-1], nk=self.k.size) + self.cleft_obj.make_ptable(D=D, + kmin=self.k[0], + kmax=self.k[-1], + nk=self.k.size) self.lpt_table.append(self.cleft_obj.pktable) self.lpt_table = np.array(self.lpt_table) + def _get_Pk_gg(self, **params_values_dict): - def _get_Pk_gg(self, Pk_mm, **params_values_dict): - - b11 = params_values_dict['b1g1']# * np.ones_like(Pk_mm) - b21 = params_values_dict['b2g1']# * np.ones_like(Pk_mm) - bs1 = params_values_dict['bsg1']# * np.ones_like(Pk_mm) - b12 = params_values_dict['b1g2']# * np.ones_like(Pk_mm) - b22 = params_values_dict['b2g2']# * np.ones_like(Pk_mm) - bs2 = params_values_dict['bsg2']# * np.ones_like(Pk_mm) + b11 = params_values_dict['b1g1'] + b21 = params_values_dict['b2g1'] + # bs1 = params_values_dict['bsg1'] + bs1 = 0.0 # ignore bs terms for now + b12 = params_values_dict['b1g2'] + b22 = params_values_dict['b2g2'] + # bs2 = params_values_dict['bsg2'] + bs2 = 0.0 # ignore bs terms for now bL11 = b11 - 1 bL12 = b12 - 1 - if self.nonlinear: - pgg = (b11 * b12) * Pk_mm - else: - Pdmdm = self.lpt_table[:, :, 1] - Pdmd1 = 0.5 * self.lpt_table[:, :, 2] - Pd1d1 = self.lpt_table[:, :, 3] - pgg = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) + Pdmdm = self.lpt_table[:, :, 1] + Pdmd1 = 0.5 * self.lpt_table[:, :, 2] + Pd1d1 = self.lpt_table[:, :, 3] + pgg = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) Pdmd2 = 0.5 * self.lpt_table[:, :, 4] Pd1d2 = 0.5 * self.lpt_table[:, :, 5] @@ -176,50 +192,47 @@ def _get_Pk_gg(self, Pk_mm, **params_values_dict): cleft_k = self.lpt_table[:, :, 0] - # pdb.set_trace() - - # spectra = {r'$(1,1)$':self.lpt_table[:,:,1],r'$(1,b_1)$':0.5*self.lpt_table[:,:,2], r'$(b_1,b_1)$': self.lpt_table[:,:,3],r'$(1,b_2)$':0.5*self.lpt_table[:,:,4], r'$(b_1,b_2)$': 0.5*self.lpt_table[:,:,5], r'$(b_2,b_2)$': self.lpt_table[:,:,6],r'$(1,b_s)$':0.5*self.lpt_table[:,:,7], r'$(b_1,b_s)$': 0.5*self.lpt_table[:,:,8], r'$(b_2,b_s)$':0.5*self.lpt_table[:,:,9], r'$(b_s,b_s)$':self.lpt_table[:,:,10],r'$(1,b_3)$':0.5*self.lpt_table[:,:,11],r'$(b_1,b_3)$': 0.5*self.lpt_table[:,:,12]} - - pgg_interpolated = interp1d(cleft_k[0], pgg, kind='cubic', fill_value='extrapolate')(self.k) + pgg_interpolated = interp1d(cleft_k[0], pgg, + kind='cubic', + fill_value='extrapolate')(self.k) return pgg_interpolated - def _get_Pk_gm(self, Pk_mm, **params_values_dict): + def _get_Pk_gm(self, **params_values_dict): - b1 = params_values_dict['b1g1']# * np.ones_like(Pk_mm) - b2 = params_values_dict['b2g1']# * np.ones_like(Pk_mm) - bs = params_values_dict['bsg1']# * np.ones_like(Pk_mm) + b1 = params_values_dict['b1g1'] + b2 = params_values_dict['b2g1'] + # bs = params_values_dict['bsg1'] + bs = 0.0 # ignore bs terms for now - bL1 = b1-1 - if self.nonlinear: - pgm = b1 * Pk_mm - else: - Pdmdm = self.lpt_table[:, :, 1] - Pdmd1 = 0.5*self.lpt_table[:, :, 2] - pgm = Pdmdm + bL1 * Pdmd1 + bL1 = b1 - 1 - Pdmd2 = 0.5*self.lpt_table[:, :, 4] - Pdms2 = 0.25*self.lpt_table[:, :, 7] + Pdmdm = self.lpt_table[:, :, 1] + Pdmd1 = 0.5 * self.lpt_table[:, :, 2] + pgm = Pdmdm + bL1 * Pdmd1 + + Pdmd2 = 0.5 * self.lpt_table[:, :, 4] + Pdms2 = 0.25 * self.lpt_table[:, :, 7] pgm += (b2 * Pdmd2 + bs * Pdms2) cleft_k = self.lpt_table[:, :, 0] - pgm_interpolated = interp1d(cleft_k[0], pgm, kind='cubic', fill_value='extrapolate')(self.k) + pgm_interpolated = interp1d(cleft_k[0], pgm, + kind='cubic', + fill_value='extrapolate')(self.k) return pgm_interpolated def calculate(self, state, want_derived=True, **params_values_dict): - Pk_mm = self._get_Pk_mm() - try: self.cleft_obj except: - self.init_cleft(k_filter=1.e-1) + self.init_cleft(k_filter=1.e-2) - self.update_lpt_table(Pk_mm) + self.update_lpt_table() - state['Pk_gg_grid'] = self._get_Pk_gg(Pk_mm, **params_values_dict) - state['Pk_gm_grid'] = self._get_Pk_gm(Pk_mm, **params_values_dict) + state['Pk_gg_grid'] = self._get_Pk_gg(**params_values_dict) + state['Pk_gm_grid'] = self._get_Pk_gm(**params_values_dict) diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 46e6c437..c1872d51 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -10,11 +10,9 @@ info = {"params": { "b_lin": 1.0, "b1g1": 1.0, - "b2g1": 0.0, - "bsg1": 0.0, + "b2g1": 1.0, "b1g2": 1.0, - "b2g2": 0.0, - "bsg2": 0.0, + "b2g2": 1.0, "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, @@ -91,40 +89,32 @@ def test_LPT_bias_model(): model = get_model(info) # noqa F841 + def test_LPT_bias_compute_grid(): from soliket.bias import LPT_bias info["theory"] = { "camb": None, - "LPT_bias": {"external": LPT_bias} + "LPT_bias": {"external": LPT_bias, + "nonlinear": True} } model = get_model(info) # noqa F841 model.add_requirements({"Pk_grid": {"z": 0., "k_max": 1., - "nonlinear": False, + "nonlinear": True, "vars_pairs": ('delta_tot', 'delta_tot'), - "extrap_kmin": 1.e-3 }, - "Pk_gg_grid": None + "Pk_gg_grid": None, + "Pk_gm_grid": None }) model.logposterior(info['params']) # force computation of model lhood = model.likelihood['one'] - k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), - nonlinear=False) - Pk_gg = lhood.provider.get_Pk_gg_grid() - # Pk_gm = lhood.provider.get_Pk_gm_grid() - - from matplotlib import pyplot as plt - # plt.ion() - # pdb.set_trace() + Pk_gm = lhood.provider.get_Pk_gm_grid() - # plt.loglog(k, Pk_mm_lin[0]) - # plt.loglog(k, Pk_gg[0]) - # plt.xlim([1.e-3, 1.e0]) - # plt.ylim([1.e3, 1.e5]) - # plt.savefig('plots/lpt-bias-pk.png', dpi=300, bbox_inches='tight') + assert np.isclose(Pk_gg.sum(), 493326018.33041435) + assert np.isclose(Pk_gm.sum(), 429303138.83088464) From 93766032d5adb1c4e886c4ffced788a4c0bff024 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 26 Apr 2022 11:41:20 +0100 Subject: [PATCH 18/29] added test skip if velocileptors import fails --- soliket/bias.py | 5 ++++- soliket/tests/test_bias.py | 2 ++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/soliket/bias.py b/soliket/bias.py index b0e56fcd..731fa6de 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -9,7 +9,10 @@ from cobaya.theory import Theory from cobaya.log import LoggedError -from velocileptors.EPT.cleft_kexpanded_resummed_fftw import RKECLEFT +try: + from velocileptors.EPT.cleft_kexpanded_resummed_fftw import RKECLEFT +except: + "velocileptors import failed, LPT_bias will be unavailable" class Bias(Theory): diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index c1872d51..a96ace5d 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -80,6 +80,7 @@ def test_linear_bias_compute_grid(): def test_LPT_bias_model(): + skip_lpt = pytest.importorskip("velocileptors") # noqa F841 from soliket.bias import LPT_bias info["theory"] = { @@ -92,6 +93,7 @@ def test_LPT_bias_model(): def test_LPT_bias_compute_grid(): + skip_lpt = pytest.importorskip("velocileptors") # noqa F841 from soliket.bias import LPT_bias info["theory"] = { From a638e5c98c59684e80d1712a84e2b325e8f6d0d5 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 26 Apr 2022 16:14:28 +0100 Subject: [PATCH 19/29] HaloFit model and interpolations --- soliket/bias.py | 71 ++++++++++++++++++++++++++++++-------- soliket/tests/test_bias.py | 27 ++++++++++++--- 2 files changed, 80 insertions(+), 18 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 731fa6de..97d39b02 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -105,6 +105,29 @@ def calculate(self, state, want_derived=True, **params_values_dict): class LPT_bias(Bias): + '''Theory object for computation of galaxy-galaxy and galaxy-matter power spectra, + based on a Lagrangian Effective Field Theory model for the galaxy bias. + + The model consists of: + - Input bias parameters in the Lagrangian definition. + - Use of the velocileptors code [1]_ to calculate power spectrum expansion terms + - Only terms up to b_1 and b_2 (i.e. no shear or third order terms) + - Use of HaloFit non-linear matter power spectrum to account for small-scales, + replacing the use of EFT counter-terms. This follows the approach of + [2]_ (equation X) and [3]_ (table Y). + - The further substitutions P_{b_1} = 2 P^{HaloFit}_{mm} + and P_{b^2_1} = P^{HaloFit}_{mm} + + Or, in mathemetics: + P_{gm}(k) = P^{HaloFit}_{mm} + b_1 P^{HaloFit}_{mm} + \frac{b_2}{2}P_{b_{2}} + P_{mm}(k) = P^{HaloFit}_{mm} + + 2 b_1 P^{HaloFit}_{mm} + b^{2}_1 P^{HaloFit}_{mm} + + b_1 b_2 P_{b_{1}b_{2}} + + b_2 P_{b_2} + b^2_2 P_{b^2_2} + where e.g. P_{b_1} is the \langle 1, \delta \rangle power spectrum and P_{b_2} is + the \langle 1, \delta^2 \rangle power spectrum. + + ''' params = {'b1g1': None, 'b2g1': None, @@ -172,10 +195,14 @@ def _get_Pk_gg(self, **params_values_dict): bL11 = b11 - 1 bL12 = b12 - 1 - Pdmdm = self.lpt_table[:, :, 1] - Pdmd1 = 0.5 * self.lpt_table[:, :, 2] - Pd1d1 = self.lpt_table[:, :, 3] - pgg = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) + if self.nonlinear: + Pk_mm = self._get_Pk_mm(update_growth=False) + pgg = Pk_mm + (bL11 + bL12) * Pk_mm + (bL11 * bL12) * Pk_mm + else: + Pdmdm = self.lpt_table[:, :, 1] + Pdmd1 = 0.5 * self.lpt_table[:, :, 2] + Pd1d1 = self.lpt_table[:, :, 3] + pgg = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) Pdmd2 = 0.5 * self.lpt_table[:, :, 4] Pd1d2 = 0.5 * self.lpt_table[:, :, 5] @@ -210,21 +237,37 @@ def _get_Pk_gm(self, **params_values_dict): bL1 = b1 - 1 - Pdmdm = self.lpt_table[:, :, 1] - Pdmd1 = 0.5 * self.lpt_table[:, :, 2] - pgm = Pdmdm + bL1 * Pdmd1 - + cleft_k = self.lpt_table[:, :, 0] Pdmd2 = 0.5 * self.lpt_table[:, :, 4] Pdms2 = 0.25 * self.lpt_table[:, :, 7] - pgm += (b2 * Pdmd2 + - bs * Pdms2) + pgm_higher_order = (b2 * Pdmd2 + + bs * Pdms2) - cleft_k = self.lpt_table[:, :, 0] + if self.nonlinear: + # here the nonlinear matter-matter power spectrum from the upstream Theory + # (i.e. HaloFit) is used for leading order terms, so only higher order terms + # come from the lpt table and hence need interpolation. + pgm_higher_order_interp = interp1d(cleft_k[0], pgm_higher_order, + kind='cubic', + fill_value='extrapolate')(self.k) - pgm_interpolated = interp1d(cleft_k[0], pgm, - kind='cubic', - fill_value='extrapolate')(self.k) + Pk_mm = self._get_Pk_mm(update_growth=False) + pgm_lin = Pk_mm + bL1 * Pk_mm + + pgm_interpolated = pgm_lin + pgm_higher_order_interp + + else: + # here the lpt table is used for all order terms, so all need interpolation. + Pdmdm = self.lpt_table[:, :, 1] + Pdmd1 = 0.5 * self.lpt_table[:, :, 2] + pgm_lin = Pdmdm + bL1 * Pdmd1 + + pgm = pgm_lin + pgm_higher_order + + pgm_interpolated = interp1d(cleft_k[0], pgm, + kind='cubic', + fill_value='extrapolate')(self.k) return pgm_interpolated diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index a96ace5d..c064f341 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -9,10 +9,10 @@ info = {"params": { "b_lin": 1.0, - "b1g1": 1.0, + "b1g1": 1.1, "b2g1": 1.0, "b1g2": 1.0, - "b2g2": 1.0, + "b2g2": 0.0, "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, @@ -118,5 +118,24 @@ def test_LPT_bias_compute_grid(): Pk_gg = lhood.provider.get_Pk_gg_grid() Pk_gm = lhood.provider.get_Pk_gm_grid() - assert np.isclose(Pk_gg.sum(), 493326018.33041435) - assert np.isclose(Pk_gm.sum(), 429303138.83088464) + k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + nonlinear=False) + k_nl, z_nl, Pk_mm_nl = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + nonlinear=True) + + from matplotlib import pyplot as plt + plt.ion() + plt.loglog(k, Pk_mm_lin[0], label='$P_{lin}(k)$') + plt.loglog(k_nl, Pk_mm_nl[0], '--', label='$P_{HaloFit}(k)$') + plt.loglog(k_nl, Pk_gm[0], label='$P_{gm}(k)$') + plt.loglog(k_nl, Pk_gg[0], label='$P_{gg}(k)$') + plt.xlabel('$k\,$[Mpc$^{-1}$]') + plt.ylabel('$P(k)$') + plt.title('velocileptors power spectra') + plt.legend() + plt.xlim([1.e-3, 1.e0]) + plt.ylim([1.e3, 1.e5]) + plt.savefig('plots/lpt-bias-pk.png', dpi=300, bbox_inches='tight') + + # assert np.isclose(Pk_gg.sum(), 493326018.33041435) + # assert np.isclose(Pk_gm.sum(), 429303138.83088464) From b5ffe3e19ce2c82f2a3da07eb2a05f5acaef7163 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 26 Apr 2022 16:19:14 +0100 Subject: [PATCH 20/29] HaloFit model and interpolations --- soliket/bias.py | 64 +++++++++++++++++++++++++++++-------------------- 1 file changed, 38 insertions(+), 26 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 97d39b02..2bfb5673 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -195,15 +195,6 @@ def _get_Pk_gg(self, **params_values_dict): bL11 = b11 - 1 bL12 = b12 - 1 - if self.nonlinear: - Pk_mm = self._get_Pk_mm(update_growth=False) - pgg = Pk_mm + (bL11 + bL12) * Pk_mm + (bL11 * bL12) * Pk_mm - else: - Pdmdm = self.lpt_table[:, :, 1] - Pdmd1 = 0.5 * self.lpt_table[:, :, 2] - Pd1d1 = self.lpt_table[:, :, 3] - pgg = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) - Pdmd2 = 0.5 * self.lpt_table[:, :, 4] Pd1d2 = 0.5 * self.lpt_table[:, :, 5] Pd2d2 = self.lpt_table[:, :, 6] * self.wk_low[None, :] @@ -212,19 +203,40 @@ def _get_Pk_gg(self, **params_values_dict): Pd2s2 = 0.25 * self.lpt_table[:, :, 9] * self.wk_low[None, :] Ps2s2 = 0.25 * self.lpt_table[:, :, 10] * self.wk_low[None, :] - pgg += ((b21 + b22) * Pdmd2 + - (bs1 + bs2) * Pdms2 + - (bL11 * b22 + bL12 * b21) * Pd1d2 + - (bL11 * bs2 + bL12 * bs1) * Pd1s2 + - (b21 * b22) * Pd2d2 + - (b21 * bs2 + b22 * bs1) * Pd2s2 + - (bs1 * bs2) * Ps2s2) + pgg_higher_order = ((b21 + b22) * Pdmd2 + + (bs1 + bs2) * Pdms2 + + (bL11 * b22 + bL12 * b21) * Pd1d2 + + (bL11 * bs2 + bL12 * bs1) * Pd1s2 + + (b21 * b22) * Pd2d2 + + (b21 * bs2 + b22 * bs1) * Pd2s2 + + (bs1 * bs2) * Ps2s2) cleft_k = self.lpt_table[:, :, 0] - pgg_interpolated = interp1d(cleft_k[0], pgg, - kind='cubic', - fill_value='extrapolate')(self.k) + if self.nonlinear: + # Nonlinear matter-matter power spectrum from the upstream Theory + # (i.e. HaloFit) is used for leading order terms, so only higher order + # terms come from the lpt table and hence need interpolation. + pgg_higher_order_interp = interp1d(cleft_k[0], pgg_higher_order, + kind='cubic', + fill_value='extrapolate')(self.k) + + Pk_mm = self._get_Pk_mm(update_growth=False) + pgg_leading_order = Pk_mm + (bL11 + bL12) * Pk_mm + (bL11 * bL12) * Pk_mm + + pgg_interpolated = pgg_leading_order + pgg_higher_order_interp + + else: + Pdmdm = self.lpt_table[:, :, 1] + Pdmd1 = 0.5 * self.lpt_table[:, :, 2] + Pd1d1 = self.lpt_table[:, :, 3] + pgg_leading_order = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) + + pgg = pgg_leading_order + pgm_higher_order + + pgg_interpolated = interp1d(cleft_k[0], pgg, + kind='cubic', + fill_value='extrapolate')(self.k) return pgg_interpolated @@ -245,25 +257,25 @@ def _get_Pk_gm(self, **params_values_dict): bs * Pdms2) if self.nonlinear: - # here the nonlinear matter-matter power spectrum from the upstream Theory - # (i.e. HaloFit) is used for leading order terms, so only higher order terms - # come from the lpt table and hence need interpolation. + # Nonlinear matter-matter power spectrum from the upstream Theory + # (i.e. HaloFit) is used for leading order terms, so only higher order + # terms come from the lpt table and hence need interpolation. pgm_higher_order_interp = interp1d(cleft_k[0], pgm_higher_order, kind='cubic', fill_value='extrapolate')(self.k) Pk_mm = self._get_Pk_mm(update_growth=False) - pgm_lin = Pk_mm + bL1 * Pk_mm + pgm_leading_order = Pk_mm + bL1 * Pk_mm - pgm_interpolated = pgm_lin + pgm_higher_order_interp + pgm_interpolated = pgm_leading_order + pgm_higher_order_interp else: # here the lpt table is used for all order terms, so all need interpolation. Pdmdm = self.lpt_table[:, :, 1] Pdmd1 = 0.5 * self.lpt_table[:, :, 2] - pgm_lin = Pdmdm + bL1 * Pdmd1 + pgm_leading_order = Pdmdm + bL1 * Pdmd1 - pgm = pgm_lin + pgm_higher_order + pgm = pgm_leading_order + pgm_higher_order pgm_interpolated = interp1d(cleft_k[0], pgm, kind='cubic', From 24ff94b2d7b4180d4a4606fb5ae6fcb71330bce7 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 26 Apr 2022 16:44:09 +0100 Subject: [PATCH 21/29] halofit or vlpt options for leading order terms --- soliket/bias.py | 2 +- soliket/tests/test_bias.py | 59 +++++++++++++++++++++----------------- 2 files changed, 34 insertions(+), 27 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 2bfb5673..f8ea998f 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -232,7 +232,7 @@ def _get_Pk_gg(self, **params_values_dict): Pd1d1 = self.lpt_table[:, :, 3] pgg_leading_order = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) - pgg = pgg_leading_order + pgm_higher_order + pgg = pgg_leading_order + pgg_higher_order pgg_interpolated = interp1d(cleft_k[0], pgg, kind='cubic', diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index c064f341..79eb4c86 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -9,10 +9,10 @@ info = {"params": { "b_lin": 1.0, - "b1g1": 1.1, + "b1g1": 1.0, "b2g1": 1.0, "b1g2": 1.0, - "b2g2": 0.0, + "b2g2": 1.0, "H0": 70., "ombh2": 0.0245, "omch2": 0.1225, @@ -90,8 +90,8 @@ def test_LPT_bias_model(): model = get_model(info) # noqa F841 - -def test_LPT_bias_compute_grid(): +@pytest.mark.parametrize("nonlinear_model", [True, False]) +def test_LPT_bias_compute_grid(nonlinear_model): skip_lpt = pytest.importorskip("velocileptors") # noqa F841 from soliket.bias import LPT_bias @@ -99,7 +99,7 @@ def test_LPT_bias_compute_grid(): info["theory"] = { "camb": None, "LPT_bias": {"external": LPT_bias, - "nonlinear": True} + "nonlinear": nonlinear_model} } model = get_model(info) # noqa F841 @@ -118,24 +118,31 @@ def test_LPT_bias_compute_grid(): Pk_gg = lhood.provider.get_Pk_gg_grid() Pk_gm = lhood.provider.get_Pk_gm_grid() - k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), - nonlinear=False) - k_nl, z_nl, Pk_mm_nl = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), - nonlinear=True) - - from matplotlib import pyplot as plt - plt.ion() - plt.loglog(k, Pk_mm_lin[0], label='$P_{lin}(k)$') - plt.loglog(k_nl, Pk_mm_nl[0], '--', label='$P_{HaloFit}(k)$') - plt.loglog(k_nl, Pk_gm[0], label='$P_{gm}(k)$') - plt.loglog(k_nl, Pk_gg[0], label='$P_{gg}(k)$') - plt.xlabel('$k\,$[Mpc$^{-1}$]') - plt.ylabel('$P(k)$') - plt.title('velocileptors power spectra') - plt.legend() - plt.xlim([1.e-3, 1.e0]) - plt.ylim([1.e3, 1.e5]) - plt.savefig('plots/lpt-bias-pk.png', dpi=300, bbox_inches='tight') - - # assert np.isclose(Pk_gg.sum(), 493326018.33041435) - # assert np.isclose(Pk_gm.sum(), 429303138.83088464) + # k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + # nonlinear=False) + # k_nl, z_nl, Pk_mm_nl = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + # nonlinear=True) + + # from matplotlib import pyplot as plt + # if nonlinear_model: + # plot_tag = 'HaloFit' + # plt.loglog(k, Pk_mm_lin[0], label='$P_{mm, lin}$') + # plt.loglog(k_nl, Pk_mm_nl[0], label='$P_{mm, HaloFit}$') + # else: + # plot_tag = 'velocileptors' + # plt.loglog(k_nl, Pk_gm[0], '-.', label='$P_{gm}$'+' ({} leading order)'.format(plot_tag)) + # plt.loglog(k_nl, Pk_gg[0], '--', label='$P_{gg}$'+' ({} leading order)'.format(plot_tag)) + # plt.xlabel('$k\\,$[Mpc$^{-1}$]') + # plt.ylabel('$P(k)$') + # plt.legend() + # plt.xlim([1.e-3, 1.e0]) + # plt.ylim([1.e3, 1.e5]) + # plt.savefig('plots/lpt-bias-pk.png', dpi=300, bbox_inches='tight') + + if nonlinear_model: + assert np.isclose(Pk_gg.sum(), 489563565.9282355) + assert np.isclose(Pk_gm.sum(), 425540686.4287062) + else: + assert np.isclose(Pk_gg.sum(), 493325841.6713596) + assert np.isclose(Pk_gm.sum(), 429302962.1718302) + From 3408ec04621837bdd0b813861e4cbcba3bfcc8df Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 26 Apr 2022 17:07:05 +0100 Subject: [PATCH 22/29] improved docs --- soliket/bias.py | 58 +++++++++++++++++++++++++++++--------- soliket/tests/test_bias.py | 13 +++++---- 2 files changed, 52 insertions(+), 19 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index f8ea998f..9c8fa34a 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -111,21 +111,51 @@ class LPT_bias(Bias): The model consists of: - Input bias parameters in the Lagrangian definition. - Use of the velocileptors code [1]_ to calculate power spectrum expansion terms - - Only terms up to b_1 and b_2 (i.e. no shear or third order terms) - - Use of HaloFit non-linear matter power spectrum to account for small-scales, - replacing the use of EFT counter-terms. This follows the approach of - [2]_ (equation X) and [3]_ (table Y). - - The further substitutions P_{b_1} = 2 P^{HaloFit}_{mm} - and P_{b^2_1} = P^{HaloFit}_{mm} - + - Only terms up to :math:`b_1` and :math:`b_2` (i.e. no shear or + third order terms) + + If nonlinear = True we make use of the following, which is the approach of + [2]_ (equation 3.1-3.7) and [3]_ (Model C of Table 1).: + - HaloFit non-linear matter power spectrum to account for small-scales, + replacing the use of EFT counter-terms in the leading order terms. + - The further substitutions :math:`P_{b_1} = 2 P^{HaloFit}_{mm}` + and :math:`P_{b^2_1} = P^{HaloFit}_{mm}`. Or, in mathemetics: + + .. math:: P_{gm}(k) = P^{HaloFit}_{mm} + b_1 P^{HaloFit}_{mm} + \frac{b_2}{2}P_{b_{2}} - P_{mm}(k) = P^{HaloFit}_{mm} + - 2 b_1 P^{HaloFit}_{mm} + b^{2}_1 P^{HaloFit}_{mm} + - b_1 b_2 P_{b_{1}b_{2}} + - b_2 P_{b_2} + b^2_2 P_{b^2_2} - where e.g. P_{b_1} is the \langle 1, \delta \rangle power spectrum and P_{b_2} is - the \langle 1, \delta^2 \rangle power spectrum. + + P_{mm}(k) =& P^{HaloFit}_{mm} + \\ + & 2 b_1 P^{HaloFit}_{mm} + b^{2}_1 P^{HaloFit}_{mm} + \\ + & b_1 b_2 P_{b_{1}b_{2}} + \\ + & b_2 P_{b_2} + b^2_2 P_{b^2_2} + + where e.g. :math:`P_{b_1}` is the :math:`\\langle 1, \\delta \\rangle` power spectrum + and :math:`P_{b_2}` is the :math:`\\langle 1, \\delta^2 \\rangle` power spectrum. + + Parameters + ---------- + + nonlinear : bool + If True, use the non-linear matter power spectrum as calculated by an upstream + Theory when calculating the leading order terms (as detailed above). Otherwise + use velocileptors LPT terms for leading order. + + b1g1 : float + Leading order Lagrangian bias value for galaxy sample 1. + b1g2 : float + Leading order Lagrangian bias value for galaxy sample 2. + b2g1 : float + Second order Lagrangian bias value for galaxy sample 1. + b2g2 : float + Second order Lagrangian bias value for galaxy sample 2. + + + References + ---------- + .. [1] https://github.com/sfschen/velocileptors + .. [2] Krolewski, Ferraro and White, 2021, arXiv:2105.03421 + .. [3] Pandey et al 2020, arXiv:2008.05991 ''' @@ -232,7 +262,7 @@ def _get_Pk_gg(self, **params_values_dict): Pd1d1 = self.lpt_table[:, :, 3] pgg_leading_order = (Pdmdm + (bL11 + bL12) * Pdmd1 + (bL11 * bL12) * Pd1d1) - pgg = pgg_leading_order + pgg_higher_order + pgg = pgg_leading_order + pgg_higher_order pgg_interpolated = interp1d(cleft_k[0], pgg, kind='cubic', diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 79eb4c86..6b165648 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -90,6 +90,7 @@ def test_LPT_bias_model(): model = get_model(info) # noqa F841 + @pytest.mark.parametrize("nonlinear_model", [True, False]) def test_LPT_bias_compute_grid(nonlinear_model): @@ -118,9 +119,10 @@ def test_LPT_bias_compute_grid(nonlinear_model): Pk_gg = lhood.provider.get_Pk_gg_grid() Pk_gm = lhood.provider.get_Pk_gm_grid() - # k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + # k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot','delta_tot'), # nonlinear=False) - # k_nl, z_nl, Pk_mm_nl = lhood.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), + # k_nl, z_nl, Pk_mm_nl = lhood.provider.get_Pk_grid(var_pair=('delta_tot', + # 'delta_tot'), # nonlinear=True) # from matplotlib import pyplot as plt @@ -130,8 +132,10 @@ def test_LPT_bias_compute_grid(nonlinear_model): # plt.loglog(k_nl, Pk_mm_nl[0], label='$P_{mm, HaloFit}$') # else: # plot_tag = 'velocileptors' - # plt.loglog(k_nl, Pk_gm[0], '-.', label='$P_{gm}$'+' ({} leading order)'.format(plot_tag)) - # plt.loglog(k_nl, Pk_gg[0], '--', label='$P_{gg}$'+' ({} leading order)'.format(plot_tag)) + # plt.loglog(k_nl, Pk_gm[0], '-.', + # label='$P_{gm}$'+' ({} leading order)'.format(plot_tag)) + # plt.loglog(k_nl, Pk_gg[0], '--', + # label='$P_{gg}$'+' ({} leading order)'.format(plot_tag)) # plt.xlabel('$k\\,$[Mpc$^{-1}$]') # plt.ylabel('$P(k)$') # plt.legend() @@ -145,4 +149,3 @@ def test_LPT_bias_compute_grid(nonlinear_model): else: assert np.isclose(Pk_gg.sum(), 493325841.6713596) assert np.isclose(Pk_gm.sum(), 429302962.1718302) - From 294691f2f04009e8771b11cb457b5099f4c92fc8 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 10 May 2022 14:00:35 +0100 Subject: [PATCH 23/29] added interpolator interface --- soliket/bias.py | 67 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/soliket/bias.py b/soliket/bias.py index 9c8fa34a..558079e9 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -9,6 +9,7 @@ from cobaya.theory import Theory from cobaya.log import LoggedError +from cobaya.theories.cosmo.boltzmannbase import PowerSpectrumInterpolator try: from velocileptors.EPT.cleft_kexpanded_resummed_fftw import RKECLEFT except: @@ -63,6 +64,66 @@ def must_provide(self, **requirements): assert len(self._var_pairs) < 2, "Bias doesn't support other Pk yet" return needs + def get_Pk_mm_interpolator(self, var_pair=("delta_tot", "delta_tot"), + nonlinear=True, extrap_kmax=None): + + return self._get_Pk_interpolator(pk_type='mm', + var_pair=var_pair, + nonlinear=nonlinear, + extrap_kmax=extrap_kmax) + + def get_Pk_gg_interpolator(self, var_pair=("delta_tot", "delta_tot"), + nonlinear=True, extrap_kmax=None): + + return self._get_Pk_interpolator(pk_type='gg', + var_pair=var_pair, + nonlinear=nonlinear, + extrap_kmax=extrap_kmax) + + def get_Pk_gm_interpolator(self, var_pair=("delta_tot", "delta_tot"), + nonlinear=True, extrap_kmax=None): + + return self._get_Pk_interpolator(pk_type='gm', + var_pair=var_pair, + nonlinear=nonlinear, + extrap_kmax=extrap_kmax) + + def _get_Pk_interpolator(self, pk_type, + var_pair=("delta_tot", "delta_tot"), nonlinear=True, + extrap_kmax=None): + + nonlinear = bool(nonlinear) + + key = ("Pk_{}_interpolator".format(pk_type), nonlinear, extrap_kmax)\ + + tuple(sorted(var_pair)) + if key in self.current_state: + return self.current_state[key] + + if pk_type == 'mm': + k, z, pk = self.get_Pk_mm(var_pair=var_pair, nonlinear=nonlinear) + elif pk_type == 'gm': + k, z, pk = self.get_Pk_gm_grid(var_pair=var_pair, nonlinear=nonlinear) + elif pk_type == 'gg': + k, z, pk = self.get_Pk_gg_grid(var_pair=var_pair, nonlinear=nonlinear) + + log_p = True + sign = 1 + if np.any(pk < 0): + if np.all(pk < 0): + sign = -1 + else: + log_p = False + if log_p: + pk = np.log(sign * pk) + elif extrap_kmax > k[-1]: + raise LoggedError(self.log, + 'Cannot do log extrapolation with zero-crossing pk ' + 'for %s, %s' % var_pair) + result = PowerSpectrumInterpolator(z, k, pk, logP=log_p, logsign=sign, + extrap_kmax=extrap_kmax) + self.current_state[key] = result + return result + def _get_Pk_mm(self, update_growth=True): @@ -76,8 +137,12 @@ def _get_Pk_mm(self, update_growth=True): nonlinear=False) if self.k is None: self.k = k + else: + assert np.allclose(k, self.k) # check we are consistent with ks if self.z is None: self.z = z + else: + assert np.allclose(z, self.z) # check we are consistent with zs if update_growth: assert(z[0] == 0) @@ -180,7 +245,7 @@ def init_cleft(self, k_filter=None): def update_lpt_table(self): k, z, Pk_mm = self.provider.get_Pk_grid(var_pair=('delta_tot', 'delta_tot'), - nonlinear=False) + nonlinear=False) # needs to be linear if self.k is None: self.k = k From b7b53c8df2d656ca45d329b6c4d7e1a80010fefb Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 10 May 2022 14:22:49 +0100 Subject: [PATCH 24/29] include velocileptors as requirement --- setup.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 9075ccc5..ddc7af4c 100644 --- a/setup.py +++ b/setup.py @@ -25,8 +25,10 @@ "cobaya", "sacc", "pyccl", + "pyfftw", "fgspectra @ git+https://github.com/simonsobs/fgspectra@act_sz_x_cib#egg=fgspectra", # noqa E501 - "mflike @ git+https://github.com/simonsobs/lat_mflike@master" + "mflike @ git+https://github.com/simonsobs/lat_mflike@master", + "velocileptors @ https://github.com/sfschen/velocileptors" ], test_suite="soliket.tests", include_package_data=True, From 68c3c6f349984505f698c165189c8cf961dfaced Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 10 May 2022 14:36:07 +0100 Subject: [PATCH 25/29] added velocileptors branch --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index ddc7af4c..198ee01d 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ "pyfftw", "fgspectra @ git+https://github.com/simonsobs/fgspectra@act_sz_x_cib#egg=fgspectra", # noqa E501 "mflike @ git+https://github.com/simonsobs/lat_mflike@master", - "velocileptors @ https://github.com/sfschen/velocileptors" + "velocileptors @ https://github.com/sfschen/velocileptors@master" ], test_suite="soliket.tests", include_package_data=True, From ca403de6c9e9755e1a6c8c2e714cfb1cbe5f5f7c Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 10 May 2022 14:47:31 +0100 Subject: [PATCH 26/29] try again with velocileptors requirement... --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 198ee01d..d1efcdc6 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ "pyfftw", "fgspectra @ git+https://github.com/simonsobs/fgspectra@act_sz_x_cib#egg=fgspectra", # noqa E501 "mflike @ git+https://github.com/simonsobs/lat_mflike@master", - "velocileptors @ https://github.com/sfschen/velocileptors@master" + "velocileptors @ git+https://github.com/sfschen/velocileptors@master" ], test_suite="soliket.tests", include_package_data=True, From 1ae811ed43b1bf309d050c07f17c3fc205c5ee21 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 10 May 2022 15:17:27 +0100 Subject: [PATCH 27/29] expanded tests to include interpolator --- soliket/bias.py | 31 +++++++++---------- soliket/tests/test_bias.py | 63 +++++++++++++++++++++++--------------- 2 files changed, 53 insertions(+), 41 deletions(-) diff --git a/soliket/bias.py b/soliket/bias.py index 558079e9..7b46ffdb 100644 --- a/soliket/bias.py +++ b/soliket/bias.py @@ -64,28 +64,25 @@ def must_provide(self, **requirements): assert len(self._var_pairs) < 2, "Bias doesn't support other Pk yet" return needs - def get_Pk_mm_interpolator(self, var_pair=("delta_tot", "delta_tot"), - nonlinear=True, extrap_kmax=None): + def get_Pk_mm_interpolator(self, extrap_kmax=None): return self._get_Pk_interpolator(pk_type='mm', - var_pair=var_pair, - nonlinear=nonlinear, + var_pair=self._var_pairs, + nonlinear=self.nonlinear, extrap_kmax=extrap_kmax) - def get_Pk_gg_interpolator(self, var_pair=("delta_tot", "delta_tot"), - nonlinear=True, extrap_kmax=None): + def get_Pk_gg_interpolator(self, extrap_kmax=None): return self._get_Pk_interpolator(pk_type='gg', - var_pair=var_pair, - nonlinear=nonlinear, + var_pair=self._var_pairs, + nonlinear=self.nonlinear, extrap_kmax=extrap_kmax) - def get_Pk_gm_interpolator(self, var_pair=("delta_tot", "delta_tot"), - nonlinear=True, extrap_kmax=None): + def get_Pk_gm_interpolator(self, extrap_kmax=None): return self._get_Pk_interpolator(pk_type='gm', - var_pair=var_pair, - nonlinear=nonlinear, + var_pair=self._var_pairs, + nonlinear=self.nonlinear, extrap_kmax=extrap_kmax) def _get_Pk_interpolator(self, pk_type, @@ -100,11 +97,11 @@ def _get_Pk_interpolator(self, pk_type, return self.current_state[key] if pk_type == 'mm': - k, z, pk = self.get_Pk_mm(var_pair=var_pair, nonlinear=nonlinear) + pk = self._get_Pk_mm(update_growth=False) elif pk_type == 'gm': - k, z, pk = self.get_Pk_gm_grid(var_pair=var_pair, nonlinear=nonlinear) + pk = self.get_Pk_gm_grid() elif pk_type == 'gg': - k, z, pk = self.get_Pk_gg_grid(var_pair=var_pair, nonlinear=nonlinear) + pk = self.get_Pk_gg_grid() log_p = True sign = 1 @@ -115,11 +112,11 @@ def _get_Pk_interpolator(self, pk_type, log_p = False if log_p: pk = np.log(sign * pk) - elif extrap_kmax > k[-1]: + elif (extrap_kmax is not None) and (extrap_kmax > self.k[-1]): raise LoggedError(self.log, 'Cannot do log extrapolation with zero-crossing pk ' 'for %s, %s' % var_pair) - result = PowerSpectrumInterpolator(z, k, pk, logP=log_p, logsign=sign, + result = PowerSpectrumInterpolator(self.z, self.k, pk, logP=log_p, logsign=sign, extrap_kmax=extrap_kmax) self.current_state[key] = result return result diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 6b165648..0708bf98 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -119,33 +119,48 @@ def test_LPT_bias_compute_grid(nonlinear_model): Pk_gg = lhood.provider.get_Pk_gg_grid() Pk_gm = lhood.provider.get_Pk_gm_grid() - # k, z, Pk_mm_lin = lhood.provider.get_Pk_grid(var_pair=('delta_tot','delta_tot'), - # nonlinear=False) - # k_nl, z_nl, Pk_mm_nl = lhood.provider.get_Pk_grid(var_pair=('delta_tot', - # 'delta_tot'), - # nonlinear=True) - - # from matplotlib import pyplot as plt - # if nonlinear_model: - # plot_tag = 'HaloFit' - # plt.loglog(k, Pk_mm_lin[0], label='$P_{mm, lin}$') - # plt.loglog(k_nl, Pk_mm_nl[0], label='$P_{mm, HaloFit}$') - # else: - # plot_tag = 'velocileptors' - # plt.loglog(k_nl, Pk_gm[0], '-.', - # label='$P_{gm}$'+' ({} leading order)'.format(plot_tag)) - # plt.loglog(k_nl, Pk_gg[0], '--', - # label='$P_{gg}$'+' ({} leading order)'.format(plot_tag)) - # plt.xlabel('$k\\,$[Mpc$^{-1}$]') - # plt.ylabel('$P(k)$') - # plt.legend() - # plt.xlim([1.e-3, 1.e0]) - # plt.ylim([1.e3, 1.e5]) - # plt.savefig('plots/lpt-bias-pk.png', dpi=300, bbox_inches='tight') - if nonlinear_model: assert np.isclose(Pk_gg.sum(), 489563565.9282355) assert np.isclose(Pk_gm.sum(), 425540686.4287062) else: assert np.isclose(Pk_gg.sum(), 493325841.6713596) assert np.isclose(Pk_gm.sum(), 429302962.1718302) + +@pytest.mark.parametrize("nonlinear_model", [True, False]) +def test_LPT_bias_compute_interpolator(nonlinear_model): + + skip_lpt = pytest.importorskip("velocileptors") # noqa F841 + from soliket.bias import LPT_bias + + info["theory"] = { + "camb": None, + "LPT_bias": {"external": LPT_bias, + "nonlinear": nonlinear_model} + } + + model = get_model(info) # noqa F841 + model.add_requirements({"Pk_grid": {"z": 0., "k_max": 1., + "nonlinear": True, + "vars_pairs": ('delta_tot', 'delta_tot'), + }, + "Pk_gg_interpolator": None, + "Pk_gm_interpolator": None, + "Pk_mm_interpolator": None + }) + + model.logposterior(info['params']) # force computation of model + + lhood = model.likelihood['one'] + + Pk_gg = lhood.provider.get_Pk_gg_interpolator().P(0.0, 1.e-2) + Pk_gm = lhood.provider.get_Pk_gm_interpolator().P(0.0, 1.e-2) + Pk_mm = lhood.provider.get_Pk_mm_interpolator().P(0.0, 1.e-2) + + if nonlinear_model: + assert np.isclose(Pk_gg, 85288.57867148393) + assert np.isclose(Pk_gm, 78672.33992060165) + assert np.isclose(Pk_mm, 78653.14602050644) + else: + assert np.isclose(Pk_gg, 85147.06509724) + assert np.isclose(Pk_gm, 78530.8263451) + assert np.isclose(Pk_mm, 78794.41444674769) From dc2e24d775a81b218128e30ce6b49b42cfd84730 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Tue, 10 May 2022 15:18:12 +0100 Subject: [PATCH 28/29] codestyle --- soliket/tests/test_bias.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/soliket/tests/test_bias.py b/soliket/tests/test_bias.py index 0708bf98..22e2037f 100644 --- a/soliket/tests/test_bias.py +++ b/soliket/tests/test_bias.py @@ -126,9 +126,10 @@ def test_LPT_bias_compute_grid(nonlinear_model): assert np.isclose(Pk_gg.sum(), 493325841.6713596) assert np.isclose(Pk_gm.sum(), 429302962.1718302) + @pytest.mark.parametrize("nonlinear_model", [True, False]) def test_LPT_bias_compute_interpolator(nonlinear_model): - + skip_lpt = pytest.importorskip("velocileptors") # noqa F841 from soliket.bias import LPT_bias From 705d3a795538524ca67169479f47a59821b74e96 Mon Sep 17 00:00:00 2001 From: Ian Harrison Date: Sun, 3 Jul 2022 15:12:15 +0100 Subject: [PATCH 29/29] use interpolators for power spectra --- soliket/xcorr/limber.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/soliket/xcorr/limber.py b/soliket/xcorr/limber.py index 153f263d..a57a378b 100644 --- a/soliket/xcorr/limber.py +++ b/soliket/xcorr/limber.py @@ -82,10 +82,14 @@ def do_limber(ell_arr, provider, dndz1, dndz2, s1, s2, pk, b1_HF, b2_HF, k_arr = (ell_arr + 0.5) / chi - p_mm_hf = pk(zatchi(chi), k_arr) - p_mm = p_mm_hf - p_gg = b1_HF * b1_HF * p_mm_hf # lets just stay at constant linear bias for now - p_gm = b1_HF * p_mm_hf + # p_mm_hf = pk(zatchi(chi), k_arr) + # p_mm = p_mm_hf + # p_gg = b1_HF * b1_HF * p_mm_hf # lets just stay at constant linear bias for now + # p_gm = b1_HF * p_mm_hf + + p_mm = self.provider.get_Pk_mm_interpolator().P(zatchi(chi), k_arr) + p_gg = self.provider.get_Pk_gg_interpolator().P(zatchi(chi), k_arr) + p_gm = self.provider.get_Pk_gm_interpolator().P(zatchi(chi), k_arr) W_g1g1 = W_g1[i_chi] * W_g1[i_chi] / (chi**2) * p_gg c_ell_g1g1[:, :, i_chi] = W_g1g1.T