Merge pull request #103 from joezuntz/fix-theta-w

Checks on theta for w models

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)

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

tests/test_cosmosis_standard_library.py

@@ -13,94 +13,123 @@ 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={
         ("2pt_shear","file"): "./shear/cl_to_corr/cl_to_corr.py",
         ("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={
         ("pk_to_cl_gg","save_kernels"):"T",
         ("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", 
                  override={("camb", "halofit_version"): "mead2020_feedback"},
                  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):
+        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")
+    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)

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)
 
 

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

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,20 +44,22 @@ 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
 
 
 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)
     omega_m = params.get('omega_m', np.nan)
     ommh2 = params.get('ommh2', np.nan)
     omega_c = params.get('omega_c', np.nan)
@@ -65,16 +69,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 +95,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 +108,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)