Fix handling of flagged detectors in noise estimation. Better support…

… for disabled templates in template matrix.

src/toast/ops/mapmaker_templates.py

@@ -211,6 +211,8 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
                     "You cannot currently initialize templates on device (please disable accel for this operator/pipeline)."
                 )
             for tmpl in self.templates:
+                if not tmpl.enabled:
+                    continue
                 if tmpl.view is None:
                     tmpl.view = self.view
                 tmpl.det_data_units = self.det_data_units
@@ -223,7 +225,8 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
         # Set the data we are using for this execution
         for tmpl in self.templates:
-            tmpl.det_data = self.det_data
+            if tmpl.enabled:
+                tmpl.det_data = self.det_data
 
         # We loop over detectors.  Internally, each template loops over observations
         # and ignores observations where the detector does not exist.
@@ -309,15 +312,16 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
             for d in all_dets:
                 for tmpl in self.templates:
-                    log.verbose(
-                        f"TemplateMatrix {d} add to signal {tmpl.name} (use_accel={use_accel})"
-                    )
-                    tmpl.add_to_signal(
-                        d,
-                        data[self.amplitudes][tmpl.name],
-                        use_accel=use_accel,
-                        **kwargs,
-                    )
+                    if tmpl.enabled:
+                        log.verbose(
+                            f"TemplateMatrix {d} add to signal {tmpl.name} (use_accel={use_accel})"
+                        )
+                        tmpl.add_to_signal(
+                            d,
+                            data[self.amplitudes][tmpl.name],
+                            use_accel=use_accel,
+                            **kwargs,
+                        )
         return
 
     def _finalize(self, data, use_accel=None, **kwargs):

src/toast/ops/noise_estimation.py

@@ -75,7 +75,7 @@ class NoiseEstim(Operator):
     )
 
     det_flag_mask = Int(
-        defaults.det_mask_nonscience,
+        defaults.det_mask_invalid,
         help="Bit mask value for optional detector flagging",
     )
 

src/toast/ops/noise_model.py

@@ -199,14 +199,21 @@ def _exec(self, data, detectors=None, **kwargs):
             nse_alpha = dict()
             nse_NET = dict()
             nse_indx = dict()
-            dets = ob.select_local_detectors(detectors, flagmask=self.det_flag_mask)
-            if len(dets) == 0:
-                # Nothing to do for this observation
-                continue
-            for det in dets:
+
+            # We are building a noise model with entries for all local detectors,
+            # even ones that are flagged.
+            for det in ob.local_detectors:
                 freqs = in_model.freq(det)
                 in_psd = in_model.psd(det)
                 cur_flag = ob.local_detector_flags[det]
+                nse_indx[det] = in_model.index(det)
+                nse_rate[det] = 2.0 * freqs[-1]
+                nse_NET[det] = 0.0 * np.sqrt(1.0 * in_psd.unit)
+                nse_fmin[det] = 0.0 * u.Hz
+                nse_fknee[det] = 0.0 * u.Hz
+                nse_alpha[det] = 0.0
+                if cur_flag & self.det_flag_mask != 0:
+                    continue
                 props = self._fit_log_psd(freqs, in_psd, guess=params)
                 if props["fit_result"].success:
                     # This was a good fit
@@ -220,8 +227,7 @@ def _exec(self, data, detectors=None, **kwargs):
                     log.verbose(msg)
                     new_flag = cur_flag | self.bad_fit_mask
                     ob.update_local_detector_flags({det: new_flag})
-                nse_indx[det] = in_model.index(det)
-                nse_rate[det] = 2.0 * freqs[-1]
+
                 nse_fmin[det] = props["fmin"]
                 nse_fknee[det] = props["fknee"]
                 nse_alpha[det] = props["alpha"]

src/toast/ops/stokes_weights/stokes_weights.py

@@ -248,7 +248,7 @@ def _exec(self, data, detectors=None, use_accel=None, **kwargs):
 
             if self.mode == "IQU":
                 det_gamma = np.zeros(len(dets), dtype=np.float64)
-                if self.hwp_angle is None:
+                if self.hwp_angle is None or self.hwp_angle not in ob.shared:
                     hwp_data = np.zeros(1, dtype=np.float64)
                 else:
                     hwp_data = ob.shared[self.hwp_angle].data

src/toast/templates/periodic.py

@@ -102,6 +102,7 @@ def _initialize(self, new_data):
         self._obs_max = list()
         self._obs_incr = list()
         self._obs_nbins = list()
+        total_bins = 0
         for ob in new_data.obs:
             if self.is_detdata_key:
                 if self.key not in ob.detdata:
@@ -147,6 +148,10 @@ def _initialize(self, new_data):
             else:
                 oincr = float(self.increment)
                 obins = int((omax - omin) / oincr)
+            if obins == 0 and ob.comm.group_rank == 0:
+                msg = f"Template {self.name}, obs {ob.name} has zero amplitude bins"
+                log.warning(msg)
+            total_bins += obins
             self._obs_nbins.append(obins)
             self._obs_incr.append(oincr)
             # Build up detector list
@@ -156,6 +161,11 @@ def _initialize(self, new_data):
 
         self._all_dets = list(all_dets.keys())
 
+        if total_bins == 0:
+            msg = f"Template {self.name} process group {new_data.comm.group}"
+            msg += f" has zero amplitude bins- change the binning size."
+            raise RuntimeError(msg)
+
         # During application of the template, we will be looping over detectors
         # in the outer loop.  So we pack the amplitudes by detector and then by
         # observation.  Compute the per-detector offsets into the amplitudes.

src/toast/tests/ops_noise_estim.py

@@ -311,15 +311,15 @@ def test_algorithm(self):
             for case, (input_freq, input_psd) in enumerate(zip(test_freq, test_psd)):
                 fit = fitter._fit_log_psd(input_freq, input_psd)
                 result = fit["fit_result"]
-                print(f"result solution = {result.x}")
-                print(f"result cost = {result.cost}")
-                print(f"result fun = {result.fun}")
-                print(f"result nfev = {result.nfev}")
-                print(f"result njev = {result.njev}")
-                print(f"result status = {result.status}")
-                print(f"result message = {result.message}")
-                print(f"result status = {result.status}")
-                print(f"result active_mask = {result.active_mask}")
+                # print(f"result solution = {result.x}")
+                # print(f"result cost = {result.cost}")
+                # print(f"result fun = {result.fun}")
+                # print(f"result nfev = {result.nfev}")
+                # print(f"result njev = {result.njev}")
+                # print(f"result status = {result.status}")
+                # print(f"result message = {result.message}")
+                # print(f"result status = {result.status}")
+                # print(f"result active_mask = {result.active_mask}")
                 fit_data = fitter._evaluate_model(
                     input_freq,
                     fit["fmin"],