Merge pull request #69 from simonsobs/v_act_dr6

Modifications for ACT DR6 analysis

.github/workflows/testing.yml

@@ -12,7 +12,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: ["3.8", "3.9", "3.10", "3.11"]
+        python-version: ["3.9", "3.10", "3.11"]
 
     steps:
     - uses: actions/checkout@v4

mflike/MFLike.yaml

@@ -191,6 +191,29 @@ params:
       max: 10.60
     proposal: 0.6
     latex: T_d
+
+  beta_s:               #beta radio
+    value: -2.5
+    latex: \beta_s    
+  alpha_s:              #alpha radio
+    value: 1.0
+    latex: \alpha_s
+  T_effd:               #effective galactic dust temperature
+    value: 19.6
+    latex: T_{\mathrm{dust},\mathrm{eff}}
+  beta_d:               #beta galactic dust
+    value: 1.5
+    latex: \beta_\mathrm{dust}
+  alpha_dT:             #galactic dust ell slope for T
+    value: -0.6
+    latex: \alpha_{\mathrm{dust},T}
+  alpha_dE:             #galactic dust ell slope for E
+    value: -0.4
+    latex: \alpha_{\mathrm{dust},E}
+  alpha_p:              #CIB poisson ell slope
+    value: 1.0
+    latex: \alpha_p
+
   # Systematics
   bandint_shift_LAT_93:
     value: 0

mflike/mflike.py

@@ -80,23 +80,49 @@ def initialize(self):
             "a_pste",
             "xi",
             "T_d",
+            "beta_s",
+            "alpha_s",
+            "T_effd",
+            "beta_d",
+            "alpha_dT",
+            "alpha_dE",
+            "alpha_p"
         ]
+
 
         self.expected_params_nuis = ["calG_all"]
         for f in self.experiments:
             self.expected_params_nuis += [
                 f"bandint_shift_{f}",
-                f"cal_{f}",
                 f"calT_{f}",
+                f"cal_{f}",
                 f"calE_{f}",
                 f"alpha_{f}",
-            ]
 
+            ]
+
         self.ThFo = TheoryForge(self)
         self.log.info("Initialized!")
 
+    def initialize_non_sampled_params(self):
+        #allowing for systematic params not used in some analysis not to be sampled
+        self.non_sampled_params = {}
+        for exp in self.experiments:
+            self.non_sampled_params.update({f"calT_{exp}": 1.0, f"alpha_{exp}": 0.0})
+
     def initialize_with_params(self):
         # Check that the parameters are the right ones
+        self.initialize_non_sampled_params()
+
+        # Remove the parameters if it appears in the input/samples ones
+        for par in self.input_params:
+            self.non_sampled_params.pop(par, None)
+
+        # Finally set the list of nuisance params
+        self.expected_params_nuis = [
+            par for par in self.expected_params_nuis if par not in self.non_sampled_params
+        ]
+
         differences = are_different_params_lists(
             self.input_params,
             self.expected_params_fg + self.expected_params_nuis,
@@ -114,9 +140,18 @@ def logp(self, **params_values):
         params_values_nocosmo = {
             k: params_values[k] for k in self.expected_params_fg + self.expected_params_nuis
         }
+
         return self.loglike(cl, **params_values_nocosmo)
 
     def loglike(self, cl, **params_values_nocosmo):
+        # This is needed if someone calls the function without initializing the likelihood
+        # (typically a call with a precomputed Cl and no cobaya initialization steps e.g
+        # test_mflike)
+        if not hasattr(self, "non_sampled_params"):
+            self.initialize_non_sampled_params()
+
+        params_values_nocosmo = self.non_sampled_params | params_values_nocosmo
+
         ps_vec = self._get_power_spectra(cl, **params_values_nocosmo)
         delta = self.data_vec - ps_vec
         logp = -0.5 * (delta @ self.inv_cov @ delta)

mflike/tests/test_mflike.py

@@ -16,7 +16,7 @@
     "tau": 0.0544,
 }
 
