Correct handling of 'pure_B' and 'tf_est_pure_B' cases

paramfiles/paramfile_SAT.yaml

@@ -77,8 +77,8 @@ general_pars:
   lmax: 600
   deltal: 10
   binning_file: "binning.npz"
-  pure_B: True
-  tf_est_pure_B: True
+  pure_B: False  # Compute data and covariance sims with purified B-modes
+  tf_est_pure_B: True  # Compute transfer function with purified B-modes
 
 ## Simulation related
 sim_pars:
@@ -99,9 +99,9 @@ sim_pars:
 
 ## Filtering transfer function related parameters
 tf_settings:
-  filtering_type: "m_filterer" # "m_filterer" or "toast"
-  m_cut: 30 # necessary only if `filtering_type` is set to "m_filterer"
-  toast: # necessary only if `filtering_type` is set to "toast"
+  filtering_type: "m_filterer"  # "m_filterer" or "toast"
+  m_cut: 30  # necessary only if `filtering_type` is set to "m_filterer"
+  toast:  # necessary only if `filtering_type` is set to "toast"
     schedule: "../outputs/schedules/schedule_sat.txt"
     instrument: "SAT1"
     band: "SAT_f090"

pipeline/mcmer.py

@@ -29,7 +29,9 @@ def mcmer(args):
     # Create dummy NaMaster fields, spin-2 is optionally purified.
     field_spin0 = nmt.NmtField(mask, None, spin=0)
     field_spin2 = nmt.NmtField(mask, None, spin=2)
-    if meta.pure_B:
+    # Whether we want to B-purify data or transfer function estimations,
+    # in either case we will need a purified mode coupling matrix.
+    if meta.pure_B or meta.tf_est_pure_B:
         field_spin2_pure = nmt.NmtField(mask, None, spin=2,
                                         purify_b=meta.pure_B)
 
@@ -45,15 +47,15 @@ def mcmer(args):
     print("Computing MCM")
     w = nmt.NmtWorkspace()
     w.compute_coupling_matrix(field_spin0, field_spin2, nmt_bins, is_teb=True)
-    if meta.pure_B:
+    if meta.pure_B or meta.tf_est_pure_B:
         w_pure = nmt.NmtWorkspace()
         w_pure.compute_coupling_matrix(field_spin0, field_spin2_pure, nmt_bins,
                                        is_teb=True)
 
     mcm = np.transpose(w.get_coupling_matrix().reshape([nl, nspec, nl, nspec]),
                        axes=[1, 0, 3, 2])
     mcm_binned = np.einsum('ij,kjlm->kilm', binner, mcm)
-    if meta.pure_B:
+    if meta.pure_B or meta.tf_est_pure_B:
         mcm_pure = np.transpose(
             w_pure.get_coupling_matrix().reshape([nl, nspec, nl, nspec]),
             axes=[1, 0, 3, 2]
@@ -87,8 +89,9 @@ def mcmer(args):
         spin2xspin2_binned=mcm_binned[3:, :, 3:, :]
     )
 
-    # If we B-purify, store also the un-beamed purified MCM
-    if meta.pure_B:
+    # If we B-purify for transfer function estimation, store the
+    # un-beamed purified mode coupling matrix
+    if meta.tf_est_pure_B:
         fname = f"{coupling_dir}/mcm_pure.npz"
         np.savez(
             fname,

pipeline/pcler.py

@@ -375,7 +375,8 @@ def pcler(args):
     if args.tf_val:
         inv_couplings = {}
         for filter_flag in ["filtered", "unfiltered"]:
-            couplings = np.load(f"{meta.coupling_directory}/couplings_{filter_flag}_pure.npz")  # noqa
+            pure_str = "_pure" if meta.tf_est_pure_B else ""
+            couplings = np.load(f"{meta.coupling_directory}/couplings_{filter_flag}{pure_str}.npz")  # noqa
             inv_couplings[filter_flag] = {
                 k1: couplings[f"inv_coupling_{k2}"].reshape([ncl*n_bins,
                                                              ncl*n_bins])

pipeline/transfer.py

@@ -61,11 +61,13 @@ def transfer(args):
         spin_pair: mcm[f"{spin_pair}_binned"] for spin_pair in spin_pairs
     }
 
-    # If we B-purify, load also the un-beamed purified mode coupling matrix
-    mcm_pure = np.load(f"{coupling_dir}/mcm_pure.npz")
-    mcms_dict_nobeam_pure = {
-        spin_pair: mcm_pure[f"{spin_pair}_binned"] for spin_pair in spin_pairs
-    }
+    # If we B-purify for transfer function estimation, load the un-beamed
+    # purified mode coupling matrix
+    if meta.tf_est_pure_B:
+        mcm_pure = np.load(f"{coupling_dir}/mcm_pure.npz")
+        mcms_dict_nobeam_pure = {
+            spin_pair: mcm_pure[f"{spin_pair}_binned"] for spin_pair in spin_pairs
+        }
 
     meta.timer.stop("Load mode-coupling matrices", verbose=True)
 
@@ -156,11 +158,15 @@ def transfer(args):
 
     for filter_tag in ["filtered", "unfiltered"]:
         couplings_nobeam = {}
-        couplings_nobeam_pure = {}
-        for couplings, mcms_dict in zip([couplings_nobeam,
-                                         couplings_nobeam_pure],
-                                        [mcms_dict_nobeam,
-                                         mcms_dict_nobeam_pure]):
+        couplings_list = [couplings_nobeam]
+        mcms_list = [mcms_dict_nobeam]
+
+        if meta.tf_est_pure_B:
+            couplings_nobeam_pure = {}
+            couplings_list.append(couplings_nobeam_pure)
+            mcms_list.append(mcms_dict_nobeam_pure)
+
+        for couplings, mcms_dict in zip(couplings_list, mcms_list):
             for spin_pair in spin_pairs:
                 if filter_tag == "filtered":
                     tbmcm = np.einsum('ijk,jklm->iklm', trans[spin_pair],
@@ -194,10 +200,11 @@ def transfer(args):
             f"{coupling_dir}/couplings_{filter_tag}.npz",
             **couplings_nobeam
         )
-        np.savez(
-            f"{coupling_dir}/couplings_{filter_tag}_pure.npz",
-            **couplings_nobeam_pure
-        )
+        if meta.tf_est_pure_B:
+            np.savez(
+                f"{coupling_dir}/couplings_{filter_tag}_pure.npz",
+                **couplings_nobeam_pure
+            )
 
     meta.timer.stop("Compute full coupling", verbose=True)