Merge branch 'joezuntz:main' into main

Author: aferte

boltzmann/camb/camb_interface.py

@@ -122,6 +122,8 @@ def setup(options):
     config['want_zdrag'] = mode != MODE_BG
     config['want_zstar'] = config['want_zdrag']
 
+    more_config['recombination_model'] = options.get_string(opt, 'recombination_model', default='recfast')
+
     more_config['want_chistar'] = options.get_bool(opt, 'want_chistar', default=True)
     more_config['n_logz'] = options.get_int(opt, 'n_logz', default=0)
     more_config['zmax_logz'] = options.get_double(opt, 'zmax_logz', default = 1100.)
@@ -157,7 +159,7 @@ def setup(options):
         raise ValueError("halofit_version must be one of : {}.  You put: {}".format(known_halofit_versions, halofit_version))
 
     more_config["accuracy_params"] = get_optional_params(options, opt, 
-                                                        ['AccuracyBoost', 'lSampleBoost', 'lAccuracyBoost', 'DoLateRadTruncation'])
+                                                        ['AccuracyBoost', 'lSampleBoost', 'lAccuracyBoost', 'DoLateRadTruncation', 'min_l_logl_sampling'])
                                                         #  'TimeStepBoost', 'BackgroundTimeStepBoost', 'IntTolBoost', 
                                                         #  'SourcekAccuracyBoost', 'IntkAccuracyBoost', 'TransferkBoost',
                                                         #  'NonFlatIntAccuracyBoost', 'BessIntBoost', 'LensingBoost',
@@ -205,35 +207,39 @@ def setup(options):
 # during the execute function
 
 def extract_recombination_params(block, config, more_config):
-    default_recomb = camb.recombination.Recfast()
- 
-    min_a_evolve_Tm = block.get_double('recfast', 'min_a_evolve_Tm', default=default_recomb.min_a_evolve_Tm)
-    RECFAST_fudge = block.get_double('recfast', 'RECFAST_fudge', default=default_recomb.RECFAST_fudge)
-    RECFAST_fudge_He = block.get_double('recfast', 'RECFAST_fudge_He', default=default_recomb.RECFAST_fudge_He)
-    RECFAST_Heswitch = block.get_int('recfast', 'RECFAST_Heswitch', default=default_recomb.RECFAST_Heswitch)
-    RECFAST_Hswitch = block.get_bool('recfast', 'RECFAST_Hswitch', default=default_recomb.RECFAST_Hswitch)
-    AGauss1 = block.get_double('recfast', 'AGauss1', default=default_recomb.AGauss1)
-    AGauss2 = block.get_double('recfast', 'AGauss2', default=default_recomb.AGauss2)
-    zGauss1 = block.get_double('recfast', 'zGauss1', default=default_recomb.zGauss1)
-    zGauss2 = block.get_double('recfast', 'zGauss2', default=default_recomb.zGauss2)
-    wGauss1 = block.get_double('recfast', 'wGauss1', default=default_recomb.wGauss1)
-    wGauss2 = block.get_double('recfast', 'wGauss2', default=default_recomb.wGauss2)
-    
-    recomb = camb.recombination.Recfast(
-        min_a_evolve_Tm = min_a_evolve_Tm, 
-        RECFAST_fudge = RECFAST_fudge, 
-        RECFAST_fudge_He = RECFAST_fudge_He, 
-        RECFAST_Heswitch = RECFAST_Heswitch, 
-        RECFAST_Hswitch = RECFAST_Hswitch, 
-        AGauss1 = AGauss1, 
-        AGauss2 = AGauss2, 
-        zGauss1 = zGauss1, 
-        zGauss2 = zGauss2, 
-        wGauss1 = wGauss1, 
-        wGauss2 = wGauss2, 
-    )
+    if more_config['recombination_model'] == 'recfast':
+        default_recomb = camb.recombination.Recfast()
+     
+        min_a_evolve_Tm = block.get_double('recfast', 'min_a_evolve_Tm', default=default_recomb.min_a_evolve_Tm)
+        RECFAST_fudge = block.get_double('recfast', 'RECFAST_fudge', default=default_recomb.RECFAST_fudge)
+        RECFAST_fudge_He = block.get_double('recfast', 'RECFAST_fudge_He', default=default_recomb.RECFAST_fudge_He)
+        RECFAST_Heswitch = block.get_int('recfast', 'RECFAST_Heswitch', default=default_recomb.RECFAST_Heswitch)
+        RECFAST_Hswitch = block.get_bool('recfast', 'RECFAST_Hswitch', default=default_recomb.RECFAST_Hswitch)
+        AGauss1 = block.get_double('recfast', 'AGauss1', default=default_recomb.AGauss1)
+        AGauss2 = block.get_double('recfast', 'AGauss2', default=default_recomb.AGauss2)
+        zGauss1 = block.get_double('recfast', 'zGauss1', default=default_recomb.zGauss1)
+        zGauss2 = block.get_double('recfast', 'zGauss2', default=default_recomb.zGauss2)
+        wGauss1 = block.get_double('recfast', 'wGauss1', default=default_recomb.wGauss1)
+        wGauss2 = block.get_double('recfast', 'wGauss2', default=default_recomb.wGauss2)
+        
+        recomb = camb.recombination.Recfast(
+            min_a_evolve_Tm = min_a_evolve_Tm, 
+            RECFAST_fudge = RECFAST_fudge, 
+            RECFAST_fudge_He = RECFAST_fudge_He, 
+            RECFAST_Heswitch = RECFAST_Heswitch, 
+            RECFAST_Hswitch = RECFAST_Hswitch, 
+            AGauss1 = AGauss1, 
+            AGauss2 = AGauss2, 
+            zGauss1 = zGauss1, 
+            zGauss2 = zGauss2, 
+            wGauss1 = wGauss1, 
+            wGauss2 = wGauss2, 
+        )
+
+    #Not yet supporting CosmoRec by default, but not too hard to compile in camb yourself if needed.
+    elif more_config['recombination_model'] == 'cosmorec':
 
