diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 875463b2..51aed389 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -48,6 +48,19 @@ jobs: run: | wget https://lambda.gsfc.nasa.gov/data/map/dr5/dcp/likelihood/wmap_likelihood_full_v5.tar.gz + - name: Cache ACT data + id: cache-act + uses: actions/cache@v3 + with: + path: ACT_dr6_likelihood_v1.1.tgz + key: act-dr6-data + enableCrossOsArchive: true + + - if: ${{ steps.cache-act.outputs.cache-hit != 'true' }} + name: Download ACT data + run: | + wget https://lambda.gsfc.nasa.gov/data/suborbital/ACT/ACT_dr6/likelihood/data/ACT_dr6_likelihood_v1.1.tgz + @@ -74,8 +87,7 @@ jobs: - name: Install dependencies with conda shell: bash -l {0} - run: mamba install -c conda-forge "cosmosis>=3.0.1" cosmosis-build-standard-library - + run: mamba install -c conda-forge "cosmosis>=3.2" cosmosis-build-standard-library pytest - name: Get Cached Planck uses: actions/cache/restore@v3 @@ -90,17 +102,31 @@ jobs: path: wmap_likelihood_full_v5.tar.gz key: wmap-9-data fail-on-cache-miss: true + + - name: Get Cached ACT + uses: actions/cache/restore@v3 + with: + path: ACT_dr6_likelihood_v1.1.tgz + key: act-dr6-data + fail-on-cache-miss: true - name: Extract data shell: bash -l {0} run: | + # WMAP pushd likelihood/wmap9/data tar -zxvf ../../../wmap_likelihood_full_v5.tar.gz mv wmap_likelihood_v5/data/* . popd + #Planck pushd likelihood/planck2018 tar -zxvf ../../COM_Likelihood_Data-baseline_R3.00.tar.gz popd + # ACT + mkdir likelihood/act-dr6-lens/data + pushd likelihood/act-dr6-lens/data + tar -zxvf ../../../ACT_dr6_likelihood_v1.1.tgz + popd - name: Build standard library shell: bash -l {0} @@ -108,170 +134,17 @@ jobs: source cosmosis-configure make - - name: Smail example + - name: Run Tests shell: bash -l {0} run: | source cosmosis-configure - cosmosis examples/various-spectra.ini -p consistency.extra_relations='omega_x=omega_c+100' - grep -e 'omega_x = 100.261' output/various-spectra/cosmological_parameters/values.txt - - - name: BAO - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/bao.ini - - - name: CAMB Planck - shell: bash -l {0} - run: | - # print out the configuration - python -m cosmosis.configure - source cosmosis-configure - export - mkdir -p output - # And that the pipeline runs afterwards - cosmosis examples/planck.ini | tee output/planck.log - grep -e 'Likelihood = -1441.30' -e 'Likelihood = -1441.46' output/planck.log - - - - name: Class Planck - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/planck_class.ini -p class.mpk=T | tee output/class.log - # This may need updating as we modify the class interface - # The settings are not optimized - # grep 'Likelihood = -5968.93' output/class.log - # This seems to give different results on different systems - - - - name: WMAP likelihood - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/wmap.ini | tee output/wmap.log - - - - name: Pantheon emcee - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/pantheon.ini -p emcee.samples=20 - cosmosis-postprocess examples/pantheon.ini -o output/pantheon - test -f output/pantheon/cosmological_parameters--omega_m.png - - - name: PantheonPlusAndShoes - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/pantheon_plus_shoes.ini -p runtime.sampler=test | tee output/pantheon_plus_shoes_log.txt - grep 'Likelihood = -738.23' output/pantheon_plus_shoes_log.txt - - - name: DES Y1 likelihood - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y1.ini | tee output/des-y1.log - # unchanged in camb 1.4 - grep 'Likelihood = 5237.3' output/des-y1.log - - - name: DES Y1 likelihood with cl_to_corr - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y1.ini -p 2pt_shear.file=./shear/cl_to_corr/cl_to_corr.py 2pt_shear.corr_type=xi | tee output/des-y1.log - # unchanged in camb 1.4 - grep 'Likelihood = 5237.3' output/des-y1.log - - - - name: DES Y3 likelihood - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y3.ini -p pk_to_cl_gg.save_kernels=T pk_to_cl.save_kernels=T | tee output/des-y3.log - # Different versions - camb changes and scipy interpolation changes both alter these - grep -e 'Likelihood = 6043.23' -e 'Likelihood = 6043.34' -e 'Likelihood = 6043.33' output/des-y3.log - - - - name: DES Y3 Class likelihood - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y3-class.ini - # class is not consistent across systems to the level needed here - - - name: DES Y3 cosmic shear likelihood - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y3-shear.ini | tee output/des-y3-shear.log - grep -e 'Likelihood = 2957.03' -e 'Likelihood = 2957.12' -e 'Likelihood = 2957.11' output/des-y3-shear.log - - - - name: DES Y3 likelihood with Mira Titan - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y3-mira-titan.ini | tee output/des-y3-mt.log - grep -e 'Likelihood = 6048.0' -e 'Likelihood = 6048.1' output/des-y3-mt.log - - - name: DES Y3 likelihood with Mead nonlinear - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y3.ini -v halo_model_parameters.logT_AGN=8.2 -p camb.halofit_version=mead2020_feedback | tee output/des-y3-mead.log - grep -e 'Likelihood = 6049.94' -e 'Likelihood = 6049.00' output/des-y3-mead.log - - - name: Joint DES Y3 and KiDS-1000 likelihood - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y3_and_kids-1000.ini | tee output/des-y3_and_kids-1000.log - grep 'Likelihood = -199.409' output/des-y3_and_kids-1000.log - - - uses: actions/cache@v2 - name: ACT Data Cache - id: cache-act - with: - path: likelihood/act-dr6-lens/data/v1.1 - key: ${{ runner.os }}-act-dr6-v1.1 - - - name: Download ACT DR6 Lensing Data - if: steps.cache-act.outputs.cache-hit != 'true' - run: | - cd likelihood/act-dr6-lens/ - ./get-act-data.sh - - - name: ACT DR6 Lensing - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/act-dr6-lens.ini | tee output/act-dr6.log - grep -e 'Likelihood = -9.89' -e 'Likelihood = -9.90' output/act-dr6.log - - - name: DES Y3 6x2pt - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/des-y3-6x2.ini - - - name: Euclid emulator - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/euclid-emulator.ini - test -f output/euclid-emulator/matter_power_nl/p_k.txt - - - name: Log w example - shell: bash -l {0} - run: | - source cosmosis-configure - cosmosis examples/w_model.ini | tee output/w_model.log + pytest -vv --durations=0 tests/test_cosmosis_standard_library.py - name: Run Demos shell: bash -l {0} run: | - .github/run-demos.sh + source cosmosis-configure + pytest -vv --durations=0 tests/test_demos.py @@ -280,7 +153,7 @@ jobs: needs: Download_Data strategy: matrix: - python-version: [3.6, 3.7, 3.8, 3.9, "3.10", "3.11"] + python-version: [3.7, 3.8, 3.9, "3.10", "3.11"] steps: - name: Checkout repository @@ -315,16 +188,30 @@ jobs: key: wmap-9-data fail-on-cache-miss: true + - name: Get Cached ACT + uses: actions/cache/restore@v3 + with: + path: ACT_dr6_likelihood_v1.1.tgz + key: act-dr6-data + fail-on-cache-miss: true + - name: Extract data shell: bash -l {0} run: | - pushd likelihood/wmap9/data - tar -zxvf ../../../wmap_likelihood_full_v5.tar.gz - mv wmap_likelihood_v5/data/* . - popd - pushd likelihood/planck2018 - tar -zxvf ../../COM_Likelihood_Data-baseline_R3.00.tar.gz - popd + # WMAP + pushd likelihood/wmap9/data + tar -zxvf ../../../wmap_likelihood_full_v5.tar.gz + mv wmap_likelihood_v5/data/* . + popd + #Planck + pushd likelihood/planck2018 + tar -zxvf ../../COM_Likelihood_Data-baseline_R3.00.tar.gz + popd + # ACT + mkdir likelihood/act-dr6-lens/data + pushd likelihood/act-dr6-lens/data + tar -zxvf ../../../ACT_dr6_likelihood_v1.1.tgz + popd - uses: actions/cache@v2 @@ -338,7 +225,7 @@ jobs: - name: Install python dependencies run: | python -m pip install --upgrade pip wheel setuptools - pip install "cosmosis>=3.2" + pip install "cosmosis>=3.2" "nautilus-sampler==0.6.*" pip install -v --no-cache-dir --no-binary=mpi4py,camb mpi4py camb pip install fitsio astropy fast-pt "Cython<3.0" jupyter @@ -365,20 +252,6 @@ jobs: grep -e "🟢 fiducial-test-only" output/campaign.log cosmosis-campaign examples/des-campaign.yml --status - - name: Demo 2 Planck 2018 Likelihood - run: | - source .github/ci-setup.sh - cosmosis demos/demo2.ini - cosmosis-postprocess demos/demo2.ini -o output/demo2 - test -f output/demo2/comoving_distance.png - - - name: Demo 9 Low-Resolution Multinest - run: | - source .github/ci-setup.sh - cosmosis demos/demo9.ini -p multinest.live_points=100 - cosmosis-postprocess output/demo9.txt -o output/demo9 - test -f output/demo9/2D_supernova_params--deltam_cosmological_parameters--omega_m.png - - name: ExampleNotebook shell: bash -l {0} run: | diff --git a/boltzmann/camb/camb_interface.py b/boltzmann/camb/camb_interface.py index f4ff7d8d..cbc590f3 100644 --- a/boltzmann/camb/camb_interface.py +++ b/boltzmann/camb/camb_interface.py @@ -1,9 +1,9 @@ from cosmosis.datablock import names, option_section as opt from cosmosis.datablock.cosmosis_py import errors import numpy as np -from scipy.interpolate import InterpolatedUnivariateSpline import warnings import traceback +import contextlib # Finally we can now import camb import camb @@ -38,6 +38,15 @@ 'v_baryon_cdm': 'baryon_cdm_relative_velocity_power', } +@contextlib.contextmanager +def be_quiet_camb(): + original_feedback_level = camb.config.FeedbackLevel + try: + camb.set_feedback_level(0) + yield + finally: + camb.set_feedback_level(original_feedback_level) + def get_optional_params(block, section, names): params = {} @@ -349,12 +358,12 @@ def extract_camb_params(block, config, more_config): warnings.warn("Parameter omega_nu and omnuh2 are being ignored. Set mnu and num_massive_neutrinos instead.") # Set h if provided, otherwise look for theta_mc - if block.has_value(cosmo, "hubble"): + if block.has_value(cosmo, "cosmomc_theta"): + cosmology_params["cosmomc_theta"] = block[cosmo, "cosmomc_theta"] / 100 + elif block.has_value(cosmo, "hubble"): cosmology_params["H0"] = block[cosmo, "hubble"] - elif block.has_value(cosmo, "h0"): - cosmology_params["H0"] = block[cosmo, "h0"]*100 else: - cosmology_params["cosmomc_theta"] = block[cosmo, "cosmomc_theta"]/100 + cosmology_params["H0"] = block[cosmo, "h0"]*100 p = camb.CAMBparams( InitPower = init_power, @@ -366,11 +375,12 @@ def extract_camb_params(block, config, more_config): **config, ) # Setting up neutrinos by hand is hard. We let CAMB deal with it instead. - p.set_cosmology(ombh2 = block[cosmo, 'ombh2'], - omch2 = block[cosmo, 'omch2'], - omk = block[cosmo, 'omega_k'], - **more_config["cosmology_params"], - **cosmology_params) + with be_quiet_camb(): + p.set_cosmology(ombh2 = block[cosmo, 'ombh2'], + omch2 = block[cosmo, 'omch2'], + omk = block[cosmo, 'omega_k'], + **more_config["cosmology_params"], + **cosmology_params) # Fix for CAMB version < 1.0.10 if np.isclose(p.omnuh2, 0) and "nnu" in cosmology_params and not np.isclose(cosmology_params["nnu"], p.num_nu_massless): @@ -642,14 +652,17 @@ def execute(block, config): more_config["n_printed_errors"] += 1 return 1 - save_derived_parameters(r, block) - save_distances(r, block, more_config) + with be_quiet_camb(): + save_derived_parameters(r, block) + save_distances(r, block, more_config) if p.WantTransfer: - save_matter_power(r, block, more_config) + with be_quiet_camb(): + save_matter_power(r, block, more_config) if p.WantCls: - save_cls(r, block) + with be_quiet_camb(): + save_cls(r, block) return 0 diff --git a/examples/act-dr6-lens-values.ini b/examples/act-dr6-lens-values.ini index 820efbeb..f431c0bb 100644 --- a/examples/act-dr6-lens-values.ini +++ b/examples/act-dr6-lens-values.ini @@ -4,7 +4,7 @@ omch2 = 0.08 0.12 0.16 ombh2 = 0.01 0.02233 0.03 log1e10As = 2.9 3.0448 3.2 n_s = 0.86 0.96 1.06 -cosmomc_theta = 0.009 0.0104 0.012 +cosmomc_theta = 0.9 1.04 1.20 omega_k = 0.0 ;spatial curvature diff --git a/examples/planck.ini b/examples/planck.ini index 3e3547f2..945a7c7a 100644 --- a/examples/planck.ini +++ b/examples/planck.ini @@ -34,13 +34,14 @@ data_3 = %(planck_path)s/low_l/simall/simall_100x143_offlike5_EE_Aplanck_B.clik ; and any other that modules in the pipeline may want (e.g. camb) [consistency] file = ./utility/consistency/consistency_interface.py +cosmomc_theta = T [camb] file = boltzmann/camb/camb_interface.py mode = cmb lmax = 2800 ;max ell to use for cmb calculation -feedback=2 ;amount of output to print +feedback=1 ;amount of output to print AccuracyBoost=1.1 ;CAMB accuracy boost parameter do_tensors = True ;include tensor modes do_lensing = true ;lensing is required w/ Planck data diff --git a/examples/planck_class.ini b/examples/planck_class.ini index 015d595e..cd0e3579 100644 --- a/examples/planck_class.ini +++ b/examples/planck_class.ini @@ -34,6 +34,7 @@ data_3 = %(planck_path)s/low_l/simall/simall_100x143_offlike5_EE_Aplanck_B.clik ; and any other that modules in the pipeline may want (e.g. camb) [consistency] file = ./utility/consistency/consistency_interface.py +cosmomc_theta = T [class] file = boltzmann/class/class_interface.py diff --git a/examples/planck_values.ini b/examples/planck_values.ini index 94c4c1f5..ed647563 100644 --- a/examples/planck_values.ini +++ b/examples/planck_values.ini @@ -2,7 +2,7 @@ a_planck = 1.0 [cosmological_parameters] -h0 = 0.6732 ;H0 (km/s/Mpc)/100.0km/s/Mpc +cosmomc_theta = 1.040909 ; this is actually 100 * theta_mc omega_k = 0.0 ;spatial curvature ombh2 = 0.022383 omch2 = 0.12011 diff --git a/examples/wmap.ini b/examples/wmap.ini index b298d792..3bde1b3e 100644 --- a/examples/wmap.ini +++ b/examples/wmap.ini @@ -27,7 +27,7 @@ file = likelihood/wmap9/wmap_interface.so ; and any other that modules in the pipeline may want (e.g. camb) [consistency] file = ./utility/consistency/consistency_interface.py - +cosmomc_theta = T [camb] file = boltzmann/camb/camb_interface.py diff --git a/tests/test_cosmosis_standard_library.py b/tests/test_cosmosis_standard_library.py new file mode 100644 index 00000000..049c0c96 --- /dev/null +++ b/tests/test_cosmosis_standard_library.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python +from cosmosis import run_cosmosis +from cosmosis.postprocessing import run_cosmosis_postprocess +import pytest +import os + +def check_likelihood(capsys, expected, *other_possible): + captured = capsys.readouterr() + expect = (expected, *other_possible) + lines = [line for line in captured.out.split("\n") if "Likelihood =" in line] + print(lines) + lines = "\n".join(lines) + 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(): + 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() + + + +def test_bao(): + run_cosmosis("examples/bao.ini") + +def test_planck(capsys): + run_cosmosis("examples/planck.ini") + check_likelihood(capsys, "-1441.14", "-1441.30", "-1441.46") + +def test_planck_class(capsys): + run_cosmosis("examples/planck_class.ini", override={("class","mpk"):"T"}) + +def test_wmap(): + 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") + +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") + +def test_pantheon_plus_shoes(capsys): + run_cosmosis("examples/pantheon_plus_shoes.ini", override={("runtime","sampler"):"test"}) + check_likelihood(capsys, "-738.23") + +def test_des_y1(capsys): + run_cosmosis("examples/des-y1.ini") + check_likelihood(capsys, "5237.3") + +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") + +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") + +def test_des_y3_class(): + 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") + +def test_des_y3_mira_titan(capsys): + run_cosmosis("examples/des-y3-mira-titan.ini") + check_likelihood(capsys, "6048.0", "6048.1", "6048.2") + +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") + +def test_act_dr6_lensing(capsys): + run_cosmosis("examples/act-dr6-lens.ini") + check_likelihood(capsys, "-9.89", "-9.86", "-9.90") + +def test_des_y3_6x2pt(): + run_cosmosis("examples/des-y3-6x2.ini") + +def test_euclid_emulator(): + run_cosmosis("examples/euclid-emulator.ini") + assert os.path.exists("output/euclid-emulator/matter_power_nl/p_k.txt") + +def test_log_w_example(): + run_cosmosis("examples/w_model.ini") + +def test_des_kids(capsys): + run_cosmosis("examples/des-y3_and_kids-1000.ini") + check_likelihood(capsys, "-199.40", "-199.41") diff --git a/tests/test_demos.py b/tests/test_demos.py new file mode 100644 index 00000000..9e3e3967 --- /dev/null +++ b/tests/test_demos.py @@ -0,0 +1,112 @@ +import os +from cosmosis import run_cosmosis +from cosmosis.postprocessing import run_cosmosis_postprocess + + +def run_demo(i, args=None, variables=None): + if args is None: + args = [] + if variables is None: + variables = [] + + print("Running demo", i) + if i == 10: + return run_demo_10() + override = {("runtime", "verbosity"): "debug"} + for arg in args: + key, val = arg.split("=") + sec, key = key.split(".") + override[sec, key] = val + + variables_arg = {} + for arg in variables: + key, val = arg.split("=") + sec, key = key.split(".") + variables_arg[sec, key] = val + + run_cosmosis(f"demos/demo{i}.ini", override=override, variables=variables_arg) + plots = run_cosmosis_postprocess([f"demos/demo{i}.ini"], outdir=f"output/plots/demo{i}") + + if plots is not None: + plots.save() + +def run_demo_10(): + override = {("grid", "nsample_dimension"): "10"} + os.environ["HALOFIT"] = "takahashi" + run_cosmosis(f"demos/demo10.ini", override=override) + os.environ["HALOFIT"] = "mead2020" + run_cosmosis(f"demos/demo10.ini", override=override) + + plots = run_cosmosis_postprocess(["output/demo10_mead2020.txt", "output/demo10_takahashi.txt"], outdir="output/plots/demo10") + plots.save() + + +def test_demo1(): + run_demo(1) + +def test_demo2(): + run_demo(2) + +def test_demo3(): + run_demo(3, ["grid.nsample_dimension=10"]) + +def test_demo4(): + run_demo(4, variables=[ + "cosmological_parameters.n_s=0.962", + "cosmological_parameters.h0=0.680", + "cosmological_parameters.ombh2=0.0191", + ], args=[ + "maxlike.tolerance=1.0"]) + +def test_demo5(): + run_demo(5, ["emcee.samples=30", "emcee.nsteps=10", "emcee.walkers=20"]) + +def test_demo6(): + run_demo(6) + +def test_demo7(): + run_demo(7) + +def test_demo8(): + run_demo(8) + +def test_demo9(): + run_demo(9, ["multinest.live_points=50"]) + +def test_demo10(): + run_demo_10() + +def test_demo11(): + run_demo(11, ["grid.nsample_dimension=10"]) + +def test_demo12(): + run_demo(12) + +def test_demo13(): + run_demo(13, ["snake.threshold=3", "snake.nsample_dimension=10"]) + +def test_demo14(): + run_demo(14, ["zeus.samples=10", "zeus.nsteps=5"]) + +def test_demo15(): + run_demo(15) + +def test_demo16(): + run_demo(16) + +def test_demo17(): + run_demo(17, + variables=[ + "cosmological_parameters.n_s=0.96", + "cosmological_parameters.omega_b=0.0468", + "cosmological_parameters.h0=0.6881", + "intrinsic_alignment_parameters.alpha=1.0", + "intrinsic_alignment_parameters.a=1.0", + ], + args=["test.fatal_errors=T"]) + +def test_demo18(): + run_demo(18) + +def test_demo19(): + run_demo(19, ["grid.nsample_dimension=20"]) diff --git a/utility/consistency/consistency.py b/utility/consistency/consistency.py index 9803bdb2..e4bd8a89 100644 --- a/utility/consistency/consistency.py +++ b/utility/consistency/consistency.py @@ -58,7 +58,11 @@ def cosmology_consistency(verbose=False, relations_file="", theta=False, extra_r extra_relations = [rel.split('=', 1) for rel in extra_relations.replace(' ', '').split(',')] relations = relations + extra_relations - return Consistency(relations, COSMOLOGY_POSSIBLE_DEFAULTS, verbose) + + # These params are needed for the stupid CMB theta parameter. + # This was a lovely elegant design before we had to support that. + extra_fixed_params = ["YHe", "nnu", "num_massive_neutrinos", "TCMB"] + return Consistency(relations, COSMOLOGY_POSSIBLE_DEFAULTS, verbose, extra_fixed_params) COSMOLOGY_CONSISTENCY_RELATIONS = [ @@ -97,8 +101,8 @@ def cosmology_consistency(verbose=False, relations_file="", theta=False, extra_r ("omega_m", "1-omega_lambda-omega_k"), ("omega_k", "1-omega_m-omega_lambda"), ("K", "-hubble*hubble*omega_k/299792.458/299792.458"), - ("omnuh2", "mnu/93.14"), - ("mnu", "omnuh2*93.14"), + ("omnuh2", "mnu * ((nnu / 3.0) ** 0.75 / 94.06410581217612 * (TCMB/2.7255)**3)"), + ("mnu", "omnuh2 / ((nnu / 3.0) ** 0.75 / 94.06410581217612 * (TCMB/2.7255)**3)"), ] THETA_RELATIONS = [ @@ -129,12 +133,13 @@ class UnderSpecifiedModel(PoorlySpecifiedModel): class Consistency(object): - def __init__(self, relations, possible_defaults, verbose=False): + def __init__(self, relations, possible_defaults, verbose=False, extra_fixed=None): self.relations = relations self.parameters = {} self.possible_defaults = possible_defaults self.verbose = verbose self.reset() + self.extra_fixed = extra_fixed or [] self.cached_defaults = None def __call__(self, parameters): diff --git a/utility/consistency/consistency_interface.py b/utility/consistency/consistency_interface.py index 87dfe263..b15a5748 100644 --- a/utility/consistency/consistency_interface.py +++ b/utility/consistency/consistency_interface.py @@ -20,7 +20,12 @@ def execute(block, config): # Create dict of all parameters that we have already known_parameters = {} - for param in cons.parameters: + + # We need TCMB and nnu to get omnuh2 from mnu. If it's not specified use the default + block.get_double("cosmological_parameters", "TCMB", 2.7255) + block.get_double("cosmological_parameters", "nnu", 3.044) + + for param in list(cons.parameters.keys()) + cons.extra_fixed: if '___' in param: section, lookup_param = param.split('___') else: diff --git a/utility/consistency/theta_h0.py b/utility/consistency/theta_h0.py index ec03125d..76e5c901 100644 --- a/utility/consistency/theta_h0.py +++ b/utility/consistency/theta_h0.py @@ -1,32 +1,23 @@ -from astropy.cosmology import LambdaCDM -from astropy import units, constants +import camb import numpy as np import scipy.integrate import scipy.optimize -def dtauda(a, cosmo): - z = 1.0 / a - 1 - H = cosmo.H(z) - return (constants.c / (H * a**2)).to("Mpc").value -def dsound_da(a, cosmo): - ombh2 = cosmo.Ob0 * cosmo.h**2 - R = 3.0e4 * a * ombh2 - c_s = 1.0 / np.sqrt(3*(1+R)) - return dtauda(a, cosmo) * c_s -def H0_to_theta(hubble, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omega_k): + +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): h = hubble / 100. if np.isnan(omnuh2): omnuh2 = omega_nu * h**2 if np.isnan(omega_m): omega_m = ommh2 / h**2 - if np.isnan(omega_b): - omega_b = ombh2 / h**2 + if np.isnan(ombh2): + ombh2 = omega_b * h**2 if np.isnan(omega_lambda): omega_lambda = omlamh2 / h**2 - if np.isnan(omega_c): - omega_c = omch2 / h**2 + if np.isnan(omch2): + omch2 = omega_c * h**2 if np.isnan(omega_m): omega_m = omega_c + omega_b + omega_nu if np.isnan(omega_m): @@ -34,20 +25,31 @@ 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, omega_b]).any(): + if np.isnan([omnuh2, hubble, omega_m, omega_lambda, ombh2]).any(): return np.nan - m_nu = 93.14 * omnuh2 * units.eV - cosmo = LambdaCDM(hubble, omega_m, omega_lambda, m_nu=m_nu/3, Ob0=omega_b, Tcmb0=2.7255, Neff=3.046) + mnu = 93.14 * omnuh2 - ombh2 = cosmo.Ob0 * cosmo.h**2 - omdmh2 = cosmo.Om0 * cosmo.h**2 + original_feedback_level = camb.config.FeedbackLevel - zstar = 1048*(1+0.00124*ombh2**(-0.738))*(1+ (0.0783*ombh2**(-0.238)/(1+39.5*ombh2**0.763)) * (omdmh2+ombh2)**(0.560/(1+21.1*ombh2**1.81))) - astar = 1 / (1+zstar) - rs = scipy.integrate.romberg(dsound_da, 1e-8, astar, rtol=5e-5, args=(cosmo,), divmax=10, vec_func=True) # in Mpc - DA = cosmo.angular_diameter_distance(zstar).to("Mpc").value / astar - return rs / DA + try: + camb.set_feedback_level(0) + p = camb.CAMBparams() + p.set_cosmology(ombh2 = ombh2, + omch2 = omch2, + omk = omega_k, + mnu = mnu, + num_massive_neutrinos=num_massive_neutrinos, + nnu=nnu, + H0=hubble) + 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): @@ -63,48 +65,46 @@ 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) + num_massive_neutrinos = params.get('num_massive_neutrinos', 1) + nnu = params.get('nnu', 3.044) - return H0_to_theta(hubble, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omega_k) - + 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) -def theta_to_H0(theta, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omega_k): - get_theta = lambda H0: H0_to_theta(H0, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omega_k) - t = get_theta(70.0) - - if np.isnan(t): +def theta_to_H0(theta, omnuh2, omch2, ombh2, omega_k, num_massive_neutrinos, nnu): + if np.isnan([ombh2, omch2, theta, omega_k, omnuh2]).any(): return np.nan + mnu = 93.14 * omnuh2 - H0_min = 20.0 - H0_max = 120.0 - - try: - theta_min = get_theta(H0_min) - except ValueError: - H0_min = 40.0 - + original_feedback_level = camb.config.FeedbackLevel try: - theta_max = get_theta(H0_max) - except ValueError: - H0_max = 100.0 + camb.set_feedback_level(0) + p = camb.CAMBparams() + p.set_cosmology(ombh2 = ombh2, + omch2 = omch2, + omk = omega_k, + mnu = mnu, + num_massive_neutrinos=num_massive_neutrinos, + nnu=nnu, + cosmomc_theta = theta) + H0 = p.h * 100 + except: + H0 = np.nan + finally: + camb.config.FeedbackLevel = original_feedback_level - f = lambda H0: get_theta(H0) - theta - H0 = scipy.optimize.brentq(f, H0_min, H0_max, rtol=5e-5) return H0 -def theta_to_H0_interface(params): - theta = params['cosmomc_theta'] - omega_nu = params.get('omega_nu', np.nan) +def theta_to_H0_interface(params): + theta = params['cosmomc_theta'] / 100 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) omch2 = params.get('omch2', np.nan) - omega_b = params.get('omega_b', np.nan) ombh2 = params.get('ombh2', np.nan) - omega_lambda = params.get('omega_lambda', np.nan) - omlamh2 = params.get('omlamh2', np.nan) omegak = params.get('omega_k', np.nan) - return theta_to_H0(theta, omega_nu, omnuh2, omega_m, ommh2, omega_c, omch2, omega_b, ombh2, omega_lambda, omlamh2, omegak) + # 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)