T and p beams

.github/workflows/wheels.yml

@@ -17,7 +17,7 @@ jobs:
 
       - name: Build sdist and wheel
         run: |
-          python -m pip install wheel
+          python -m pip install wheel setuptools
           python setup.py sdist bdist_wheel
 
       - uses: actions/upload-artifact@v3

pspipe_utils/best_fits.py

@@ -4,6 +4,7 @@
 import numpy as np
 from mflike import theoryforge as th_mflike
 from pspy import pspy_utils, so_spectra
+from pspipe_utils import misc
 
 def cmb_dict_from_file(f_name_cmb, lmax, spectra, lmin=2):
     """
@@ -31,18 +32,16 @@ def cmb_dict_from_file(f_name_cmb, lmax, spectra, lmin=2):
     return l_cmb, cmb_dict
 
 
-def fg_dict_from_files(f_name_fg, array_list, lmax, spectra, lmin=2, f_name_cmb=None):
+def fg_dict_from_files(f_name_fg, map_set_list, lmax, spectra, lmin=2, f_name_cmb=None):
     """
     create a fg power spectrum dict from files
 
     Parameters
     __________
     f_name_fg: string
       a template for the name of the fg power spectra files
-    array_list: list
-      list of all arrays
-    spectra: list
-      the list of spectra ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
+    map_set_list: list
+      list of all map set we will consider, format is {survey}_{array}
     lmax: integer
       the maximum multipole to consider (not inclusive)
     spectra: list