-nuisance_params = {
+nuis_params = {
     "a_tSZ": 3.30,
     "a_kSZ": 1.60,
     "a_p": 6.90,
@@ -31,6 +31,13 @@
     "a_psee": 0,
     "a_pste": 0,
     "xi": 0.10,
+    "beta_s": -2.5,    
+    "alpha_s": 1,      
+    "T_effd": 19.6,    
+    "beta_d": 1.5,     
+    "alpha_dT": -0.6,  
+    "alpha_dE": -0.4,  
+    "alpha_p": 1, 
     "bandint_shift_LAT_93": 0,
     "bandint_shift_LAT_145": 0,
     "bandint_shift_LAT_225": 0,
@@ -94,7 +101,8 @@ def test_mflike(self):
                     },
                 }
             )
-            loglike = my_mflike.loglike(cl_dict, **nuisance_params)
+
+            loglike = my_mflike.loglike(cl_dict,  **nuis_params)
             self.assertAlmostEqual(-2 * (loglike - my_mflike.logp_const), chi2, 2)
 
     def test_cobaya(self):
@@ -113,14 +121,14 @@ def test_cobaya(self):
 
         model = get_model(info)
         my_mflike = model.likelihood["mflike.MFLike"]
-        chi2 = -2 * (model.loglikes(nuisance_params)[0] - my_mflike.logp_const)
+        chi2 = -2 * (model.loglikes(nuis_params)[0] - my_mflike.logp_const)
         self.assertAlmostEqual(chi2[0], chi2s["tt-te-et-ee"], 2)
 
     def test_top_hat_bandpasses(self):
         from copy import deepcopy
 
         # Let's vary values of bandint_shift parameters
-        params = deepcopy(nuisance_params)
+        params = deepcopy(nuis_params)
         params.update(
             {
                 k: {"prior": dict(min=0.9 * v, max=1.1 * v)}
@@ -142,7 +150,7 @@ def get_model(nsteps, bandwidth):
                     }
                 },
                 "theory": {"camb": {"extra_args": {"lens_potential_accuracy": 1}}},
-                "params": {**cosmo_params, **params},
+                "params": {**cosmo_params,  **params},
                 "packages_path": packages_path,
             }
 

mflike/theoryforge.py

@@ -4,7 +4,6 @@
 import numpy as np
 from cobaya.log import LoggedError
 
-
 # Converts from cmb temperature to differential source intensity
 # (see eq. 8 of https://arxiv.org/abs/1303.5070).
 # The bandpass transmission needs to be divided by
@@ -81,6 +80,7 @@ def __init__(self, mflike=None):
                         self.log, "One band has width = 0, set a positive width and run again"
                     )
 
+
     # Takes care of the bandpass construction. It returns a list of nu-transmittance
     # for each frequency or an array with the effective freqs.
     # bandpasses saved in the sacc file have to be divided by nu^2
@@ -139,7 +139,7 @@ def _bandpass_construction(self, **params):
     def get_modified_theory(self, Dls, **params):
         fg_params = {k: params[k] for k in self.expected_params_fg}
         nuis_params = {k: params[k] for k in self.expected_params_nuis}
-
+        
         # compute bandpasses at each step only if bandint_shift params are not null
         # and bandint_freqs has been computed at least once
         if np.all(
@@ -245,11 +245,11 @@ def _get_foreground_model(self, ell=None, freqs_order=None, **fg_params):
                 "temp": fg_params["T_d"],
                 "beta": fg_params["beta_p"],
             },
