Merge branch 'joezuntz:main' into main

.github/workflows/ci.yml

@@ -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
@@ -137,7 +137,7 @@ 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"
 
       - name: Run Tests
         shell: bash -l {0}
@@ -158,7 +158,7 @@ jobs:
     needs: Download_Data
     strategy:
       matrix:
-        python-version: [3.7, 3.8, 3.9, "3.10", "3.11"]
+        python-version: [3.7, 3.8, 3.9, "3.10", "3.11", "3.12"]
 
     steps:
     - name: Checkout repository
@@ -230,9 +230,9 @@ 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
+        pip install fitsio astropy fast-pt "Cython>=3.0" jupyter sacc "git+https://github.com/carronj/planck_PR4_lensing"
 
     - name: Install likelihood python dependencies
       run: |

boltzmann/camb/camb_interface.py

@@ -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
@@ -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
@@ -565,6 +569,16 @@ def save_matter_power(r, block, more_config):
         # Save the linear
         section_name = matter_power_section_names[transfer_type] + "_lin"
         block.put_grid(section_name, "z", z, "k_h", k, "p_k", p_k)
+
+        # Save the transfer function
+        # 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 ...
+        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")
+        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)
 
         # Now if requested we also save the linear version
         if p.NonLinear is not camb.model.NonLinear_none:

boltzmann/camb/module.yaml

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

boltzmann/class/class_v3.2.0/python/classy.pyx

@@ -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
 
 
 
@@ -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
@@ -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 = {}

examples/desi-figure-2-campaign.yaml

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

examples/desi-figure-2.py

@@ -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")

examples/desi-values.ini

@@ -0,0 +1,25 @@
+[cosmological_parameters]
+; Fix these to fiducial values
+h0 = 0.7
+ombh2 = 0.02236
+omega_m      =  0.01 0.3  0.99
+rdh          = 10.0  100.0  1000.0
+
+
+; hubble       =  20.  70.  120.
+; ombh2        =  0.005  0.02  0.1
+
+; These are needed to run camb but will have no
+; effect on the observables or likelihood
+n_s          =  0.97
+A_s          =  2.19e-9
+w            = -1.0
+wa           =  0.0
+
+; New parametrization and names:
+mnu = 0.0
+num_massive_neutrinos = 0
+nnu = 3.046
+
+omega_k      =  0.0
+tau          =  0.0697186

examples/desi.ini

@@ -0,0 +1,42 @@
+[runtime]
+sampler = test
+verbosity = quiet
+resume = T
+
+[pipeline]
+modules =  consistency  camb sample_rdh  desi
+timing=F
+extra_output = distances/rs_zdrag distances/h0rd
+values = examples/desi-values.ini
+
+[maxlike]
+max_posterior = T
+
+[test]
+save_dir=output/desi
+fatal_errors=T
+
+[output]
+filename = output/desi.txt
+
+[sample_rdh]
+file = utility/rescale_distances_rdh/rescale_distances_rdh.py
+
+[consistency]
+file = utility/consistency/consistency_interface.py
+
+[camb]
+file = boltzmann/camb/camb_interface.py
+mode = background
+feedback = 0
+AccuracyBoost = 1.0
+zmin_background = 0.
+zmax_background = 3.
+nz_background = 201
+use_ppf_w = T
+want_chistar = F
+
+
+[desi]
+file = likelihood/bao/desi1-dr1-arxiv/desi1_dr1_arxiv.py
+desi_data_sets = BGS,LRG1

examples/hsc-y3-shear-priors.ini

@@ -0,0 +1,15 @@
+[wl_photoz_errors]
+bias_1  =   gaussian 0. 0.024
+bias_2  =   gaussian 0. 0.022
+
+[shear_calibration_parameters]
+m1      =   gaussian 0. 0.01
+m2      =   gaussian 0. 0.01
+m3      =   gaussian 0. 0.01
+m4      =   gaussian 0. 0.01
+
+[psf_parameters]
+psf_alpha2 = gaussian 0.0 1.0
+psf_beta2 = gaussian 0.0 1.0
+psf_alpha4 = gaussian 0.0 1.0
+psf_beta4 = gaussian 0.0 1.0

examples/hsc-y3-shear-real-priors.ini

@@ -0,0 +1,15 @@
+[nz_sample_errors]
+bias_1  =   gaussian 0. 0.024
+bias_2  =   gaussian 0. 0.022
+
+[shear_calibration_parameters]
+m1      =   gaussian 0. 0.01
+m2      =   gaussian 0. 0.01
+m3      =   gaussian 0. 0.01
+m4      =   gaussian 0. 0.01
+
+[psf_systematics_parameters]
+psf_cor1_z1 = gaussian 0. 1.
+psf_cor2_z1 = gaussian 0. 1.
+psf_cor3_z1 = gaussian 0. 1.
+psf_cor4_z1 = gaussian 0. 1.