-    #Not yet supporting CosmoRec, but not too hard if needed.
+        recomb = camb.recombination.CosmoRec()
 
     return recomb
 
@@ -579,9 +585,11 @@ def save_matter_power(r, block, more_config):
         # Note that the transfer function has different k range and precision as
         # there is no interpolating function for it in CAMB.
         # We might consider creating an interpolator ...
+        # Note here, that like in the returned sigma_8 and fsigma_8, the redshifts here are reversed
+        # thus for a redshift=0 transfer function we ask for the last element
         section_name_transfer = matter_power_section_names[transfer_type] + "_transfer_func"
-        transfer_func = r.get_matter_transfer_data().transfer_z(tt)
-        k_transfer_func = r.get_matter_transfer_data().transfer_z("k/h")
+        transfer_func = r.get_matter_transfer_data().transfer_z(tt, -1)
+        k_transfer_func = r.get_matter_transfer_data().transfer_z("k/h", -1)
         block.put_double_array_1d(section_name_transfer, "k_h", k_transfer_func)
         block.put_double_array_1d(section_name_transfer, "t_k", transfer_func)
 

boltzmann/camb/module.yaml

@@ -219,6 +219,10 @@ params:
         meaning: Whether to use the tau parametrization of reionization instead of z.
         type: bool
         default: true
+    recombination_model:
+        meaning: Which model to use for recombination physics. Choose from recfast, cosmorec
+        type: str
+        default: "recfast"
     DoLateRadTruncation:
         meaning: Whether to use smooth approx to radiation perturbations after decoupling
             on small scales, saving evolution of irrelevant osciallatory multipole

examples/act-dr6-lite.ini

@@ -9,16 +9,24 @@ save_dir=output/act
 fatal_errors=T
 
 [output]
-filename = output/act-lite.txt
+filename = output/act-lite-default.txt
 
 [maxlike]
 maxiter = 1000
 tolerance = 0.01
 max_posterior = T
+output_ini = output/act-lite-default-maxlike.ini
 
 [nautilus]
 n_live = 2000
 
+[emcee]
+walkers = 64
+samples = 2000
+nsteps = 10
+
+[grid]
+nsample_dimension = 30
 
 [pipeline]
 ; these names refer to sections later in the file:
@@ -27,6 +35,7 @@ values = examples/act-lite-values.ini
 priors = examples/act-lite-priors.ini
 debug=T
 timing=F
+extra_output = cosmological_parameters/sigma_8 cosmological_parameters/H_0
 
 
 [act]
@@ -43,11 +52,23 @@ cosmomc_theta = T
 
 
 [camb]
+; These accuracy settings are required to reproduce the published result
+; NOTE: cosmorec is also required for accurate reproduction, but not included in
+; supplied version of CAMB (so you should compile your own following instructions
+; in the CAMB documentation).
 file = boltzmann/camb/camb_interface.py
 mode = cmb
+kmax = 10.
+k_per_logint = 130
+nonlinear = true
+lens_potential_accuracy = 8
+lens_margin = 2050
+lAccuracyBoost = 1.2
+min_l_logl_sampling =  6000
+DoLateRadTruncation =  false
+recombination_model =  recfast ; cosmorec required for accurate reproduction
 lmax = 8500          ;max ell to use for cmb calculation
 feedback=0         ;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
 NonLinear = lens

