Skip to content

Commit

Permalink
Merge pull request #86 from tilmantroester/better_sample_sigma_8
Browse files Browse the repository at this point in the history 
Better sampling of sigma_8
  • Loading branch information
@joezuntz
joezuntz authored Jul 17, 2023
2 parents 8ae31eb + ee6df23 commit 54f641f
Show file tree
Hide file tree
Showing 3 changed files with 58 additions and 70 deletions.
121 changes: 52 additions & 69 deletions boltzmann/camb/camb_interface.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 = {
Expand Down Expand Up @@ -301,14 +304,21 @@ 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'):
warnings.warn("Parameter sigma8_input will be deprecated in favour of sigma_8.")
block[cosmo, 'sigma_8'] = block[cosmo, 'sigma_8_input']

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.
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

Expand Down Expand Up @@ -338,8 +348,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"]
Expand All @@ -348,13 +356,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,
Expand Down Expand Up @@ -400,9 +401,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():
Expand Down Expand Up @@ -437,7 +437,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(
Expand Down Expand Up @@ -492,6 +493,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,
Expand All @@ -506,9 +508,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
Expand All @@ -526,8 +527,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)
Expand All @@ -543,8 +546,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)
Expand Down Expand Up @@ -577,7 +582,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])
Expand All @@ -587,7 +592,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.
Expand All @@ -598,44 +603,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']
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
Expand All @@ -649,6 +629,9 @@ def execute(block, config):
else:
# other modes
r = camb.get_results(p)
if block.has_value(cosmo, 'sigma_8'):
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")
Expand All @@ -659,14 +642,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

Expand Down
5 changes: 5 additions & 0 deletions boltzmann/camb/module.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion examples/des-y1-values.ini
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 54f641f

Please sign in to comment.