Skip to content

Commit

Permalink
Merge branch 'main' into two-point-one-point
Browse files Browse the repository at this point in the history
  • Loading branch information
andrejdvornik committed Sep 19, 2024
2 parents 0af6bc9 + 6096352 commit 68089c8
Show file tree
Hide file tree
Showing 32 changed files with 1,295 additions and 29 deletions.
22 changes: 15 additions & 7 deletions .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ jobs:

- name: Install dependencies with conda
shell: bash -l {0}
run: mamba install -c conda-forge "cosmosis>=3.2" cosmosis-build-standard-library pytest
run: mamba install -c conda-forge "cosmosis==3.9.2" cosmosis-build-standard-library pytest

- name: Get Cached Planck
uses: actions/cache/restore@v3
Expand Down Expand Up @@ -137,8 +137,14 @@ jobs:
- name: Install likelihood python dependencies from pip
shell: bash -l {0}
run: |
pip install "act_dr6_lenslike>=1.0.2"
pip install "act_dr6_lenslike>=1.0.2" "git+https://github.com/carronj/planck_PR4_lensing"
# Not working - don't know why
# - name: Install recent python dependencies
# if: ${{ (matrix.pyversion != '3.7') && (matrix.pyversion != '3.8') }}
# run: |
# pip install "candl-like"

- name: Run Tests
shell: bash -l {0}
run: |
Expand Down Expand Up @@ -230,13 +236,15 @@ jobs:
- name: Install python dependencies
run: |
python -m pip install --upgrade pip wheel setuptools
pip install "cosmosis>=3.2" "nautilus-sampler==0.6.*"
pip install "cosmosis==3.9.2" "nautilus-sampler==1.0.1" "scipy<1.14"
pip install -v --no-cache-dir --no-binary=mpi4py,camb mpi4py camb
pip install fitsio astropy fast-pt "Cython<3.0" jupyter sacc
- name: Install likelihood python dependencies
pip install fitsio astropy fast-pt "Cython>=3.0" jupyter sacc "act_dr6_lenslike>=1.0.2" "git+https://github.com/carronj/planck_PR4_lensing"
- name: Install recent python dependencies
if: ${{ (matrix.python-version != '3.7') && (matrix.python-version != '3.8') }}
run: |
pip install "act_dr6_lenslike>=1.0.2"
pip install "candl-like"
- name: Build
run: |
Expand Down
4 changes: 4 additions & 0 deletions boltzmann/camb/camb_interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@

DEFAULT_A_S = 2.1e-9

C_KMS = 299792.458

# See this table for description:
#https://camb.readthedocs.io/en/latest/transfer_variables.html#transfer-variables
Expand Down Expand Up @@ -430,6 +431,9 @@ def save_derived_parameters(r, block):
for k, v in derived.items():
block[names.distances, k] = v
block[names.distances, 'rs_zdrag'] = block[names.distances, 'rdrag']
zstar = derived['zstar']
shift = r.angular_diameter_distance(zstar) * (1 + zstar) * (p.omegam * p.H0**2)**0.5 / C_KMS
block[names.distances, "cmbshift"] = shift

