20240819 satp3 catchup (#984)

* updates for 20240819 satp3 catchup job

* remove changes added after catchup job finished

* add update to site_pipeline

* remove profiling calls

* adjust fill_glitches to address #740

* add error message for empty iir_params to partially address #969

* adjustment to new t2p leakage model handling

* remove repeated calculation get_gap_fill

* add lmfit to conf.py

* add lmfit to setup.py

* updates from satp1 site pipeline job

* remove lingering profile call

* fix processes.py

* adjustments to t2p, flags, gapfill

* fix docstring

* update glitch fill test

* update pca for compatibility

---------

Co-authored-by: Michael McCrackan <[email protected]>

docs/conf.py

@@ -103,7 +103,7 @@
                 'skyfield', 'h5py', 'pyfftw', 'scipy',
                 'toast', 'pixell', 'scikit', 'skimage', 'numdifftools',
                 'traitlets', 'ephem', 'influxdb', 'megham', 'detmap',
-                'sodetlib'):
+                'sodetlib', 'lmfit'):
     try:
         foo = import_module(missing)
     except ImportError:

setup.py

@@ -64,6 +64,7 @@
     'pyfftw',
     'numdifftools',
     'psycopg2-binary',
+    'lmfit',
 ]
 setup_opts["extras_require"] = {
     "site_pipeline": [

sotodlib/preprocess/processes.py

@@ -98,7 +98,7 @@ class Trends(_FracFlaggedMixIn, _Preprocess):
           signal: "signal" # optional
           calc:
             max_trend: 2.5
-            n_pieces: 10
+            t_piece: 100
           save: True
           plot: True
           select:
@@ -742,6 +742,8 @@ class GlitchFill(_Preprocess):
           nbuf: 10
           use_pca: False
           modes: 1
+          in_place: True
+          wrap: None
 
     .. autofunction:: sotodlib.tod_ops.gapfill.fill_glitches
     """
@@ -1164,6 +1166,7 @@ def plot(self, aman, proc_aman, filename):
                 band_aman = proc_aman[self.run_name].restrict('dets', aman.dets.vals[proc_aman[self.run_name][f'{band}_idx']], in_place=False)
                 plot_pcabounds(pca_aman, band_aman, filename=filename.replace('{name}', f'{ufm}_{band}_pca'), signal=self.signal, band=band, plot_ds_factor=self.plot_cfgs.get('plot_ds_factor', 20))
 
+
 class PTPFlags(_Preprocess):
     """Find detectors with anomalous peak-to-peak signal.
 
@@ -1256,6 +1259,7 @@ class EstimateT2P(_Preprocess):
             T_sig_name: 'dsT'
             Q_sig_name: 'demodQ'
             U_sig_name: 'demodU'
+            joint_fit: True
             trim_samps: 2000
             lpf_cfgs:
               type: 'sine2'

sotodlib/tod_ops/fft_ops.py

@@ -1,5 +1,6 @@
 """FFTs and related operations
 """
+import sys
 import numdifftools as ndt
 import numpy as np
 import pyfftw
@@ -151,6 +152,7 @@ def build_rfft_object(n_det, n, direction="FFTW_FORWARD", **kwargs):
 
     a = pyfftw.empty_aligned((n_det, n), dtype="float32")
     b = pyfftw.empty_aligned((n_det, (n + 2) // 2), dtype="complex64")
+
     if direction == "FFTW_FORWARD":
         t_fun = pyfftw.FFTW(a, b, direction=direction, **fftargs)
     elif direction == "FFTW_BACKWARD":
@@ -426,7 +428,8 @@ def fit_noise_model(
         pfit = np.polyfit(np.log10(f[f < lowf]), np.log10(p[f < lowf]), 1)
         fidx = np.argmin(np.abs(10 ** np.polyval(pfit, np.log10(f)) - wnest))
         p0 = [f[fidx], wnest, -pfit[0]]
-        res = minimize(neglnlike, p0, args=(f, p), method="Nelder-Mead")
+        bounds = [(0, None), (sys.float_info.min, None), (None, None)]
+        res = minimize(neglnlike, p0, args=(f, p), bounds=bounds, method="Nelder-Mead")
         try:
             Hfun = ndt.Hessian(lambda params: neglnlike(params, f, p), full_output=True)
             hessian_ndt, _ = Hfun(res["x"])

sotodlib/tod_ops/filters.py

@@ -521,6 +521,9 @@ def iir_filter(freqs, tod, b=None, a=None, fscale=1., iir_params=None,
                                              sub_iir_params['fscale'] != _fscale,])):
                             raise ValueError('iir parameters are not uniform.')
             iir_params = sub_iir_params
+            # check if iir_params from axis manager are None
+            if iir_params['a'] is None or iir_params['b'] is None:
+                raise ValueError('axis manager iir parameters are empty')
         try:
             a = iir_params['a']
             b = iir_params['b']

sotodlib/tod_ops/flags.py

@@ -16,7 +16,7 @@
 from . import fourier_filter 
 
 def get_det_bias_flags(aman, detcal=None, rfrac_range=(0.1, 0.7),
-                       psat_range=(0, 15), rn_range=None, si_nan=False,
+                       psat_range=None, rn_range=None, si_nan=False,
                        merge=True, overwrite=True,
                        name='det_bias_flags', full_output=False):
     """
@@ -34,6 +34,7 @@ def get_det_bias_flags(aman, detcal=None, rfrac_range=(0.1, 0.7),
     psat_range : Tuple
         Tuple (lower_bound, upper_bound) for P_SAT from IV analysis.
         P_SAT in the IV analysis is the bias power at 90% Rn in pW.
+        If None, no flags are not applied from P_SAT. 
     rn_range : Tuple
         Tuple (lower_bound, upper_bound) for r_n det selection.
     si_nan : bool
@@ -73,15 +74,17 @@ def get_det_bias_flags(aman, detcal=None, rfrac_range=(0.1, 0.7),
     ranges = [detcal.bg >= 0,
               detcal.r_tes > 0,
               detcal.r_frac >= rfrac_range[0],
-              detcal.r_frac <= rfrac_range[1],
-              detcal.p_sat*1e12 >= psat_range[0],
-              detcal.p_sat*1e12 <= psat_range[1]]
+              detcal.r_frac <= rfrac_range[1]
+             ]    
+    if psat_range is not None:
+        ranges.append(detcal.p_sat*1e12 >= psat_range[0])
+        ranges.append(detcal.p_sat*1e12 <= psat_range[1])
     if rn_range is not None:
         ranges.append(detcal.r_n >= rn_range[0])
         ranges.append(detcal.r_n <= rn_range[1])
     if si_nan:
         ranges.append(np.isnan(detcal.s_i) == False)
-    
+
     msk = ~(np.all(ranges, axis=0))
     # Expand mask to ndets x nsamps RangesMatrix
     if 'samps' in aman:
@@ -108,16 +111,22 @@ def get_det_bias_flags(aman, detcal=None, rfrac_range=(0.1, 0.7),
         ranges = [detcal.bg >= 0,
                   detcal.r_tes > 0,
                   detcal.r_frac >= rfrac_range[0],
-                  detcal.r_frac <= rfrac_range[1],
-                  detcal.p_sat*1e12 >= psat_range[0],
-                  detcal.p_sat*1e12 <= psat_range[1]]
+                  detcal.r_frac <= rfrac_range[1]
+                 ]
+        if psat_range is not None:
+            ranges.append(detcal.p_sat*1e12 >= psat_range[0])
+            ranges.append(detcal.p_sat*1e12 <= psat_range[1])
+
         for range in ranges:
             msk = ~(np.all([range], axis=0))
             msks.append(RangesMatrix([Ranges.ones_like(x) if Y
                                       else Ranges.zeros_like(x) for Y in msk]))
 
-        msk_names = ['bg', 'r_tes', 'r_frac_gt', 'r_frac_lt', 'p_sat_gt', 'p_sat_lt']
-
+        msk_names = ['bg', 'r_tes', 'r_frac_gt', 'r_frac_lt']
+
+        if psat_range is not None:
+            msk_names.extend(['p_sat_gt', 'p_sat_lt'])
+
         for i, msk in enumerate(msks):
             if 'samps' in aman:
                 msk_aman.wrap(f'{msk_names[i]}_flags', msk, [(0, 'dets'), (1, 'samps')])
@@ -401,7 +410,7 @@ def get_glitch_flags(aman,
 
 def get_trending_flags(aman,
                        max_trend=1.2,
-                       n_pieces=1,
+                       t_piece=500,
                        max_samples=500,
                        signal=None,
                        timestamps=None,
@@ -420,8 +429,8 @@ def get_trending_flags(aman,
         The tod
     max_trend : float
         Slope at which detectors are unlocked. The default is for use with phase units.
-    n_pieces : int
-        Number of pieces to cut the timestream in to to look for trends.
+    t_piece : float
+        Duration in seconds of each pieces to cut the timestream in to to look for trends
     max_samples : int
         Maximum samples to compute the slope with.
     signal : array
@@ -456,14 +465,21 @@ def get_trending_flags(aman,
     if timestamps is None:
         timestamps = aman.timestamps
     assert len(timestamps) == signal.shape[1]
-
+    
     # This helps with floating point precision
     # Not modifying inplace since we don't want to touch aman.timestamps
     timestamps = timestamps - timestamps[0]
 
     slopes = np.zeros((len(signal), 0))
     cut = np.zeros((len(signal), 0), dtype=bool)
     samp_edges = [0]
+
+    # Get sampling rate
+    fs = 1. / np.nanmedian(np.diff(timestamps))
+    n_samples_per_piece = int(t_piece * fs)
+    # How many pieces can timestamps be divided into
+    n_pieces = len(timestamps) // n_samples_per_piece
+
     for t, s in zip(
         np.array_split(timestamps, n_pieces), np.array_split(signal, n_pieces, 1)
     ):

sotodlib/tod_ops/gapfill.py

@@ -385,7 +385,7 @@ def get_contaminated_ranges(good_flags, bad_flags):
 
 
 def fill_glitches(aman, nbuf=10, use_pca=False, modes=3, signal=None,
-                  glitch_flags=None, wrap=True):
+                  glitch_flags=None, in_place=True, wrap=None):
     """
     This function fills pre-computed glitches provided by the caller in
     time-ordered data using either a polynomial (default) or PCA-based
@@ -407,9 +407,11 @@ def fill_glitches(aman, nbuf=10, use_pca=False, modes=3, signal=None,
     glitch_flags : RangesMatrix or None
         RangesMatrix containing flags to use for gap filling. If None then
         uses ``aman.flags.glitches``.
-    wrap : bool or str
-        If True wraps new field called ``gap_filled``, if False returns the
-        gap filled array, if a string wraps new field with provided name.
+    in_place : bool
+        If False it makes a copy of signal before gap filling and returns
+        the copy.
+    wrap : str or None
+        If not None, wrap the gap filled data into tod with this name.
 
     Returns
     -------
@@ -418,9 +420,15 @@ def fill_glitches(aman, nbuf=10, use_pca=False, modes=3, signal=None,
     """
     # Process Args
     if signal is None:
-        sig = np.copy(aman.signal)
+        if in_place:
+            sig = aman.signal
+        else:
+            sig = np.copy(aman.signal)
     else:
-        sig = np.copy(signal)
+        if in_place:
+            sig = signal
+        else:
+            sig = np.copy(signal)
 
     if glitch_flags is None:
         glitch_flags = aman.flags.glitches
@@ -437,26 +445,16 @@ def fill_glitches(aman, nbuf=10, use_pca=False, modes=3, signal=None,
                            f'{aman.dets.count}, setting modes = number of ' +
                            'detectors')
             modes = aman.dets.count
-        # fill with poly fill before PCA
-        gaps = get_gap_fill(aman, nbuf=nbuf, flags=glitch_flags,
-                            signal=np.float32(sig))
-        sig = gaps.swap(aman, signal=sig)
         # PCA fill
         mod = pca.get_pca_model(tod=aman, n_modes=modes,
                                 signal=sig)
         gfill = get_gap_model(tod=aman, model=mod, flags=glitch_flags)
         sig = gfill.swap(aman, signal=sig)
 
     # Wrap and Return
-    if isinstance(wrap, str):
+    if wrap is not None:
         if wrap in aman._assignments:
             aman.move(wrap, None)
         aman.wrap(wrap, sig, [(0, 'dets'), (1, 'samps')])
-        return sig
-    elif wrap:
-        if 'gap_filled' in aman._assignments:
-            aman.move('gap_filled', None)
-        aman.wrap('gap_filled', sig, [(0, 'dets'), (1, 'samps')])
-        return sig
-    else:
-        return sig
+
+    return sig

sotodlib/tod_ops/pca.py

@@ -364,5 +364,4 @@ def get_common_mode(
         raise ValueError("method flag must be median or average")
     if wrap is not None:
         tod.wrap(wrap, common_mode, [(0, 'samps')])
-    return common_mode
-
+    return common_mode