From 9c192b91fc30b6d5540bb0f428921b7e88a88220 Mon Sep 17 00:00:00 2001 From: Tilman Troester Date: Mon, 10 Jul 2023 20:28:47 +0200 Subject: [PATCH 1/5] don't run camb twice, reuse transfer function instead --- boltzmann/camb/camb_interface.py | 76 ++++++++++---------------------- 1 file changed, 24 insertions(+), 52 deletions(-) diff --git a/boltzmann/camb/camb_interface.py b/boltzmann/camb/camb_interface.py index 711bfe3e..491b4d12 100644 --- a/boltzmann/camb/camb_interface.py +++ b/boltzmann/camb/camb_interface.py @@ -18,6 +18,9 @@ MODE_ALL = "all" MODES = [MODE_BG, MODE_THERM, MODE_CMB, MODE_POWER, MODE_ALL] +DEFAULT_A_S = 2.1e-9 + + # See this table for description: #https://camb.readthedocs.io/en/latest/transfer_variables.html#transfer-variables matter_power_section_names = { @@ -301,14 +304,14 @@ def extract_camb_params(block, config, more_config): want_perturbations = more_config['mode'] not in [MODE_BG, MODE_THERM] want_thermal = more_config['mode'] != MODE_BG - # JMedit - check for input sigma8 - samplesig8 = block.has_value(cosmo, 'sigma_8_input') + if block.has_value(cosmo, 'sigma_8_input'): + if not block.has_value(cosmo, 'A_s'): + block[cosmo, 'A_s'] = DEFAULT_A_S # if want_perturbations: - if not samplesig8: #JMedit; if we're sampling sigma8, wait til we have A_s - init_power = extract_initial_power_params(block, config, more_config) + init_power = extract_initial_power_params(block, config, more_config) nonlinear = extract_nonlinear_params(block, config, more_config) -# else: + # else: # init_power = None # nonlinear = None @@ -338,8 +341,6 @@ def extract_camb_params(block, config, more_config): if (block.has_value(cosmo, "omega_nu") or block.has_value(cosmo, "omnuh2")) and not (block.has_value(cosmo, "mnu")): 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"): cosmology_params["H0"] = block[cosmo, "hubble"] @@ -348,13 +349,6 @@ def extract_camb_params(block, config, more_config): else: cosmology_params["cosmomc_theta"] = block[cosmo, "cosmomc_theta"]/100 - #JMedit - if samplesig8: - # compute linear matter power spec to figure out what A_s should be - # for desired input sigma8, add it to the datblock - sigma8_to_As(block, config, more_config, cosmology_params, dark_energy, reion) - init_power = extract_initial_power_params(block, config, more_config) - p = camb.CAMBparams( InitPower = init_power, Recomb = recomb, @@ -598,44 +592,19 @@ def save_cls(r, p, block): block[names.cmb_cl, "PE"] = cl[2:,2]*(ell*(ell+1))/(2*np.pi) -# JMedit: new function here -def sigma8_to_As(block, config, more_config, cosmology_params, dark_energy, reion): - """ - If input parameters include sigma_8_input, convert that to A_s. - - This function will run CAMB once to compute the linear matter power spectrum - - This function is adapted from the sigma8toAs module in the - KIDS KCAP repoistory written by by Tilman Troester. - """ - sigma_8_input = block[cosmo,'sigma_8_input'] - temp_As = 2.1e-9 - block[cosmo,'A_s'] = temp_As - init_power_temp = extract_initial_power_params(block, config, more_config) - - # do nothing except get linear power spectrum - p_temp = camb.CAMBparams(WantTransfer=True, - Want_CMB=False, Want_CMB_lensing=False, DoLensing=False, - NonLinear="NonLinear_none", - WantTensors=False, WantVectors=False, WantCls=False, - WantDerivedParameters=False, - want_zdrag=False, want_zstar=False,\ - DarkEnergy=dark_energy, - InitPower = init_power_temp,\ - ) - # making these choices match main setup - p_temp.set_accuracy(**more_config["accuracy_params"]) - p_temp.set_cosmology(ombh2 = block[cosmo, 'ombh2'], - omch2 = block[cosmo, 'omch2'], - omk = block[cosmo, 'omega_k'], - **more_config["cosmology_params"], - **cosmology_params) - p_temp.set_matter_power(redshifts=[0.], nonlinear=False, **more_config["transfer_params"]) - p_temp.Reion = reion - r_temp = camb.get_results(p_temp) - temp_sig8 = r_temp.get_sigma8()[-1] #what sigma8 comes out from using temp_As? - As = temp_As*(sigma_8_input/temp_sig8)**2 - block[cosmo,'A_s'] = As +def recompute_for_sigma8(camb_results, block, config, more_config): + # Recompute the power spectra using a new initial power spectrum + # This reuses the transfer function, which is what takes time to compute + sigma8 = camb_results.get_sigma8_0() + + sigma_8_target = block[cosmo,'sigma_8_input'] + block[cosmo,'A_s'] = block[cosmo,'A_s'] * sigma_8_target**2/sigma8**2 + init_power_scaled = extract_initial_power_params(block, config, more_config) + + camb_results.Params.set_initial_power(init_power_scaled) + camb_results.calc_power_spectra() + return camb_results + def execute(block, config): config, more_config = config @@ -649,6 +618,9 @@ def execute(block, config): else: # other modes r = camb.get_results(p) + if block.has_value(cosmo, 'sigma_8_input'): + r = recompute_for_sigma8(r, block, config, more_config) + except camb.CAMBError: if more_config["n_printed_errors"] <= more_config["max_printed_errors"]: print("CAMB error caught: for these parameters") From 0b9c088f03133758dd6252d45f1038609429ef1c Mon Sep 17 00:00:00 2001 From: Tilman Troester Date: Mon, 10 Jul 2023 21:27:04 +0200 Subject: [PATCH 2/5] fix subtle bug of get_BAO recalculating things --- boltzmann/camb/camb_interface.py | 38 ++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 17 deletions(-) diff --git a/boltzmann/camb/camb_interface.py b/boltzmann/camb/camb_interface.py index 491b4d12..cc871a93 100644 --- a/boltzmann/camb/camb_interface.py +++ b/boltzmann/camb/camb_interface.py @@ -394,9 +394,8 @@ def extract_camb_params(block, config, more_config): return p - - -def save_derived_parameters(r, p, block): +def save_derived_parameters(r, block): + p = r.Params # Write the default derived parameters to distance section derived = r.get_derived_params() for k, v in derived.items(): @@ -431,7 +430,8 @@ def save_derived_parameters(r, p, block): block[names.cosmological_parameters, cosmosis_name] = CAMB_value -def save_distances(r, p, block, more_config): +def save_distances(r, block, more_config): + p = r.Params # Evaluate z on a different grid than the spectra, so we can easily extend it further z_background = np.linspace( @@ -486,6 +486,7 @@ def save_distances(r, p, block, more_config): chistar = (r.conformal_time(0)- r.tau_maxvis) block[names.distances, "CHISTAR"] = chistar + def compute_growth_factor(r, block, P_tot, k, z, more_config): if P_tot is None: # If we don't have it already, get the default matter power interpolator, @@ -500,9 +501,8 @@ def compute_growth_factor(r, block, P_tot, k, z, more_config): return D - -def save_matter_power(r, p, block, more_config): - +def save_matter_power(r, block, more_config): + p = r.Params # Grids in k, z on which to save matter power. # There are two kmax values - the max one calculated directly, # and the max one extrapolated out too. We output to the larger @@ -520,8 +520,10 @@ def save_matter_power(r, p, block, more_config): # Get an interpolator. By default bicubic if size is large enough, # otherwise it drops down to linear. # First we do the linear version of the spectrum - P, zcalc, kcalc = r.get_matter_power_interpolator(nonlinear=False, var1=tt, var2=tt, return_z_k=True, - extrap_kmax=more_config['kmax_extrapolate']) + P, zcalc, kcalc = r.get_matter_power_interpolator( + nonlinear=False, var1=tt, var2=tt, return_z_k=True, + extrap_kmax=more_config['kmax_extrapolate'] + ) assert P.islog # P.P evaluates at k instead of logk p_k = P.P(z, k, grid=True) @@ -537,8 +539,10 @@ def save_matter_power(r, p, block, more_config): # Now if requested we also save the linear version if p.NonLinear is not camb.model.NonLinear_none: # Exact same process as before - P = r.get_matter_power_interpolator(nonlinear=True, var1=tt, var2=tt, - extrap_kmax=more_config['kmax_extrapolate']) + P = r.get_matter_power_interpolator( + nonlinear=True, var1=tt, var2=tt, + extrap_kmax=more_config['kmax_extrapolate'] + ) p_k = P.P(z, k, grid=True) section_name = matter_power_section_names[transfer_type] + "_nl" block.put_grid(section_name, "z", z, "k_h", k, "p_k", p_k) @@ -571,7 +575,7 @@ def save_matter_power(r, p, block, more_config): block[names.cosmological_parameters, "S_8"] = sigma_8[0]*np.sqrt(p.omegam/0.3) -def save_cls(r, p, block): +def save_cls(r, block): # Get total (scalar + tensor) lensed CMB Cls cl = r.get_total_cls(raw_cl=False, CMB_unit="muK") ell = np.arange(2,cl.shape[0]) @@ -581,7 +585,7 @@ def save_cls(r, p, block): block[names.cmb_cl, "BB"] = cl[2:,2] block[names.cmb_cl, "TE"] = cl[2:,3] - if p.DoLensing: + if r.Params.DoLensing: # Get CMB lensing potential # The cosmosis-standard-library clik interface expects ell(ell+1)/2 pi Cl # for all angular power spectra, including the lensing potential. @@ -631,14 +635,14 @@ def execute(block, config): more_config["n_printed_errors"] += 1 return 1 - save_derived_parameters(r, p, block) - save_distances(r, p, block, more_config) + save_derived_parameters(r, block) + save_distances(r, block, more_config) if p.WantTransfer: - save_matter_power(r, p, block, more_config) + save_matter_power(r, block, more_config) if p.WantCls: - save_cls(r, p, block) + save_cls(r, block) return 0 From 690be61e9a1d483465998f538c4726b67ead061a Mon Sep 17 00:00:00 2001 From: Tilman Troester Date: Mon, 10 Jul 2023 21:30:27 +0200 Subject: [PATCH 3/5] Allow sigma_8 as input --- boltzmann/camb/camb_interface.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/boltzmann/camb/camb_interface.py b/boltzmann/camb/camb_interface.py index cc871a93..f538da37 100644 --- a/boltzmann/camb/camb_interface.py +++ b/boltzmann/camb/camb_interface.py @@ -304,9 +304,9 @@ def extract_camb_params(block, config, more_config): want_perturbations = more_config['mode'] not in [MODE_BG, MODE_THERM] want_thermal = more_config['mode'] != MODE_BG - if block.has_value(cosmo, 'sigma_8_input'): - if not block.has_value(cosmo, 'A_s'): - block[cosmo, 'A_s'] = DEFAULT_A_S + # Set A_s for now, this gets rescaled later if sigma_8 is provided. + if not block.has_value(cosmo, 'A_s'): + block[cosmo, 'A_s'] = DEFAULT_A_S # if want_perturbations: init_power = extract_initial_power_params(block, config, more_config) @@ -601,7 +601,10 @@ def recompute_for_sigma8(camb_results, block, config, more_config): # This reuses the transfer function, which is what takes time to compute sigma8 = camb_results.get_sigma8_0() - sigma_8_target = block[cosmo,'sigma_8_input'] + if block.has_value(cosmo, 'sigma_8'): + sigma_8_target = block[cosmo,'sigma_8'] + else: + sigma_8_target = block[cosmo,'sigma_8_input'] block[cosmo,'A_s'] = block[cosmo,'A_s'] * sigma_8_target**2/sigma8**2 init_power_scaled = extract_initial_power_params(block, config, more_config) @@ -622,7 +625,8 @@ def execute(block, config): else: # other modes r = camb.get_results(p) - if block.has_value(cosmo, 'sigma_8_input'): + if (block.has_value(cosmo, 'sigma_8_input') + or block.has_value(cosmo, 'sigma_8')): r = recompute_for_sigma8(r, block, config, more_config) except camb.CAMBError: From a52ec3760670181a6efcf665c7cacff750dfd0f1 Mon Sep 17 00:00:00 2001 From: Tilman Troester Date: Wed, 12 Jul 2023 16:31:18 +0200 Subject: [PATCH 4/5] Checks for A_s and sigma_8, deprecation warning for sigma_8_input --- boltzmann/camb/camb_interface.py | 15 +++++++++------ boltzmann/camb/module.yaml | 5 +++++ examples/des-y1-values.ini | 2 +- 3 files changed, 15 insertions(+), 7 deletions(-) diff --git a/boltzmann/camb/camb_interface.py b/boltzmann/camb/camb_interface.py index f538da37..6af8cd94 100644 --- a/boltzmann/camb/camb_interface.py +++ b/boltzmann/camb/camb_interface.py @@ -304,6 +304,13 @@ def extract_camb_params(block, config, more_config): want_perturbations = more_config['mode'] not in [MODE_BG, MODE_THERM] want_thermal = more_config['mode'] != MODE_BG + if block.has_value('sigma_8_input'): + warnings.warn("Parameter sigma8_input will be deprecated in favour of sigma_8.") + block['sigma_8'] = block['sigma_8_input'] + + if block.has_value('A_s') and block.has_value('sigma_8'): + warnings.warn("Parameter A_s is being ignored in favour of sigma_8") + # Set A_s for now, this gets rescaled later if sigma_8 is provided. if not block.has_value(cosmo, 'A_s'): block[cosmo, 'A_s'] = DEFAULT_A_S @@ -601,10 +608,7 @@ def recompute_for_sigma8(camb_results, block, config, more_config): # This reuses the transfer function, which is what takes time to compute sigma8 = camb_results.get_sigma8_0() - if block.has_value(cosmo, 'sigma_8'): - sigma_8_target = block[cosmo,'sigma_8'] - else: - sigma_8_target = block[cosmo,'sigma_8_input'] + sigma_8_target = block[cosmo,'sigma_8'] block[cosmo,'A_s'] = block[cosmo,'A_s'] * sigma_8_target**2/sigma8**2 init_power_scaled = extract_initial_power_params(block, config, more_config) @@ -625,8 +629,7 @@ def execute(block, config): else: # other modes r = camb.get_results(p) - if (block.has_value(cosmo, 'sigma_8_input') - or block.has_value(cosmo, 'sigma_8')): + if block.has_value(cosmo, 'sigma_8'): r = recompute_for_sigma8(r, block, config, more_config) except camb.CAMBError: diff --git a/boltzmann/camb/module.yaml b/boltzmann/camb/module.yaml index 46091e28..d3689659 100644 --- a/boltzmann/camb/module.yaml +++ b/boltzmann/camb/module.yaml @@ -264,6 +264,11 @@ inputs: meaning: Primordial scalar spectral amplitude type: real default: + sigma_8: + meaning: Amplitude of linear perturbations at z=0 on scales of 8 Mpc/h. + Only A_s or sigma_8 should be provided. + type: real + default: hubble: meaning: Hubble parameter in km/s/Mpc type: real diff --git a/examples/des-y1-values.ini b/examples/des-y1-values.ini index e8a6a704..96dba977 100644 --- a/examples/des-y1-values.ini +++ b/examples/des-y1-values.ini @@ -7,7 +7,7 @@ n_s = 0.87 0.9676 1.07 ; the likelihood will be the same, but the prior will ; be different. ; A_s = 0.5e-09 2.260574e-09 5.0e-09 -sigma_8_input = 0.5 0.8341371714072584 1.0 +sigma_8 = 0.5 0.8341371714072584 1.0 w = -1.0 omega_k = 0.0 tau = 0.08 From ee6df23dafd1209a0de72b641d87c5c04260b644 Mon Sep 17 00:00:00 2001 From: Tilman Troester Date: Wed, 12 Jul 2023 16:58:36 +0200 Subject: [PATCH 5/5] should have tested locally... --- boltzmann/camb/camb_interface.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/boltzmann/camb/camb_interface.py b/boltzmann/camb/camb_interface.py index 6af8cd94..f4ff7d8d 100644 --- a/boltzmann/camb/camb_interface.py +++ b/boltzmann/camb/camb_interface.py @@ -304,11 +304,11 @@ def extract_camb_params(block, config, more_config): want_perturbations = more_config['mode'] not in [MODE_BG, MODE_THERM] want_thermal = more_config['mode'] != MODE_BG - if block.has_value('sigma_8_input'): + if block.has_value(cosmo, 'sigma_8_input'): warnings.warn("Parameter sigma8_input will be deprecated in favour of sigma_8.") - block['sigma_8'] = block['sigma_8_input'] + block[cosmo, 'sigma_8'] = block[cosmo, 'sigma_8_input'] - if block.has_value('A_s') and block.has_value('sigma_8'): + if block.has_value(cosmo, 'A_s') and block.has_value(cosmo, 'sigma_8'): warnings.warn("Parameter A_s is being ignored in favour of sigma_8") # Set A_s for now, this gets rescaled later if sigma_8 is provided.