fix(flagging): modify variance to add a significant boost during sola…

…r transit

draco/analysis/flagging.py

@@ -1037,21 +1037,20 @@ class RFIStokesIMask(transform.ReduceVar):
     frac_samples : float, optional
         Fraction of flagged samples in map space above which the entire
         time sample will be flagged. Default is 0.01.
-    include_sumthreshold : bool, optional
-        If True, apply a sumthreshold mask over the output `power` metric.
-        Default is True.
     max_m : int, optional
-        Maximum size of the SumThreshold window. Default is 16.
+        Maximum size of the SumThreshold window. Default is 64.
     nsigma : float, optional
         Initial threshold for SumThreshold. Default is 5.0.
+    solar_var_boost : float, optional
+        Variance boost during solar transit. Default is 1e4.
     bg_win_size : list, optional
         The size of the window used to estimate the background sky, provided
         as (number of frequency channels, number of time samples).
         Default is [11, 3].
     var_win_size : list, optional
         The size of the window used when estimating the variance, provided
         as (number of frequency channels, number of time samples).
-        Default is [3, 101].
+        Default is [3, 31].
     lowpass_cutoff : float, optional
         Angular cutoff of the ra lowpass filter. Default is 7.5, which
         corresponds to about 30 minutes of observation time.
@@ -1063,11 +1062,11 @@ class RFIStokesIMask(transform.ReduceVar):
     sigma_low = config.Property(proptype=float, default=2.0)
     frac_samples = config.Property(proptype=float, default=0.01)
 
-    include_sumthreshold = config.Property(proptype=bool, default=True)
-    max_m = config.Property(proptype=int, default=16)
+    max_m = config.Property(proptype=int, default=64)
     nsigma = config.Property(proptype=float, default=5.0)
+    solar_var_boost = config.Property(proptype=float, default=1e4)
     bg_win_size = config.list_type(int, length=2, default=[11, 3])
-    var_win_size = config.list_type(int, length=2, default=[3, 101])
+    var_win_size = config.list_type(int, length=2, default=[3, 31])
     lowpass_cutoff = config.Property(proptype=float, default=7.5)
 
     def setup(self, telescope):
@@ -1138,10 +1137,7 @@ def process(self, stream):
         power = mpiarray.MPIArray.wrap(power, axis=0).allgather()
 
         # Mask high power across frequencies
-        if self.include_sumthreshold:
-            summask = self.mask_multi_channel(power, mask, times)
-        else:
-            summask = np.zeros_like(mask)
+        summask = self.mask_multi_channel(power, mask, times)
 
         # Create the output containers
         if "ra" in stream.index_map:
@@ -1220,14 +1216,14 @@ def mask_single_channel(self, vis, weight, mask, freq, baselines, ra):
 
     def mask_multi_channel(self, power, mask, times):
         """Mask slow-moving narrow-band RFI."""
-        # Find times where there are bright sources transiting
+        # Find times where there are bright sources in the sky
+        # which should be treated differently
         source_flag = self._source_flag_hook(times)
-        # Avoid bright parts of the sky in the variance estimation
-        weight = ~mask & ~source_flag[np.newaxis]
+        sun_flag = self._solar_transit_hook(times)
 
-        # Calculate the weighted variance over time, excluding daytime
-        # and times where bright sources are transiting
-        wvar, ws = self.reduction(power, weight, axis=1)
+        # Calculate the weighted variance over time, excluding times
+        # flagged to have higher than normal variance
+        wvar, ws = self.reduction(power, ~mask & ~source_flag[np.newaxis], axis=1)
         # Get a smoothed estimate of the per-frequency variance
         wvar = tools.arPLS_1d(wvar, ws == 0, lam=1e1)[:, np.newaxis]
 
@@ -1241,7 +1237,10 @@ def mask_multi_channel(self, power, mask, times):
         # variance estimate
         med = weighted_median.weighted_median(p_med, (~mask).astype(p_med.dtype))
         rmed = medfilt(p_med, mask, size=self.var_win_size)
+        # Get the initial full variance using the lower variance estimate.
+        # Increase the variance estimate during solar transit
         var = wvar * rmed * tools.invert_no_zero(med)[:, np.newaxis]
+        var[:, sun_flag] *= self.solar_var_boost
 
         # Generate an RFI mask from the background-subtracted data
         summask = rfi.sumthreshold(
@@ -1326,6 +1325,21 @@ def _source_flag_hook(self, times):
         """
         return np.zeros_like(times, dtype=bool)
 
+    def _solar_transit_hook(self, times):
+        """Override to flag solar transit times.
+
+        Parameters
+        ----------
+        times : np.ndarray[float]
+            Array of timestamps.
+
+        Returns
+        -------
+        mask : np.ndarray[float]
+            Mask array. True will mask out a time sample.
+        """
+        return np.zeros_like(times, dtype=bool)
+
 
 class RFIMaskChisqHighDelay(task.SingleTask):
     """Mask frequencies and times with anomalous chi-squared test statistic.