examples/act-dr6.ini

@@ -9,7 +9,7 @@ save_dir=output/act
 fatal_errors=T
 
 [output]
-filename = output/act.txt
+filename = output/act-default.txt
 
 [maxlike]
 maxiter = 1000
@@ -19,6 +19,15 @@ max_posterior = T
 [nautilus]
 n_live = 2000
 
+[maxlike]
+maxiter = 1000
+tolerance = 0.01
+max_posterior = T
+
+[emcee]
+walkers = 128
+samples = 2000
+nsteps = 100
 
 [pipeline]
 ; these names refer to sections later in the file:
@@ -27,6 +36,7 @@ values = examples/act-values.ini
 priors = examples/act-priors.ini
 debug=T
 timing=F
+extra_output = cosmological_parameters/sigma_8 cosmological_parameters/H_0
 
 
 [act]
@@ -43,14 +53,26 @@ cosmomc_theta = T
 
 
 [camb]
+
+; These accuracy settings are required to reproduce the published result
+; NOTE: cosmorec is also required for accurate reproduction, but not included in
+; supplied version of CAMB (so you should compile your own following instructions
+; in the CAMB documentation).
 file = boltzmann/camb/camb_interface.py
 mode = cmb
+kmax = 10.
+k_per_logint = 130
+nonlinear = true
+lens_potential_accuracy = 8
+lens_margin = 2050
+lAccuracyBoost = 1.2
+min_l_logl_sampling =  6000
+DoLateRadTruncation =  false
+recombination_model =  recfast ; cosmorec required for accurate reproduction
 lmax = 8500          ;max ell to use for cmb calculation
 feedback=0         ;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
 NonLinear = lens
 accurate_massive_neutrino_transfers = T
-theta_H0_range = "20 100"
-
+theta_H0_range = "20 100"

examples/act-lite-priors.ini

@@ -1,3 +1,6 @@
 [act_params]
-A_act = gaussian 1.0 0.1 
-P_act = gaussian 1.0 0.1
+A_act = gaussian 1.0 0.003 
+P_act = gaussian 1.0 0.1
+
+[cosmological_parameters]
+tau = gaussian 0.0544 0.0073

examples/act-priors.ini

@@ -17,21 +17,21 @@ cal_dr6_pa6_f150 = gaussian 1.0 0.0024
 
 [cmb_foreground_params]
 ; TT
-a_tSZ = uniform 3.0 3.6
-a_kSZ = uniform 1.4 1.8
-a_p = uniform 6.2 7.6
-beta_p = uniform 1.8 2.2
-a_c = uniform 4.4 5.4
-beta_c = uniform 2.0 2.4
-a_s = uniform 2.8 3.4
-a_gtt = gaussian 2.79 0.45
+a_tSZ = uniform 0 inf
+a_kSZ = uniform 0 inf
+a_p = uniform 0 inf
+;beta_p = uniform 1.8 2.2
+a_c = uniform 0 inf
+;beta_c = uniform 2.0 2.4
+a_s = uniform 0 inf
+a_gtt = gaussian 7.95 0.32
 xi = uniform 0 0.2
-T_d = uniform 8.6 10.6
+;T_d = uniform 8.6 10.6
 
 ; TE
-a_gte = gaussian 0.36 0.04
-a_pste = uniform -1 1
+a_gte = gaussian 0.423 0.03
+;a_pste = uniform -1 1
 
 ; EE
-a_gee = gaussian 0.13 0.03
+a_gee = gaussian 0.1681 0.017
 a_psee = uniform 0 inf

examples/act-values.ini

@@ -28,42 +28,42 @@ w = -1.0   ;equation of state of dark energy
 wa = 0.0   ;equation of state of dark energy (redshift dependency)
 
 [act_params]
-bandint_shift_dr6_pa4_f220 = -14.4  0.0  14.4
-bandint_shift_dr6_pa5_f090 = -4.0  0.0  4.0
-bandint_shift_dr6_pa5_f150 = -5.2  0.0  5.2
-bandint_shift_dr6_pa6_f090 = -4.8  0.0  4.8
-bandint_shift_dr6_pa6_f150 = -4.4  0.0  4.4
-calG_all = 0.988  1.0  1.012
+bandint_shift_dr6_pa4_f220 = -25.0  0.0  25.0
+bandint_shift_dr6_pa5_f090 = -10.0  0.0  10.0
+bandint_shift_dr6_pa5_f150 = -10.0  0.0  10.0
+bandint_shift_dr6_pa6_f090 = -10.0  0.0  10.0
+bandint_shift_dr6_pa6_f150 = -10.0  0.0  10.0
+calG_all = 0.9  1.0  1.1
 calE_dr6_pa4_f220 = 1.0
 calE_dr6_pa5_f090 = 0.9  1.0  1.1
 calE_dr6_pa5_f150 = 0.9  1.0  1.1
 calE_dr6_pa6_f090 = 0.9  1.0  1.1
 calE_dr6_pa6_f150 = 0.9  1.0  1.1