p.omegal = 1 - p.omegam - p.omk
p.ommh2 = p.omegam * p.h**2
Expand Down
16 changes: 15 additions & 1 deletion boltzmann/camb/module.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -464,11 +464,14 @@ outputs:
d_l:
meaning: Luminosity distance in Mpc
type: real 1d
d_v:
meaning: Angular average of comoving angular diameter and line of sight distance in Mpc
type: real 1d
mu:
meaning: Distance modulus
type: real 1d
h:
meaning: Hubble parameter with in units of Mpc
meaning: Hubble parameter as function of redshift in units of Mpc
type: real 1d
age:
meaning: Age of universe in GYr
Expand Down Expand Up @@ -513,6 +516,9 @@ outputs:
thetarseq:
meaning: Angle 100 r_s(eq)/DA(zstar)
type: real
cmbshift:
meaning: CMB shift parameter equation 69 of Komatsu et al 2009
type: real
growth_parameters:
z:
meaning: Redshift samples of other values in this section, (all if mode=power
Expand Down Expand Up @@ -602,3 +608,11 @@ outputs:
p_k:
meaning: Non-linear power spectrum at samples in (Mpc/h)^-3.
type: real 2d
transfer_func:
k_h:
meaning: Wavenumbers k of samples in Mpc/h.
Different than the k of power spectra!
type: real 1d
t_k:
meaning: Linear transfer function at z=0.
type: real 1d
9 changes: 8 additions & 1 deletion boltzmann/class/class_interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -202,7 +202,8 @@ def get_class_outputs(block, c, config):
block[distances, 'nz'] = nz

# Save distance samples
block[distances, 'd_l'] = np.array([c.luminosity_distance(zi) for zi in z])
d_l = np.array([c.luminosity_distance(zi) for zi in z])
block[distances, 'd_l'] = d_l
d_a = np.array([c.angular_distance(zi) for zi in z])
block[distances, 'd_a'] = d_a
block[distances, 'd_m'] = d_a * (1 + z)
Expand All @@ -215,6 +216,12 @@ def get_class_outputs(block, c, config):
h_z = np.array([c.Hubble(zi) for zi in z])
block[distances, 'H'] = h_z

mu = np.zeros(z.size)
mu[0] = -np.inf
mu[1:] = 5.0 * np.log10(d_l[1:]) + 25.0
block[distances, 'mu'] = mu


##
# Now the CMB C_ell
##
Expand Down
12 changes: 6 additions & 6 deletions boltzmann/class/class_v3.2.0/python/classy.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ def viewdictitems(d):
return d.viewitems()

ctypedef np.float_t DTYPE_t
ctypedef np.int_t DTYPE_i
# ctypedef np.int_t DTYPE_i



Expand Down Expand Up @@ -100,10 +100,10 @@ cdef class Class:
cdef distortions sd
cdef file_content fc

cpdef int computed # Flag to see if classy has already computed with the given pars
cpdef int allocated # Flag to see if classy structs are allocated already
cpdef object _pars # Dictionary of the parameters
cpdef object ncp # Keeps track of the structures initialized, in view of cleaning.
cdef int computed # Flag to see if classy has already computed with the given pars
cdef int allocated # Flag to see if classy structs are allocated already
cdef object _pars # Dictionary of the parameters
cdef object ncp # Keeps track of the structures initialized, in view of cleaning.

# Defining two new properties to recover, respectively, the parameters used
# or the age (set after computation). Follow this syntax if you want to
Expand All @@ -128,7 +128,7 @@ cdef class Class:
self.set(**_pars)

def __cinit__(self, default=False):
cpdef char* dumc
cdef char* dumc
self.allocated = False
self.computed = False
self._pars = {}
Expand Down
64 changes: 64 additions & 0 deletions examples/candl_test.ini
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
[runtime]
sampler = test
root = ${PWD}
resume = T
verbosity = noisy

[test]
save_dir = output/spt3g_2018_lensing
fatal_errors = T

[nautilus]
n_live = 1500
verbose = T

[emcee]
walkers = 64
samples = 1000
nsteps = 5


[pipeline]
; these names refer to sections later in the file:
modules = consistency camb spt3g_2018_lensing
values = examples/candl_test_values.ini
priors = examples/candl_test_priors.ini
debug = F
timing = F
extra_output = cosmological_parameters/sigma_8 cosmological_parameters/omega_m

[spt3g_2018_lensing]
file = ./likelihood/candl/candl_cosmosis_interface.py ; Location of interface code - change depending on the location of your .ini file
data_set = 'candl.data.SPT3G_2018_Lens' ; Data set or path to .yaml file
variant = 'use_CMB' ; Select a variant of the data set if pointing to an index file
lensing = T ; Switch on for lensing likelihoods
feedback = T ; Switch on to request feedback from candl initialisation
data_selection = "..." ; Select a subset of the data set
clear_1d_internal_priors = F ; Switch off to use candl internal 1d priors
clear_nd_internal_priors = F ; Switch on to ignore candl internal higher dimensional priors. Careful: higher-dimensional priors are not implemented in CosmoSIS itself.
force_ignore_transformations = '' ; Backdoor if you want to ignore certain transformations in the data model.

[output]
filename = output/spt3g-18-lens.txt

; The consistency module translates between our chosen parameterization
; 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
feedback = 0 ; verbosity of output
mode = cmb ; mode to run camb. For CMB lensing only, cmb is sufficient
lmax = 4000 ; max ell to use for cmb calculation
lens_margin = 1250 ; Lmax
AccuracyBoost = 1.0 ; CAMB accuracy boost parameter
lSampleBoost = 1.0 ; CAMB lsample boost parameter
lAccuracyBoost = 1.0 ; CAMB lAccuracy boost parameter
lens_potential_accuracy = 4 ; CAMB lens_potential accuracy paramater
do_tensors = T ;include tensor modes
do_lensing = T ;lensing is required w/ Planck data
NonLinear = lens ; Non-linear calculation
theta_H0_range = "20 100" ; Set bounds in H0
halofit_version = takahashi ; Halofit version
3 changes: 3 additions & 0 deletions examples/candl_test_priors.ini
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
[cosmological_parameters]
ombh2 = normal 0.02233 0.00036
n_s = normal 0.96 0.02
16 changes: 16 additions & 0 deletions examples/candl_test_values.ini
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@

[cosmological_parameters]
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.9 1.04 1.20
omega_k = 0.0
mnu = 0.06
nnu = 3.046
yhe = 0.245341
tau = 0.0543
num_massive_neutrinos = 3

[nuisance_parameters]
A_fg = 0.0 1.0 2.0
60 changes: 60 additions & 0 deletions examples/desi-figure-2-campaign.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
# This reproduces the chains that make up figure 2 of the DESI paper
# 2404.03002

output_dir: "./output/desi-campaign"

runs:
- name: base
base: examples/desi.ini
params:
- sampler = metropolis
- metropolis.samples = 20000
- metropolis.use_cobaya = T
- metropolis.tuning_frequency = 100
- metropolis.tuning_grace = 500
- metropolis.tuning_end = 2000
- emcee.walkers = 32
- emcee.samples = 1000
- emcee.nsteps = 20

- name: BGS
parent: base
params:
- desi.desi_data_sets = BGS
- metropolis.samples = 30000

- name: LRG1
parent: base
params:
- desi.desi_data_sets = LRG1

- name: LRG2
parent: base
params:
- desi.desi_data_sets = LRG2

- name: LRG3+ELG1
parent: base
params:
- desi.desi_data_sets = LRG3+ELG1

- name: ELG2
parent: base
params:
- desi.desi_data_sets = ELG2

- name: QSO
parent: base
params:
- desi.desi_data_sets = QSO

- name: Lya QSO
parent: base
params:
- desi.desi_data_sets = Lya QSO

- name: All
parent: base
params:
# - pipeline.timing=T
- desi.desi_data_sets = BGS,LRG1,LRG2,LRG3+ELG1,ELG2,QSO,Lya QSO
89 changes: 89 additions & 0 deletions examples/desi-figure-2.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,89 @@
import cosmosis.campaign
from cosmosis.postprocessing.v2 import Chain
import getdist.plots
import matplotlib.pyplot as plt
plt.rcParams['text.usetex'] = True

# The campaign file defines all the chains that we want to generate
campaign, _ = cosmosis.campaign.parse_yaml_run_file("./examples/desi-figure-2-campaign.yaml")


# First make all the chains, by looping through the campaign.
# The next major version of CosmoSIS will do this for you more easily.
for name, run in campaign.items():
if name != "base":
cosmosis.campaign.launch_run(run)

# These colors match those used in the DESI paper as
# best I could
colors = {
'All': [(0.0, 0.0, 0.0), (0.0, 0.0, 0.0)],
'BGS': [(0.729, 0.811, 0.552), (0.517, 0.745, 0.012)],
'ELG2': [(0.584, 0.725, 0.772), (0.219, 0.447, 0.6)],
'LRG1': [(0.996, 0.854, 0.596), (0.956, 0.627, 0.086)],
'LRG2': [(0.960, 0.678, 0.580), (0.929, 0.262, 0.094)],
'LRG3+ELG1': [(0.733, 0.603, 0.741), (0.341, 0.211, 0.701)],
'Lya QSO': [(0.733, 0.556, 0.749), (0.439, 0.0, 0.474)],
'QSO': [(0.662, 0.815, 0.725), (0.247, 0.552, 0.360)]
}


# This is the order we want the contours to be plotted in. It's not the same
# as the legend order. This is because we don't want the big contours to
# cover up the small ones
order = ["BGS", "LRG1", "QSO", "ELG2", "Lya QSO", "LRG3+ELG1", "LRG2", "All"]


# The chains are complete now so load them all in
chain_data = {}
for i, (name, run) in enumerate(campaign.items()):
# The base parent run it not used, it's just a template.
# We will soon be able to specify this in the yaml file
if name == "base" :
continue

# Load the chain using a chain object. This is in the upcoming
# v2 api for cosmosis postprocessing, which will soon replace the
# main version. It uses GetDist to do everything.
chain_file = campaign[name]['params']['output', 'filename']
chain = Chain.load(chain_file, name=name, burn=0.2)

# Store the GetDist MCSamples object in a dictionary
chain_data[name] = chain.mcsamples


# split the chains into two groups, the single-data set chains
# and the combined one. This is because the contours are all filled
# for the single data sets, but not for the combined one.
combined_sample = chain_data["All"]
single_data_samples = {name: chain_data[name] for name in order if name != "All"}

# This is the order we want the labels to appear in.
new_order = [order.index(name) for name in colors.keys()]

# The rest of this is standard GetDist calls.
plotter = getdist.plots.get_single_plotter()


plotter.plot_2d(list(single_data_samples.values()),
"DISTANCES--H0RD",
"cosmological_parameters--omega_m",
filled=True,
colors=[colors[name] for name in single_data_samples.keys()],
add_legend_proxy=True,
lims=[70, 130, 0.1, 0.7],)

plotter.plot_2d(combined_sample,
"DISTANCES--H0RD",
"cosmological_parameters--omega_m",
filled=False,
colors=["k", "k"],
lims=[70, 130, 0.1, 0.7],
add_legend_proxy=True)


ax=plotter.get_axes()
plotter.add_legend(list(single_data_samples.keys()) + ["All"], label_order=new_order, legend_loc="upper right")
ax.set_xlabel(r"$\mathrm{H}_0 r_d \, [100 \mathrm{km}\, \mathrm{s}^{-1}]$")
ax.set_ylabel(r"$\Omega_m$")
plotter.export("output/desi.pdf")
Loading

0 comments on commit 68089c8

Please sign in to comment.