33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407 | def get_theory_cls(
theta_dict,
spectra,
path2wsp,
ellmin=0,
ellmax=6500,
binning=True,
precise=False,
extra_class_args={},
path_desi=None,
return_sigma8=False,
matter_power_spectrum="halofit",
verbose=False,
):
"""
Compute the theory Cls for the given cosmological parameters.
:param theta_dict: Dictionary containing the cosmological parameters
:param spectra: List of spectra to calculate
:param path2wsp: Path to the NaMaster workspace for binning
:param ellmin: Minimum ell value
:param ellmax: Maximum ell value
:param binning: Whether to apply binning to the Cls
:param extra_class_args: Extra arguments passed to Class if desired
"""
path_boss = utils.rel2realpath(
"../../combined_probes_data/survey_data/data_boss_dr12/"
)
path_kids = utils.rel2realpath(
"../../combined_probes_data/survey_data/data_kids_1000/"
)
# Define default values for the parameters
defaults = {
# "h": 0.67,
"omega_b": 0.022,
"omega_cdm": 0.12,
# "sigma8": 0.8,
"n_s": 0.96,
"tau_reio": 0.055,
"m_ncdm": 0.02,
"del_z1": 0.0,
"del_z2": 0.0,
"del_z3": 0.0,
"del_z4": 0.0,
"del_z5": 0.0,
"bl": 1.837,
"bc": 2.086,
"A_IA": 0.0,
"A_L": 1.0,
"bias_d1": 1.8,
"bias_d2": 2.0,
"bias_d3": 2.2,
"bias_d4": 2.4,
"w0_fld": -1.0,
"wa_fld": 0.0,
# "fEDE": 0.0,
# "thetai_scf": 2.83,
# "log10z_c": 3.562,
"fraction_axion_ac": 0.0,
"log10_axion_ac": -3.56,
"theta_i": 2.72,
}
for param in defaults.keys():
if param not in theta_dict:
if verbose:
LOGGER.warning(
f"{param} not found in theta_dict, "
f"using default value: {defaults[param]}"
)
theta_dict[param] = theta_dict.get(param, defaults[param])
if "logT_AGN" in theta_dict:
if verbose:
LOGGER.info("Using HMcode Baryon feedback model with KiDS-1000 TT setup")
extra_parameters = {
"camb": {
"halofit_version": "mead2020_feedback",
"HMCode_logT_AGN": theta_dict["logT_AGN"],
"dark_energy_model": "ppf",
},
}
else:
extra_parameters = None
cl_dict = {}
ells = np.arange(ellmin, ellmax)
class_computed = False
lss_computed = False
for spec in spectra:
if spec in np.array(["tt", "te", "ee"]):
if not class_computed:
# precision parameters from https://arxiv.org/abs/2503.14454
precision_parameters = {
"modes": "s",
"delta_l_max": 1800,
"l_logstep": 1.025,
"l_linstep": 20,
"perturbations_sampling_stepsize": 0.05,
"l_switch_limber": 30.0,
"hyper_sampling_flat": 32,
"l_max_g": 40,
"l_max_ur": 35,
"l_max_pol_g": 60,
"ur_fluid_approximation": 2,
"ur_fluid_trigger_tau_over_tau_k": 130.0,
"radiation_streaming_approximation": 2,
"radiation_streaming_trigger_tau_over_tau_k": 240.0,
"hyper_flat_approximation_nu": 7000.0,
"transfer_neglect_delta_k_S_t0": 0.17,
"transfer_neglect_delta_k_S_t1": 0.05,
"transfer_neglect_delta_k_S_t2": 0.17,
"transfer_neglect_delta_k_S_e": 0.17,
"accurate_lensing": 1,
"start_small_k_at_tau_c_over_tau_h": 0.0004,
"start_large_k_at_tau_h_over_tau_k": 0.05,
"tight_coupling_trigger_tau_c_over_tau_h": 0.005,
"tight_coupling_trigger_tau_c_over_tau_k": 0.008,
"start_sources_at_tau_c_over_tau_h": 0.006,
"l_max_ncdm": 30,
"tol_ncdm_synchronous": 1e-6,
}
classy_params = {
# "h": theta_dict["h"],
"omega_b": theta_dict["omega_b"],
"omega_cdm": theta_dict["omega_cdm"],
# "sigma8": theta_dict["sigma8"],
"n_s": theta_dict["n_s"],
"tau_reio": theta_dict["tau_reio"],
"m_ncdm": theta_dict["m_ncdm"],
"A_L": theta_dict["A_L"],
"w0_fld": theta_dict["w0_fld"],
"wa_fld": theta_dict["wa_fld"],
"Omega_Lambda": 0.0, # required for DE model
"N_ur": 0.00641,
"N_ncdm": 1,
"deg_ncdm": 3,
"lensing": "yes",
"output": "tCl,pCl,lCl",
"l_max_scalars": ellmax,
}
M = Class()
if theta_dict.get("fraction_axion_ac", 0.0) > 0.0:
# --- build AxiCLASS block ---
classy_params.update(
{
"scf_potential": "axion",
"n_axion": 3,
"log10_axion_ac": theta_dict["log10_axion_ac"],
"fraction_axion_ac": theta_dict["fraction_axion_ac"],
"scf_parameters": f"{theta_dict['theta_i']},0.0",
"scf_evolve_as_fluid": False,
"scf_evolve_like_axionCAMB": False,
"scf_has_perturbations": True,
"attractor_ic_scf": False,
"compute_phase_shift": False,
"include_scf_in_delta_m": True,
"include_scf_in_delta_cb": True,
"loop_over_background_for_closure_relation": "yes",
"do_shooting": "yes",
"do_shooting_scf": "yes",
}
)
# remove w0_fld and wa_fld and Omega_Lambda being set to 0.0
classy_params.pop("Omega_Lambda")
classy_params.pop("w0_fld")
classy_params.pop("wa_fld")
M = axiClass()
# switch to modrec set up (check if not using axiclassy)
elif any(f"q_{i}" in theta_dict for i in range(1, 8)):
# Check if any modrec parameters are provided
LOGGER.info("ModRec model detected")
# Set up modrec parameters according to CLASS requirements
modrec_model = {"zmin": 533.333333, "zmax": 1600, "nfree": 7}
# Create pivot points (evenly spaced between zmin and zmax)
pivots = np.linspace(
modrec_model["zmin"],
modrec_model["zmax"],
modrec_model["nfree"] + 2,
)
pivots_str = ",".join(f"{p:.4f}" for p in pivots)
# Collect control points (q_1 through q_7)
list_of_control_points = []
for i in range(1, 8):
q_key = f"q_{i}"
if q_key in theta_dict:
list_of_control_points.append(theta_dict[q_key])
else:
list_of_control_points.append(0.0)
# Create control points string: 0.0, q_1, q_2, ..., q_7, 0.0
cp_list = [0.0] + list_of_control_points + [0.0]
cp_str = ",".join(f"{c:.4e}" for c in cp_list)
# Add modrec recombination settings
recombination_settings = {
"xe_pert_type": "control",
"xe_pert_num": modrec_model["nfree"] + 2,
"xe_control_pivots": pivots_str,
"zmin_pert": modrec_model["zmin"],
"zmax_pert": modrec_model["zmax"],
"xe_interp_type": "cubic",
"xe_control_is_bounded": "yes",
"xe_control_points": cp_str,
}
classy_params.update(recombination_settings)
# Keep w0_fld and wa_fld for modrec model
if return_sigma8:
classy_params["output"] = "tCl,pCl,lCl,mPk"
if "sigma8" in theta_dict.keys():
classy_params["sigma8"] = theta_dict["sigma8"]
elif "ln10^{10}A_s" in theta_dict.keys():
classy_params["ln10^{10}A_s"] = theta_dict["ln10^{10}A_s"]
else:
raise ValueError(
"sigma8 or ln10^{10}A_s must be provided in theta_dict"
)
# if both ln10^{10}A_s and sigma8 are provided, use ln10^{10}A_s
if (
"ln10^{10}A_s" in theta_dict.keys()
and "sigma8" in theta_dict.keys()
):
classy_params.pop("sigma8")
if "h" in theta_dict.keys():
classy_params["h"] = theta_dict["h"]
elif "theta_s_100" in theta_dict.keys():
classy_params["theta_s_100"] = theta_dict["theta_s_100"]
else:
raise ValueError("h or theta_s_100 must be provided in theta_dict")
M.set(classy_params)
if precise:
M.set(precision_parameters)
M.compute()
class_computed = True
if return_sigma8:
sigma8 = M.sigma8()
cl_lensed_class = M.lensed_cl()
Tcmb0 = M.T_cmb()
if spec == "tt":
cl_dict["tt"] = cl_lensed_class["tt"][ellmin:ellmax] * Tcmb0**2
elif spec == "ee":
cl_dict["ee"] = cl_lensed_class["ee"][ellmin:ellmax] * Tcmb0**2
elif spec == "te":
cl_dict["te"] = cl_lensed_class["te"][ellmin:ellmax] * Tcmb0**2
# switch to kk being computed by CLASS for the auto spectrum!!
# elif spec == "kk":
# ells = cl_lensed_class["ell"]
# cl_dict["kk"] = (
# cl_lensed_class["pp"] * (1.0 / 4.0)
# * ells**2.0 * (ells + 1.0) ** 2.0
# )
else:
if not lss_computed:
# switch to axiclassy set up
if theta_dict.get("fraction_axion_ac", 0.0) > 0.0:
cosmo_ccl = ccl.Cosmology(
Omega_c=theta_dict["omega_cdm"] / theta_dict["h"] ** 2,
Omega_b=theta_dict["omega_b"] / theta_dict["h"] ** 2,
h=theta_dict["h"],
sigma8=theta_dict["sigma8"],
n_s=theta_dict["n_s"],
m_nu=theta_dict["m_ncdm"]
* 3, # bugfix! This needs to be the total neutrino mass!
mass_split="equal", # this is likely my issue!!--
# does this split the mass provided
Neff=3.044,
transfer_function="boltzmann_class",
extra_parameters=None,
baryonic_effects=None,
matter_power_spectrum=matter_power_spectrum,
fraction_axion_ac=theta_dict["fraction_axion_ac"],
log10_axion_ac=theta_dict["log10_axion_ac"],
scf_parameters=f"{theta_dict['theta_i']},0.0",
)
# switch to modrec set up
elif any(f"q_{i}" in theta_dict for i in range(1, 8)):
# Check if any modrec parameters are provided
LOGGER.info("ModRec model detected for LSS")
cosmo_ccl = ccl.Cosmology(
Omega_c=theta_dict["omega_cdm"] / theta_dict["h"] ** 2,
Omega_b=theta_dict["omega_b"] / theta_dict["h"] ** 2,
h=theta_dict["h"],
sigma8=theta_dict["sigma8"],
n_s=theta_dict["n_s"],
m_nu=theta_dict["m_ncdm"]
* 3, # bugfix! This needs to be the total neutrino mass!
mass_split="equal", # this is likely my issue!!--
# does this split the mass provided
Neff=3.044,
w0=theta_dict.get("w0_fld", -1.0),
wa=theta_dict.get("wa_fld", 0.0),
transfer_function="boltzmann_class",
extra_parameters=None,
baryonic_effects=None,
matter_power_spectrum=matter_power_spectrum,
q_1=theta_dict.get("q_1"),
q_2=theta_dict.get("q_2"),
q_3=theta_dict.get("q_3"),
q_4=theta_dict.get("q_4"),
q_5=theta_dict.get("q_5"),
q_6=theta_dict.get("q_6"),
q_7=theta_dict.get("q_7"),
)
else:
cosmo_ccl = ccl.Cosmology(
Omega_c=theta_dict["omega_cdm"] / theta_dict["h"] ** 2,
Omega_b=theta_dict["omega_b"] / theta_dict["h"] ** 2,
h=theta_dict["h"],
sigma8=theta_dict["sigma8"],
n_s=theta_dict["n_s"],
m_nu=theta_dict["m_ncdm"]
* 3, # bugfix! This needs to be the total neutrino mass!
mass_split="equal", # this is likely my issue!!--
# does this split the mass provided
Neff=3.044,
w0=theta_dict["w0_fld"],
wa=theta_dict["wa_fld"],
transfer_function="boltzmann_camb",
extra_parameters=extra_parameters,
baryonic_effects=None,
matter_power_spectrum=matter_power_spectrum,
)
cl_dict = theory_utils.get_lss_spectra(
cosmo_ccl,
theta_dict,
spectra,
cl_dict,
ells,
path_boss,
path_kids,
path_desi,
)
lss_computed = True
if binning:
cl_binned_dict = {
spec: theory_utils.apply_binning(spec, path2wsp, cl_dict)
for spec in spectra
}
if return_sigma8:
return cl_binned_dict, cl_dict, sigma8
else:
return cl_binned_dict, cl_dict
else:
return cl_dict
|