-cal_dr6_pa4_f220 = 0.96  1.0  1.04
-cal_dr6_pa5_f090 = 0.96  1.0  1.04
-cal_dr6_pa5_f150 = 0.96  1.0  1.04
-cal_dr6_pa6_f090 = 0.96  1.0  1.04
-cal_dr6_pa6_f150 = 0.96  1.0  1.04
+cal_dr6_pa4_f220 = 0.9  1.0  1.1
+cal_dr6_pa5_f090 = 0.9  1.0  1.1
+cal_dr6_pa5_f150 = 0.9  1.0  1.1
+cal_dr6_pa6_f090 = 0.9  1.0  1.1
+cal_dr6_pa6_f150 = 0.9  1.0  1.1
 
 [cmb_foreground_params]
-a_tSZ = 3.0  3.3  3.6
-a_kSZ = 1.4  1.6  1.8
-a_p = 6.2  6.9  7.6
-beta_p = 1.8  2.0  2.2
-a_c = 4.4  4.9  5.4
-beta_c = 2.0  2.2  2.4
-a_s = 2.8  3.0999999999999996  3.4
-a_gtt = 0.99  2.79  4.59
-xi = 0  0.1  0.2
-T_d = 8.6  9.6  10.6
+a_tSZ = 0.0 	3.4		10.0
+a_kSZ = 0.0  	1.5  	10.0
+a_p = 	0.0  	7.7 	50.0
+beta_p = 0.0 	1.9 	5.0
+a_c = 	0.0 	3.6  	50.0
+beta_c = 0.0  	1.9  	5.0 ; lambda beta_p: beta_p
+a_s = 	0.0 	2.8  	50.0
+a_gtt = 0.0  	8.0  	50.0
+xi = 	0.0  	0.09  	0.2
+T_d = 	9.6
 alpha_dT = -0.6
 alpha_p = 1.0
-alpha_tSZ = 0.0
-a_gte = 0.19999999999999998  0.36  0.52
-a_pste = -1  0.0  1
+alpha_tSZ = -5.0 -0.6 	5.0
+a_gte = 0.0  	0.42  	1.0
+a_pste = -1.0  -0.026  	1.0
 alpha_dE = -0.4
-a_gee = 0.010000000000000009  0.13  0.25
-a_psee = 0  0.1  0.2
+a_gee = 0.0  0.17  1.0
+a_psee = 0  0.02  1.0
 alpha_dE = -0.4
 
 

examples/planck_lite.ini

@@ -22,6 +22,10 @@ use_low_ell_bins = T
 spectra = TTTEEE
 year = 2018
 
+; you can use these parameters to set the maximum ell for each spectrum:
+; ell_max_tt = 500
+; ell_max_te = 500
+; ell_max_ee = 500
 
 ; The consistency module translates between our chosen parameterization
 ; and any other that modules in the pipeline may want (e.g. camb)

likelihood/act-dr6-lite/act_dr6_lite_interface.py

@@ -1,7 +1,7 @@
 try:
     import act_dr6_cmbonly
 except ImportError:
-    raise RuntimeError('The act_dr6_cmbonly python module is required for the act_dr6_lite likelihood.')
+    raise RuntimeError('The act_dr6_cmbonly python module is required for the act_dr6_lite likelihood. It can be obtained from https://github.com/ACTCollaboration/DR6-ACT-lite')
 from cosmosis.datablock import names
 cosmo = names.cosmological_parameters
 import numpy as np
@@ -25,9 +25,9 @@ def execute(block, config):
     act = config
 
     cl_dict = {
-        "tt": block[names.cmb_cl, 'tt'],
-        "te": block[names.cmb_cl, 'te'],
-        "ee": block[names.cmb_cl, 'ee'],
+        "tt": np.append(np.array([0,0]), block[names.cmb_cl, 'tt']),
+        "te": np.append(np.array([0,0]), block[names.cmb_cl, 'te']),
+        "ee": np.append(np.array([0,0]), block[names.cmb_cl, 'ee']),
     }
 
     nuisance = {}