diff --git a/draco/analysis/flagging.py b/draco/analysis/flagging.py index 673b71eb3..1baa1c05b 100644 --- a/draco/analysis/flagging.py +++ b/draco/analysis/flagging.py @@ -1037,13 +1037,12 @@ 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). @@ -1051,7 +1050,7 @@ class RFIStokesIMask(transform.ReduceVar): 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.