-            {"ell": ell_clp, "ell_0": ell_0clp, "alpha": 1},
+            {"ell": ell_clp, "ell_0": ell_0clp, "alpha": fg_params["alpha_p"]},
         )
         model["tt", "radio"] = fg_params["a_s"] * self.radio(
-            {"nu": self.bandint_freqs, "nu_0": nu_0, "beta": -0.5 - 2.0},
-            {"ell": ell_clp, "ell_0": ell_0clp, "alpha": 1},
+            {"nu": self.bandint_freqs, "nu_0": nu_0, "beta": fg_params["beta_s"]},
+            {"ell": ell_clp, "ell_0": ell_0clp, "alpha": fg_params["alpha_s"]},
         )
         model["tt", "tSZ"] = fg_params["a_tSZ"] * self.tsz(
             {"nu": self.bandint_freqs, "nu_0": nu_0},
@@ -265,8 +265,9 @@ def _get_foreground_model(self, ell=None, freqs_order=None, **fg_params):
             {"ell": ell, "ell_0": ell_0},
         )
         model["tt", "dust"] = fg_params["a_gtt"] * self.dust(
-            {"nu": self.bandint_freqs, "nu_0": nu_0, "temp": 19.6, "beta": 1.5},
-            {"ell": ell, "ell_0": 500.0, "alpha": -0.6},
+            {"nu": self.bandint_freqs, "nu_0": nu_0, 
+             "temp": fg_params["T_effd"], "beta": fg_params["beta_d"]},
+            {"ell": ell, "ell_0": 500.0, "alpha": fg_params["alpha_dT"]},
         )
         model["tt", "tSZ_and_CIB"] = self.tSZ_and_CIB(
             {
@@ -294,21 +295,23 @@ def _get_foreground_model(self, ell=None, freqs_order=None, **fg_params):
         )
 
         model["ee", "radio"] = fg_params["a_psee"] * self.radio(
-            {"nu": self.bandint_freqs, "nu_0": nu_0, "beta": -0.5 - 2.0},
-            {"ell": ell_clp, "ell_0": ell_0clp, "alpha": 1},
+            {"nu": self.bandint_freqs, "nu_0": nu_0, "beta": fg_params["beta_s"]},
+            {"ell": ell_clp, "ell_0": ell_0clp, "alpha": fg_params["alpha_s"]},
         )
         model["ee", "dust"] = fg_params["a_gee"] * self.dust(
-            {"nu": self.bandint_freqs, "nu_0": nu_0, "temp": 19.6, "beta": 1.5},
-            {"ell": ell, "ell_0": 500.0, "alpha": -0.4},
+            {"nu": self.bandint_freqs, "nu_0": nu_0, 
+             "temp": fg_params["T_effd"], "beta": fg_params["beta_d"]},
+            {"ell": ell, "ell_0": 500.0, "alpha": fg_params["alpha_dE"]},
         )
 
         model["te", "radio"] = fg_params["a_pste"] * self.radio(
-            {"nu": self.bandint_freqs, "nu_0": nu_0, "beta": -0.5 - 2.0},
-            {"ell": ell_clp, "ell_0": ell_0clp, "alpha": 1},
+            {"nu": self.bandint_freqs, "nu_0": nu_0, "beta": fg_params["beta_s"]},
+            {"ell": ell_clp, "ell_0": ell_0clp, "alpha": fg_params["alpha_s"]},
         )
         model["te", "dust"] = fg_params["a_gte"] * self.dust(
-            {"nu": self.bandint_freqs, "nu_0": nu_0, "temp": 19.6, "beta": 1.5},
-            {"ell": ell, "ell_0": 500.0, "alpha": -0.4},
+            {"nu": self.bandint_freqs, "nu_0": nu_0, 
+             "temp": fg_params["T_effd"], "beta": fg_params["beta_d"]},
+            {"ell": ell, "ell_0": 500.0, "alpha": fg_params["alpha_dE"]},
         )
 
         fg_dict = {}
@@ -348,10 +351,12 @@ def _get_foreground_model(self, ell=None, freqs_order=None, **fg_params):
     def _get_calibrated_spectra(self, dls_dict, **nuis_params):
         from syslibrary import syslib_mflike as syl
 
+        #allowing for not having calT_{exp} in the yaml
+
         cal_pars = {}
         if "tt" in self.requested_cls or "te" in self.requested_cls:
             cal = nuis_params["calG_all"] * np.array(
-                [nuis_params[f"cal_{exp}"] * nuis_params[f"calT_{exp}"] for exp in self.experiments]
+                [nuis_params[f"cal_{exp}"] * nuis_params.get(f"calT_{exp}",1) for exp in self.experiments]
             )
             cal_pars["t"] = 1 / cal
 
@@ -373,7 +378,9 @@ def _get_calibrated_spectra(self, dls_dict, **nuis_params):
     def _get_rotated_spectra(self, dls_dict, **nuis_params):
         from syslibrary import syslib_mflike as syl
 
-        rot_pars = [nuis_params[f"alpha_{exp}"] for exp in self.experiments]
+        #allowing for not having polarization angles in the yaml
+
+        rot_pars = [nuis_params.get(f"alpha_{exp}", 0) for exp in self.experiments]
 
         rot = syl.Rotation_alm(ell=self.l_bpws, spectra=dls_dict)