@@ -57,20 +56,20 @@ def fg_dict_from_files(f_name_fg, array_list, lmax, spectra, lmin=2, f_name_cmb=
         l_cmb, cmb_dict = cmb_dict_from_file(f_name_cmb, lmax, spectra, lmin)
 
     fg_dict = {}
-    for i, ar1 in enumerate(array_list):
-        for j, ar2 in enumerate(array_list):
-            if i > j: array_tuple = (ar2, ar1)
-            else: array_tuple = (ar1, ar2)
+    for i, ms_1 in enumerate(map_set_list):
+        for j, ms_2 in enumerate(map_set_list):
+            if i > j: ms_tuple = (ms_2, ms_1)
+            else: ms_tuple = (ms_1, ms_2)
 
-            l_fg, fg = so_spectra.read_ps(f_name_fg.format(*array_tuple), spectra=spectra)
+            l_fg, fg = so_spectra.read_ps(f_name_fg.format(*ms_tuple), spectra=spectra)
             id_fg = np.where((l_fg >= lmin) & (l_fg < lmax))
-            fg_dict[ar1, ar2] = {}
+            fg_dict[ms_1, ms_2] = {}
             for spec in spectra:
                 if i > j:
                     spec = spec[::-1]
-                fg_dict[ar1, ar2][spec] = fg[spec][id_fg]
+                fg_dict[ms_1, ms_2][spec] = fg[spec][id_fg]
                 if f_name_cmb is not None:
-                    fg_dict[ar1, ar2][spec] += cmb_dict[spec]
+                    fg_dict[ms_1, ms_2][spec] += cmb_dict[spec]
 
     l_fg = l_fg[id_fg]
 
@@ -120,14 +119,16 @@ def noise_dict_from_files(f_name_noise, sv_list, arrays, lmax, spectra, n_splits
     return l_noise, nl_dict
 
 
-def beam_dict_from_files(f_name_beam, sv_list, arrays, lmax, lmin=2):
+def beam_dict_from_files(f_name_beam_T, f_name_beam_pol, sv_list, arrays, lmax, lmin=2):
     """
     create a beam dict from files
 
     Parameters
     __________
-    f_name_beam: string
-        a template for the name of the beam files
+    f_name_beam_T: string
+        a template for the name of the temperature beam files
+    f_name_beam_pol: string
+        a template for the name of the polarisation beam files
     sv_list: list
         list of the surveys to consider
     arrays: dict
@@ -141,16 +142,22 @@ def beam_dict_from_files(f_name_beam, sv_list, arrays, lmax, lmin=2):
     bl_dict = {}
     for sv in sv_list:
         for ar in arrays[sv]:
-            l_beam, bl = pspy_utils.read_beam_file(f_name_beam.format(sv, ar))
+
+            l_beam, bl = misc.read_beams(f_name_beam_T.format(sv, ar),
+                                         f_name_beam_pol.format(sv, ar))
+
             id_beam = np.where((l_beam >= lmin) & (l_beam < lmax))
-            bl_dict[sv, ar] = bl[id_beam]
+
+            bl_dict[sv, ar] = {}
+            for field in ["T", "E", "B"]:
+                bl_dict[sv, ar][field] = bl[field][id_beam]
 
     l_beam = l_beam[id_beam]
 
     return l_beam, bl_dict
 
 
-def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, nl_dict=None, bl_dict=None):
+def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, delimiter="&", nl_dict=None, bl_dict=None):
     """
     This function prepare all best fit corresponding to the spec_name_list.
     the ps_all_th and nl_all_th are in particular useful for the analytical covariance computation
@@ -166,6 +173,13 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, nl_dict=N
         the cmb ps (format is [spec]
     fg_dict: dict of dict
         the fg ps (format is [sv1_ar1,sv2_ar2][spec])
+    spectra: list
+      the list of spectra ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
+    delimiter: string
+        serve to split the map_set into survey and array
+        a bit annoying to have to keep this, in principle we could work only
+        with map_set, but noise for different survey should be set to zero since
+        this entry doesn't exist
     nl_dict: dict of dict
         the noise ps (format is [sv, ar1, ar2][spec])
     bl_dict: dict
@@ -175,33 +189,37 @@ def get_all_best_fit(spec_name_list, l_th, cmb_dict, fg_dict, spectra, nl_dict=N
     ps_all_th, nl_all_th = {}, {}
 
     for spec_name in spec_name_list:
-        na, nb = spec_name.split("x")
-        if len(na.split("&")) == 2:
-            sv_a, ar_a = na.split("&")
-            sv_b, ar_b = nb.split("&")
+        ms_a, ms_b = spec_name.split("x")
+        if len(ms_a.split(delimiter)) == 2:
+            sv_a, ar_a = ms_a.split(delimiter)
+            sv_b, ar_b = ms_b.split(delimiter)
             noise_key_a = ar_a
             noise_key_b = ar_b
-        elif len(na.split("&")) == 3:
-            sv_a, ar_a, split_a = na.split("&")
-            sv_b, ar_b, split_b = nb.split("&")
+        elif len(ms_a.split(delimiter)) == 3:
+            sv_a, ar_a, split_a = ms_a.split(delimiter)
+            sv_b, ar_b, split_b = ms_b.split(delimiter)
             noise_key_a = f"{ar_a}_{split_a}"
             noise_key_b = f"{ar_b}_{split_b}"
 
         for spec in spectra:
-            ps_all_th[na, nb, spec] = cmb_dict[spec] + fg_dict[f"{sv_a}_{ar_a}", f"{sv_b}_{ar_b}"][spec]
+
+            ps_all_th[ms_a, ms_b, spec] = cmb_dict[spec] + fg_dict[f"{sv_a}_{ar_a}", f"{sv_b}_{ar_b}"][spec]
+            ps_all_th[ms_b, ms_a, spec] = ps_all_th[ms_a, ms_b, spec].copy()
 
             if bl_dict is not None:
-                ps_all_th[na, nb, spec] *= bl_dict[sv_a, ar_a] * bl_dict[sv_b, ar_b]
-
-            ps_all_th[nb, na, spec] = ps_all_th[na, nb, spec]
+                X, Y = spec
+                ps_all_th[ms_a, ms_b, spec] *=  bl_dict[sv_a, ar_a][X] * bl_dict[sv_b, ar_b][Y]
+                if ms_a != ms_b:
+                    # the if avoid a repetition in the case ms_a == ms_b
+                    ps_all_th[ms_b, ms_a, spec] *= bl_dict[sv_b, ar_b][X] * bl_dict[sv_a, ar_a][Y]
 
             if nl_dict is not None:
                 if sv_a == sv_b:
-                    nl_all_th[na, nb, spec] = nl_dict[sv_a, noise_key_a, noise_key_b][spec]
+                    nl_all_th[ms_a, ms_b, spec] = nl_dict[sv_a, noise_key_a, noise_key_b][spec]
                 else:
-                    nl_all_th[na, nb, spec] = 0
+                    nl_all_th[ms_a, ms_b, spec] = 0.0
 
-                nl_all_th[nb, na, spec] = nl_all_th[na, nb, spec]
+                nl_all_th[ms_b, ms_a, spec] = nl_all_th[ms_a, ms_b, spec]
 
     if nl_dict is not None:
         return l_th, ps_all_th, nl_all_th

pspipe_utils/misc.py

@@ -2,11 +2,60 @@
 Some general utility functions that do not belong to other files.
 """
 
+from pspy import pspy_utils
+from pixell import curvedsky
 
 def str_replace(my_str, old, new):
-    my_new_str = my_str
+    """
+    just like replace but check that the replacement actually happened
+
+    Parameters
+    __________
+    my_str: string
+        the string in which the replacment will happen
+    old: string
+        old part of the string to be replaced
+    new: string
+        what will replace old
+    """
+
     my_new_str = my_str.replace(old, new)
     if my_new_str == my_str:
         error = f" the name '{my_str}' does not contain '{old}' so I can't replace '{old}' by '{new}'"
         raise NameError(error)
     return my_new_str
+
+def read_beams(f_name_beam_T, f_name_beam_pol, lmax=None):
+    """
+    read T and pol beams and return a beam dictionnary with entry T, E, B
+
+    Parameters
+    __________
+    f_name_beam_T: string
+        the filename of the temperature beam file
+    f_name_beam_pol: string
+        the filename of the polarisation beam file
+    lmax : integer
+        the maximum multipole to consider (note that usually beam file start at l=0)
+    """
+
+    bl = {}
+    l, bl["T"] = pspy_utils.read_beam_file(f_name_beam_T, lmax=lmax)
+    l, bl["E"] = pspy_utils.read_beam_file(f_name_beam_pol, lmax=lmax)
+    bl["B"] = bl["E"]
+    return l, bl
+
+def apply_beams(alms, bl):
+    """
+    apply T and pol beams to alms
+
+    Parameters
+    __________
+    alms: 2d array
+        array of alms, alms[0]=alm_T, alms[1]=alm_E, alms[2]=alm_B
+    bl: dict
+        dictionnary containing T and E,B beams
+    """
+    for i, f in enumerate(["T", "E", "B"]):
+        alms[i] = curvedsky.almxfl(alms[i], bl[f])
+    return alms

pspipe_utils/pspipe_list.py

@@ -86,14 +86,14 @@ def get_covariances_list(dict):
 
     return ncovs, na_list, nb_list, nc_list, nd_list
 
-def get_spec_name_list(dict, char="&", kind=None, freq_pair=None, remove_same_ar_and_sv=False, return_nu_tag=False):
+def get_spec_name_list(dict, delimiter="&", kind=None, freq_pair=None, remove_same_ar_and_sv=False, return_nu_tag=False):
     """This function creates a list with the name of all spectra we consider
 
     Parameters
     ----------
     dict : dict
         the global dictionnary file used in pspipe
-    char: str
+    delimiter: str
         a character that separate the suvey and array name
     kind : str
         if "noise" or "auto" won't return
@@ -136,7 +136,7 @@ def get_spec_name_list(dict, char="&", kind=None, freq_pair=None, remove_same_ar
                     if remove_same_ar_and_sv == True:
                         if (sv1 == sv2) & (ar1 == ar2): continue
 
-                    spec_name_list += [f"{sv1}{char}{ar1}x{sv2}{char}{ar2}"]
+                    spec_name_list += [f"{sv1}{delimiter}{ar1}x{sv2}{delimiter}{ar2}"]
                     nu_tag_list += [(nu_tag1, nu_tag2)]
 
     if return_nu_tag == False:

setup.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.1.3
+current_version = 0.1.4
 commit = True
 tag = True
 

tutorials/tuto_covariances.py

@@ -2,10 +2,11 @@
 This script generate covariance matrices
 """
 import numpy as np
-from pspipe_utils import best_fits, pspipe_list
+from pspipe_utils import pspipe_list, best_fits
 from pspy import pspy_utils, so_cov, so_mcm, so_dict
 import time, sys
 
+
 d = so_dict.so_dict()
 d.read_from_file(sys.argv[1])
 
@@ -37,12 +38,13 @@
 f_name_cmb = tuto_data_dir + "/cmb.dat"
 f_name_noise = tuto_data_dir + "/mean_{}x{}_{}_noise.dat"
 f_name_fg = tuto_data_dir + "/fg_{}x{}.dat"
-f_name_beam = tuto_data_dir + "/beam_{}_{}.dat"
+f_name_beam_T = tuto_data_dir + "/beam_{}_{}_T.dat"
+f_name_beam_pol = tuto_data_dir + "/beam_{}_{}_pol.dat"
 
 l_cmb, cmb_dict = best_fits.cmb_dict_from_file(f_name_cmb, lmax, spectra)
 l_fg, fg_dict = best_fits.fg_dict_from_files(f_name_fg, arrays_list, lmax, spectra)
 l_noise, nl_dict = best_fits.noise_dict_from_files(f_name_noise,  surveys, arrays, lmax, spectra, n_splits=n_splits)
-l_beam, bl_dict = best_fits.beam_dict_from_files(f_name_beam, surveys, arrays, lmax)
+l_beam, bl_dict = best_fits.beam_dict_from_files(f_name_beam_T, f_name_beam_pol, surveys, arrays, lmax)
 
 l_cmb, ps_all_th, nl_all_th = best_fits.get_all_best_fit(spec_name_list,
                                                          l_cmb,