Mss003 tuning

sotodlib/toast/ops/mlmapmaker.py

@@ -461,6 +461,21 @@ def _wrap_obs(self, ob, dets, passinfo):
 
         return axobs
 
+    @function_timer
+    def _add_obs(self, mapmaker, name, axobs, nmat, signal_estimate):
+        """ Add data to the mapmaker
+
+        Singled out for more granular timing
+        """
+        mapmaker.add_obs(
+            name,
+            axobs,
+            deslope=self.deslope,
+            noise_model=nmat,
+            signal_estimate=signal_estimate,
+        )
+        return
+
     @function_timer
     def _init_mapmaker(
             self, mapmaker, signal_map, mapmaker_prev, eval_prev, comm, gcomm, prefix,
@@ -690,13 +705,8 @@ def _exec(self, data, detectors=None, **kwargs):
                 else:
                     signal_estimate = None
 
-                mapmaker.add_obs(
-                    ob.name,
-                    axobs,
-                    deslope=self.deslope,
-                    noise_model=nmat,
-                    signal_estimate=signal_estimate,
-                )
+                self._add_obs(mapmaker, ob.name, axobs, nmat, signal_estimate)
+
                 del axobs
                 del signal_estimate
 

sotodlib/toast/ops/splits.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2024-2024 Simons Observatory.
+# Copyright (c) 2024-2025 Simons Observatory.
 # Full license can be found in the top level "LICENSE" file.
 
 import os
@@ -12,7 +12,7 @@
 from toast.timing import function_timer
 from toast.intervals import IntervalList
 from toast.observation import default_values as defaults
-from toast.ops import Operator
+from toast.ops import Operator, FlagIntervals
 from toast.pixels_io_healpix import write_healpix_fits, write_healpix_hdf5
 from toast.pixels_io_wcs import write_wcs_fits
 from toast import qarray as qa
@@ -435,6 +435,47 @@ class Splits(Operator):
 
     splits = List([], help="The list of named splits to apply")
 
+    splits_as_flags = Bool(
+        False,
+        help="If True, splits are applied through flagging rather than views. "
+        "This can be more efficient if there are only a few samples in each "
+        "interval.",
+    )
+
+    shared_flags = Unicode(
+        defaults.shared_flags,
+        allow_none=True,
+        help="Observation shared key for telescope flags to use when "
+        "splits_as_flags is True",
+    )
+
+    shared_flag_mask = Int(
+        128,
+        help="Bit mask value for flagging when splits_as_flags is True",
+    )
+
+    write_hits = Bool(True, help="If True, write the hits map")
+
+    write_cov = Bool(
+        False, help="If True, write the white noise covariance matrices."
+    )
+
+    write_invcov = Bool(
+        True,
+        help="If True, write the inverse white noise covariance matrices.",
+    )
+
+    write_rcond = Bool(
+        False, help="If True, write the reciprocal condition numbers."
+    )
+
+    write_map = Bool(True, help="If True, write the filtered/destriped map")
+
+    write_noiseweighted_map = Bool(
+        False,
+        help="If True, write the noise-weighted map (for fast co-add)",
+    )
+
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
         return
@@ -516,6 +557,7 @@ def _exec(self, data, detectors=None, **kwargs):
         else:
             map_binner = self.mapmaker.binning
         pointing_view_save = map_binner.pixel_pointing.view
+        shared_flag_mask_save = map_binner.shared_flag_mask
 
         # If the pixel distribution does yet exit, we create it once prior to
         # doing any splits.
@@ -547,8 +589,23 @@ def _exec(self, data, detectors=None, **kwargs):
             for ob in data.obs:
                 spl.create_split(ob)
 
-            # Set mapmaking tools to use the current split interval list
-            map_binner.pixel_pointing.view = spl.split_intervals
+            if self.splits_as_flags:
+                # Split is applied through flagging, not as a view
+                # Flag samples outside the intervals by prefixing '~'
+                # to the view name
+                FlagIntervals(
+                    shared_flags=self.shared_flags,
+                    shared_flag_bytes=1,
+                    view_mask=[
+                        (f"~{spl.split_intervals}", np.uint8(self.shared_flag_mask))
+                    ],
+                    reset=True,
+                ).apply(data)
+                map_binner.shared_flag_mask |= self.shared_flag_mask
+            else:
+                # Set mapmaking tools to use the current split interval list
+                map_binner.pixel_pointing.view = spl.split_intervals
+
             if not map_binner.full_pointing:
                 # We are not using full pointing and so we clear the
                 # residual pointing for this split
@@ -574,17 +631,22 @@ def _exec(self, data, detectors=None, **kwargs):
         for k, v in mapmaker_save_traits.items():
             setattr(self.mapmaker, k, v)
         map_binner.pixel_pointing.view = pointing_view_save
+        map_binner.pixel_pointing.shared_flag_mask = shared_flag_mask_save
 
     def write_splits(self, data, split_name=None):
         """Write out all split products."""
         if not hasattr(self, "_split_obj"):
             msg = "No splits have been created yet, cannot write"
             raise RuntimeError(msg)
 
-        is_pix_wcs = hasattr(self.mapmaker.map_binning.pixel_pointing, "wcs")
+        if hasattr(self.mapmaker, "map_binning"):
+            pixel_pointing = self.mapmaker.map_binning.pixel_pointing
+        else:
+            pixel_pointing = self.mapmaker.binning.pixel_pointing
+        is_pix_wcs = hasattr(pixel_pointing, "wcs")
         is_hpix_nest = None
         if not is_pix_wcs:
-            is_hpix_nest = self.mapmaker.map_binning.pixel_pointing.nest
+            is_hpix_nest = pixel_pointing.nest
 
         if split_name is None:
             to_write = dict(self._split_obj)
@@ -593,8 +655,26 @@ def write_splits(self, data, split_name=None):
 
         for spname, spl in to_write.items():
             mname = f"{self.name}_{split_name}"
-            for prod in ["hits", "map", "invcov", "noiseweighted_map"]:
-                mkey = f"{mname}_{prod}"
+            for prod, binner_key, write in [
+                    ("hits", None, self.write_hits),
+                    ("cov", None, self.write_cov),
+                    ("invcov", None, self.write_invcov),
+                    ("rcond", None, self.write_rcond),
+                    ("map", "binned", self.write_map),
+                    ("noiseweighted_map", "noiseweighted", self.write_noiseweighted_map),
+            ]:
+                if not write:
+                    continue
+                if binner_key is not None:
+                    # get the product name from BinMap
+                    mkey = getattr(self.mapmaker.binning, binner_key)
+                else:
+                    # hits and covariance are not made by BinMap.
+                    # Try synthesizing the product name
+                    mkey = f"{mname}_{prod}"
+                if mkey not in data:
+                    msg = f"'{mkey}' not found in data.  Available keys are {data.keys()}"
+                    raise RuntimeError(msg)
                 if is_pix_wcs:
                     fname = os.path.join(
                         self.output_dir, f"{self.name}_{split_name}_{prod}.fits"

sotodlib/toast/workflows/proc_mapmaker.py

@@ -182,6 +182,11 @@ def mapmaker_select_noise_and_binner(job, otherargs, runargs, data):
 
         if job_tmpls.baselines.enabled:
             job_tmpls.noise_model = noise_model
+    elif job_ops.filterbin.enabled:
+        job_ops.filterbin.binning = job_ops.binner_final
+        job_ops.filterbin.binning.full_pointing = otherargs.full_pointing
+        job_ops.filterbin.det_data = job_ops.sim_noise.det_data
+        job_ops.filterbin.output_dir = otherargs.out_dir
 
 
 @workflow_timer

sotodlib/toast/workflows/sim_observe.py

@@ -98,6 +98,13 @@ def setup_simulate_observing(parser, operators):
         help="Realization index",
         type=int,
     )
+    parser.add_argument(
+        "--pwv_limit",
+        required=False,
+        type=float,
+        help="If set, discard observations with simulated PWV "
+        "higher than the limit [mm]",
+    )
 
     operators.append(
         toast.ops.SimGround(
@@ -230,6 +237,34 @@ def simulate_observing(job, otherargs, runargs, comm):
     job_ops.mem_count.prefix = "After Scan Simulation"
     job_ops.mem_count.apply(data)
 
+    if otherargs.pwv_limit is not None:
+        iobs = 0
+        ngood = 0
+        nbad = 0
+        while iobs < len(data.obs):
+            pwv = data.obs[iobs].telescope.site.weather.pwv.to_value(u.mm)
+            if pwv <= otherargs.pwv_limit:
+                ngood += 1
+                iobs += 1
+            else:
+                nbad += 1
+                del data.obs[iobs]
+                if len(data.obs) == 0:
+                    msg = (
+                        f"PWV limit = {otherargs.pwv_limit} mm rejected all "
+                        f"{nbad} observations assigned to this process"
+                    )
+                    raise RuntimeError(msg)
+        if toast_comm.comm_group_rank is not None:
+            nbad = toast_comm.comm_group_rank.allreduce(nbad)
+            ngood = toast_comm.comm_group_rank.allreduce(ngood)
+        log.info_rank(
+            f"  Discarded {nbad} / {ngood + nbad} observations "
+            f"with PWV > {otherargs.pwv_limit} mm in",
+            comm=comm,
+            timer=timer,
+        )
+
     # Apply LAT co-rotation
     if job_ops.corotate_lat.enabled:
         log.info_rank("  Running simulated LAT corotation...", comm=comm)

sotodlib/toast/workflows/sim_sky.py

@@ -1,4 +1,4 @@
-# Copyright (c) 2023-2024 Simons Observatory.
+# Copyright (c) 2023-2025 Simons Observatory.
 # Full license can be found in the top level "LICENSE" file.
 """Simulated detector response to sky signal.
 """
@@ -28,6 +28,7 @@ def setup_simulate_sky_map_signal(operators):
 
     """
     operators.append(toast.ops.ScanHealpixMap(name="scan_map", enabled=False))
+    operators.append(toast.ops.ScanHealpixDetectorMap(name="scan_detector_map", enabled=False))
     operators.append(
         toast.ops.StokesWeights(
             name="scan_map_weights",
@@ -89,54 +90,65 @@ def simulate_sky_map_signal(job, otherargs, runargs, data):
     # Configured operators for this job
     job_ops = job.operators
 
-    if job_ops.scan_map.enabled and job_ops.scan_wcs_map.enabled:
+    scan_healpix = job_ops.scan_map.enabled or job_ops.scan_detector_map.enabled
+    scan_wcs = job_ops.scan_wcs_map.enabled
+    scan_map = scan_healpix or scan_wcs
+
+    # Check for conflicting options
+
+    if scan_healpix and scan_wcs:
         msg = "Cannot scan from both healpix and WCS maps"
         log.error(msg)
         raise RuntimeError(msg)
 
     if job_ops.det_pointing_radec_sim is not job_ops.det_pointing_radec:
-        if job_ops.scan_map.enabled and not job_ops.scan_map_pixels.enabled:
+        if scan_healpix and not job_ops.scan_map_pixels.enabled:
             msg = "Simulation pointing is different from data reduction pointing. " \
                 " You must enable scan_map_pixels"
             log.error(msg)
             raise RuntimeError(msg)
-        if job_ops.scan_wcs_map.enabled and not job_ops.scan_wcs_map_pixels.enabled:
+        if scan_wcs and not job_ops.scan_wcs_map_pixels.enabled:
             msg = "Simulation pointing is different from data reduction pointing. " \
                 " You must enable scan_wcs_map_pixels"
             log.error(msg)
             raise RuntimeError(msg)
 
-    if job_ops.scan_map.enabled:
-        if job_ops.scan_map_pixels.enabled:
-            # We are using a custom pointing matrix
-            job_ops.scan_map_pixels.detector_pointing = job_ops.det_pointing_radec_sim
-            job_ops.scan_map.pixel_dist = "scan_map_pixel_dist"
-            job_ops.scan_map.pixel_pointing = job_ops.scan_map_pixels
-        else:
-            # We are using the same pointing matrix as the mapmaking
-            job_ops.scan_map.pixel_dist = job_ops.binner_final.pixel_dist
-            job_ops.scan_map.pixel_pointing = job.pixels_final
-        if job_ops.scan_map_weights.enabled:
-            job_ops.scan_map_weights.detector_pointing = job_ops.det_pointing_radec_sim
-            job_ops.scan_map.stokes_weights = job_ops.scan_map_weights
-        else:
-            job_ops.scan_map.stokes_weights = job_ops.weights_radec
-        job_ops.scan_map.save_pointing = otherargs.full_pointing
-        job_ops.scan_map.apply(data)
-        if job_ops.scan_map_pixels.enabled:
-            # Clean up our custom pointing
-            toast.ops.Delete(detdata=[
-                job_ops.scan_map_pixels.pixels,
-                job_ops.det_pointing_radec_sim.quats,
-            ]).apply(data)
-        if job_ops.scan_map_weights.enabled:
-            # Clean up our custom pointing
-            toast.ops.Delete(detdata=[
-                job_ops.scan_map_weights.weights,
-            ]).apply(data)
-        data.info()
-
-    if job_ops.scan_wcs_map.enabled:
+    # Scan signal from HEALPix maps
+
+    for scan_op in job_ops.scan_map, job_ops.scan_detector_map:
+        if scan_op.enabled:
+            if job_ops.scan_map_pixels.enabled:
+                # We are using a custom pointing matrix
+                job_ops.scan_map_pixels.detector_pointing = job_ops.det_pointing_radec_sim
+                scan_op.pixel_dist = "scan_map_pixel_dist"
+                scan_op.pixel_pointing = job_ops.scan_map_pixels
+            else:
+                # We are using the same pointing matrix as the mapmaking
+                scan_op.pixel_dist = job_ops.binner_final.pixel_dist
+                scan_op.pixel_pointing = job.pixels_final
+            if job_ops.scan_map_weights.enabled:
+                job_ops.scan_map_weights.detector_pointing = job_ops.det_pointing_radec_sim
+                scan_op.stokes_weights = job_ops.scan_map_weights
+            else:
+                scan_op.stokes_weights = job_ops.weights_radec
+            scan_op.save_pointing = otherargs.full_pointing
+            scan_op.apply(data)
+            if job_ops.scan_map_pixels.enabled:
+                # Clean up our custom pointing
+                toast.ops.Delete(detdata=[
+                    job_ops.scan_map_pixels.pixels,
+                    job_ops.det_pointing_radec_sim.quats,
+                ]).apply(data)
+            if job_ops.scan_map_weights.enabled:
+                # Clean up our custom pointing
+                toast.ops.Delete(detdata=[
+                    job_ops.scan_map_weights.weights,
+                ]).apply(data)
+            data.info()
+
+    # Scan signal from WCS maps
+
+    if scan_wcs:
         if job_ops.scan_wcs_map_pixels.enabled:
             # We are using a custom pointing matrix
             job_ops.scan_wcs_map_pixels.detector_pointing = job_ops.det_pointing_radec_sim