From da453007b57017924e41c371d801d6f5d52954a7 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Wed, 30 Aug 2023 12:22:03 +0100 Subject: [PATCH 1/4] Checks on theta for w models --- boltzmann/camb/camb_interface.py | 5 +++ tests/test_cosmosis_standard_library.py | 4 +++ utility/consistency/consistency.py | 6 ++-- utility/consistency/consistency_interface.py | 13 ++++++-- utility/consistency/theta_h0.py | 34 ++++++++++++++------ 5 files changed, 48 insertions(+), 14 deletions(-) diff --git a/boltzmann/camb/camb_interface.py b/boltzmann/camb/camb_interface.py index cbc590f3..83fec043 100644 --- a/boltzmann/camb/camb_interface.py +++ b/boltzmann/camb/camb_interface.py @@ -251,6 +251,11 @@ def extract_dark_energy_params(block, config, more_config): dark_energy = de_class() if more_config['use_tabulated_w']: + if block.has_value(cosmo, "consistency_module_was_used") and block.has_value(cosmo, "cosmomc_theta"): + raise RuntimeError("You used the consistency module with cosmomc_theta=T but are also" + "using a tabulated w(a) in camb. The theta-H0 relation as implemeted" + "in the consistency module will not work for models other than w0-wa" + ) a = block[names.de_equation_of_state, 'a'] w = block[names.de_equation_of_state, 'w'] dark_energy.set_w_a_table(a, w) diff --git a/tests/test_cosmosis_standard_library.py b/tests/test_cosmosis_standard_library.py index d0588702..2bc6a886 100644 --- a/tests/test_cosmosis_standard_library.py +++ b/tests/test_cosmosis_standard_library.py @@ -96,6 +96,10 @@ def test_euclid_emulator(): def test_log_w_example(): run_cosmosis("examples/w_model.ini") +def test_theta_warning(): + with pytest.raises(RuntimeError): + run_cosmosis("examples/w_model.ini", override={("consistency","cosmomc_theta"):"T"}) + def test_des_kids(capsys): run_cosmosis("examples/des-y3_and_kids-1000.ini") check_likelihood(capsys, "-199.40", "-199.41") diff --git a/utility/consistency/consistency.py b/utility/consistency/consistency.py index e4bd8a89..1a5e3e1b 100644 --- a/utility/consistency/consistency.py +++ b/utility/consistency/consistency.py @@ -59,9 +59,11 @@ def cosmology_consistency(verbose=False, relations_file="", theta=False, extra_r for rel in extra_relations.replace(' ', '').split(',')] relations = relations + extra_relations - # These params are needed for the stupid CMB theta parameter. - # This was a lovely elegant design before we had to support that. + # To relate mnu and omnuh2 we need these parameters. extra_fixed_params = ["YHe", "nnu", "num_massive_neutrinos", "TCMB"] + # If we are relating H0 and the CosmoMC theta parameter then we need these too. + if theta: + extra_fixed_params += ["w", "wa"] return Consistency(relations, COSMOLOGY_POSSIBLE_DEFAULTS, verbose, extra_fixed_params) diff --git a/utility/consistency/consistency_interface.py b/utility/consistency/consistency_interface.py index b15a5748..e4b27e9b 100644 --- a/utility/consistency/consistency_interface.py +++ b/utility/consistency/consistency_interface.py @@ -11,11 +11,11 @@ def setup(options): option_section, "relations_file", default="") extra_relations = options.get_string(option_section, "extra_relations", default="") cons = consistency.cosmology_consistency(verbose, relations_file, cosmomc_theta, extra_relations) - return cons + return cons, cosmomc_theta def execute(block, config): - cons = config + cons, cosmomc_theta = config cosmo = names.cosmological_parameters # Create dict of all parameters that we have already @@ -25,6 +25,12 @@ def execute(block, config): block.get_double("cosmological_parameters", "TCMB", 2.7255) block.get_double("cosmological_parameters", "nnu", 3.044) + # We need these to get cosmomc_theta from H0 or vice versa, so only get them + # if cosmomc_theta mode is active + if cosmomc_theta: + block.get_double("cosmological_parameters", "w", -1.0) + block.get_double("cosmological_parameters", "wa", 0.0) + for param in list(cons.parameters.keys()) + cons.extra_fixed: if '___' in param: section, lookup_param = param.split('___') @@ -89,4 +95,7 @@ def execute(block, config): block[cosmo, "sigma_8"] = sigma_8 + # Add a marker to the consistency + block[cosmo, "consistency_module_was_used"] = True + return 0 diff --git a/utility/consistency/theta_h0.py b/utility/consistency/theta_h0.py index 76e5c901..8c0d3a99 100644 --- a/utility/consistency/theta_h0.py +++ b/utility/consistency/theta_h0.py @@ -6,7 +6,7 @@ -def H0_to_theta(hubble, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omega_k, num_massive_neutrinos, nnu): +def H0_to_theta(hubble, omega_nu, omnuh2, mnu, TCMB, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omega_k, num_massive_neutrinos, nnu, w, wa): h = hubble / 100. if np.isnan(omnuh2): omnuh2 = omega_nu * h**2 @@ -25,13 +25,15 @@ def H0_to_theta(hubble, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_ if np.isnan(omega_lambda): omega_lambda = 1 - omega_k - omega_m - if np.isnan([omnuh2, hubble, omega_m, omega_lambda, ombh2]).any(): + if np.isnan(mnu): + mnu = omnuh2 / ((nnu / 3.0) ** 0.75 / 94.06410581217612 * (TCMB/2.7255)**3) + + if np.isnan([hubble, omega_m, omega_lambda, ombh2, mnu]).any(): return np.nan - - mnu = 93.14 * omnuh2 original_feedback_level = camb.config.FeedbackLevel + try: camb.set_feedback_level(0) p = camb.CAMBparams() @@ -42,13 +44,13 @@ def H0_to_theta(hubble, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_ num_massive_neutrinos=num_massive_neutrinos, nnu=nnu, H0=hubble) + p.set_dark_energy(w=w, wa=wa, dark_energy_model='ppf') r = camb.get_background(p) theta = r.cosmomc_theta() except: theta = np.nan finally: camb.config.FeedbackLevel = original_feedback_level - return theta @@ -56,6 +58,9 @@ def H0_to_theta_interface(params): hubble = params['hubble'] omega_nu = params.get('omega_nu', np.nan) omnuh2 = params.get('omnuh2', np.nan) + mnu = params.get('mnu', np.nan) + TCMB = params.get('TCMB', np.nan) + omnuh2 = params.get('omnuh2', np.nan) omega_m = params.get('omega_m', np.nan) ommh2 = params.get('ommh2', np.nan) omega_c = params.get('omega_c', np.nan) @@ -65,16 +70,20 @@ def H0_to_theta_interface(params): omega_lambda = params.get('omega_lambda', np.nan) omlamh2 = params.get('omlamh2', np.nan) omega_k = params.get('omega_k', np.nan) + w = params.get('w', -1) + wa = params.get('wa', 0) num_massive_neutrinos = params.get('num_massive_neutrinos', 1) nnu = params.get('nnu', 3.044) - return 100 * H0_to_theta(hubble, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omega_k, num_massive_neutrinos, nnu) + return 100 * H0_to_theta(hubble, omega_nu, omnuh2, mnu, TCMB, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omega_k, num_massive_neutrinos, nnu, w, wa) + +def theta_to_H0(theta, omnuh2, mnu, TCMB, omch2, ombh2, omega_k, num_massive_neutrinos, nnu, w, wa): + if np.isnan(mnu): + mnu = omnuh2 / ((nnu / 3.0) ** 0.75 / 94.06410581217612 * (TCMB/2.7255)**3) -def theta_to_H0(theta, omnuh2, omch2, ombh2, omega_k, num_massive_neutrinos, nnu): - if np.isnan([ombh2, omch2, theta, omega_k, omnuh2]).any(): + if np.isnan([ombh2, omch2, theta, omega_k, mnu]).any(): return np.nan - mnu = 93.14 * omnuh2 original_feedback_level = camb.config.FeedbackLevel try: @@ -87,6 +96,7 @@ def theta_to_H0(theta, omnuh2, omch2, ombh2, omega_k, num_massive_neutrinos, nnu num_massive_neutrinos=num_massive_neutrinos, nnu=nnu, cosmomc_theta = theta) + p.set_dark_energy(w=w, wa=wa, dark_energy_model='ppf') H0 = p.h * 100 except: H0 = np.nan @@ -99,12 +109,16 @@ def theta_to_H0(theta, omnuh2, omch2, ombh2, omega_k, num_massive_neutrinos, nnu def theta_to_H0_interface(params): theta = params['cosmomc_theta'] / 100 omnuh2 = params.get('omnuh2', np.nan) + mnu = params.get('mnu', np.nan) + TCMB = params.get('TCMB', np.nan) omch2 = params.get('omch2', np.nan) ombh2 = params.get('ombh2', np.nan) omegak = params.get('omega_k', np.nan) + w = params.get('w', -1.0) + wa = params.get('wa', 0.0) # These are the camb defaults num_massive_neutrinos = params.get('num_massive_neutrinos', 1) nnu = params.get('nnu', 3.044) - return theta_to_H0(theta, omnuh2, omch2, ombh2, omegak, num_massive_neutrinos, nnu) + return theta_to_H0(theta, omnuh2, mnu, TCMB, omch2, ombh2, omegak, num_massive_neutrinos, nnu, w, wa) From a5a7ff13ae6b4000cd0e32581fe00ba9a95a2db5 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Wed, 30 Aug 2023 16:34:49 +0100 Subject: [PATCH 2/4] remove duplicate param --- examples/planck_values.ini | 3 --- 1 file changed, 3 deletions(-) diff --git a/examples/planck_values.ini b/examples/planck_values.ini index ed647563..e39a85ce 100644 --- a/examples/planck_values.ini +++ b/examples/planck_values.ini @@ -1,6 +1,3 @@ -[planck] -a_planck = 1.0 - [cosmological_parameters] cosmomc_theta = 1.040909 ; this is actually 100 * theta_mc omega_k = 0.0 ;spatial curvature From 55fee545e20702f53cf901a2f6c62c9321a1eadd Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Thu, 31 Aug 2023 13:10:19 +0100 Subject: [PATCH 3/4] remove double line --- utility/consistency/theta_h0.py | 1 - 1 file changed, 1 deletion(-) diff --git a/utility/consistency/theta_h0.py b/utility/consistency/theta_h0.py index 8c0d3a99..89c2835e 100644 --- a/utility/consistency/theta_h0.py +++ b/utility/consistency/theta_h0.py @@ -60,7 +60,6 @@ def H0_to_theta_interface(params): omnuh2 = params.get('omnuh2', np.nan) mnu = params.get('mnu', np.nan) TCMB = params.get('TCMB', np.nan) - omnuh2 = params.get('omnuh2', np.nan) omega_m = params.get('omega_m', np.nan) ommh2 = params.get('ommh2', np.nan) omega_c = params.get('omega_c', np.nan) From 97a3d3d196802d87f2d2e00e450012952727d526 Mon Sep 17 00:00:00 2001 From: Joe Zuntz Date: Thu, 31 Aug 2023 13:21:00 +0100 Subject: [PATCH 4/4] more testing --- tests/test_cosmosis_standard_library.py | 39 ++++++++++++++++++++----- 1 file changed, 32 insertions(+), 7 deletions(-) diff --git a/tests/test_cosmosis_standard_library.py b/tests/test_cosmosis_standard_library.py index 2bc6a886..b215a188 100644 --- a/tests/test_cosmosis_standard_library.py +++ b/tests/test_cosmosis_standard_library.py @@ -13,41 +13,55 @@ def check_likelihood(capsys, expected, *other_possible): msg = f"Likelihood was expected to be one of {expect} but this was not found. Found these lines: \n{lines}" assert any([f"Likelihood = {val}" in lines for val in (expected, *other_possible)]), msg -def test_projection(): +def check_no_camb_warnings(capsys): + captured = capsys.readouterr() + lines = [line for line in captured.out.split("\n") if "UserWarning: Parameter" in line] + assert len(lines)==0, f"Found some warnings from CAMB: {lines}" + + +def test_projection(capsys): run_cosmosis("examples/various-spectra.ini", override={("consistency","extra_relations"):"omega_x=omega_c+100"}) with open("output/various-spectra/cosmological_parameters/values.txt") as f: assert "omega_x = 100.261" in f.read() + check_no_camb_warnings(capsys) -def test_bao(): +def test_bao(capsys): run_cosmosis("examples/bao.ini") + check_no_camb_warnings(capsys) def test_planck(capsys): run_cosmosis("examples/planck.ini") check_likelihood(capsys, "-1441.14", "-1441.30", "-1441.46") + check_no_camb_warnings(capsys) def test_planck_class(capsys): run_cosmosis("examples/planck_class.ini", override={("class","mpk"):"T"}) + check_no_camb_warnings(capsys) -def test_wmap(): +def test_wmap(capsys): if not os.path.exists("likelihood/wmap9/data/lowlP/mask_r3_p06_jarosik.fits"): pytest.skip("WMAP data not found") run_cosmosis("examples/wmap.ini") + check_no_camb_warnings(capsys) def test_pantheon_emcee(capsys): run_cosmosis("examples/pantheon.ini", override={("emcee","samples"):"20"}) plots = run_cosmosis_postprocess(["examples/pantheon.ini"], outdir="output/pantheon") plots.save() assert os.path.exists("output/pantheon/cosmological_parameters--omega_m.png") + check_no_camb_warnings(capsys) def test_pantheon_plus_shoes(capsys): run_cosmosis("examples/pantheon_plus_shoes.ini", override={("runtime","sampler"):"test"}) check_likelihood(capsys, "-738.23") + check_no_camb_warnings(capsys) def test_des_y1(capsys): run_cosmosis("examples/des-y1.ini") check_likelihood(capsys, "5237.3") + check_no_camb_warnings(capsys) def test_des_y1_cl_to_corr(capsys): run_cosmosis("examples/des-y1.ini", override={ @@ -55,6 +69,7 @@ def test_des_y1_cl_to_corr(capsys): ("2pt_shear","corr_type"): "xi" }) check_likelihood(capsys, "5237.3") + check_no_camb_warnings(capsys) def test_des_y3(capsys): run_cosmosis("examples/des-y3.ini", override={ @@ -62,18 +77,21 @@ def test_des_y3(capsys): ("pk_to_cl","save_kernels"):"T" }) check_likelihood(capsys, "6043.23", "6043.34", "6043.37", "6043.33") + check_no_camb_warnings(capsys) -def test_des_y3_class(): +def test_des_y3_class(capsys): run_cosmosis("examples/des-y3-class.ini") # class is not consistent across systems to the level needed here def test_des_y3_shear(capsys): run_cosmosis("examples/des-y3-shear.ini") check_likelihood(capsys, "2957.03", "2957.12", "2957.11", "2957.13") + check_no_camb_warnings(capsys) def test_des_y3_mira_titan(capsys): run_cosmosis("examples/des-y3-mira-titan.ini") check_likelihood(capsys, "6048.0", "6048.1", "6048.2") + check_no_camb_warnings(capsys) def test_des_y3_mead(capsys): run_cosmosis("examples/des-y3.ini", @@ -81,20 +99,25 @@ def test_des_y3_mead(capsys): variables={("halo_model_parameters", "logT_AGN"): "8.2"} ) check_likelihood(capsys, "6049.94", "6049.00", "6049.03", "6049.04") + check_no_camb_warnings(capsys) def test_act_dr6_lensing(capsys): run_cosmosis("examples/act-dr6-lens.ini") check_likelihood(capsys, "-9.89", "-9.86", "-9.90") + check_no_camb_warnings(capsys) -def test_des_y3_6x2pt(): +def test_des_y3_6x2pt(capsys): run_cosmosis("examples/des-y3-6x2.ini") + check_no_camb_warnings(capsys) -def test_euclid_emulator(): +def test_euclid_emulator(capsys): run_cosmosis("examples/euclid-emulator.ini") assert os.path.exists("output/euclid-emulator/matter_power_nl/p_k.txt") + check_no_camb_warnings(capsys) -def test_log_w_example(): +def test_log_w_example(capsys): run_cosmosis("examples/w_model.ini") + check_no_camb_warnings(capsys) def test_theta_warning(): with pytest.raises(RuntimeError): @@ -103,8 +126,10 @@ def test_theta_warning(): def test_des_kids(capsys): run_cosmosis("examples/des-y3_and_kids-1000.ini") check_likelihood(capsys, "-199.40", "-199.41") + check_no_camb_warnings(capsys) def test_kids(capsys): run_cosmosis("examples/kids-1000.ini") check_likelihood(capsys, "-47.6") + check_no_camb_warnings(capsys)