From b9c4e516b51ca9859c3ccf8f4cdcaeeae3ff203b Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 23:25:14 +0000 Subject: [PATCH 1/4] Add piecewise (per-segment) MAD scaling to surgical cleaner The surgical cleaner detrends its per-cell statistics piecewise in the frequency direction (subint_scaler) but then normalized them to sigma with a single global MAD across the whole band. On wideband receivers (e.g. Parkes UWL, 704-4032 MHz) the off-pulse RMS varies several-fold with Tsys across the band, so a global MAD inflates the apparent significance of channels in high-Tsys sub-bands and systematically over-flags them (and under-flags real RFI in quiet sub-bands). Compute the MAD per segment, using the same np.linspace segmentation as the piecewise detrend, so each channel is judged against the noise floor of its own sub-band: - clean_utils.py: new _piecewise_mad_scale() helper and a piecewise path in subint_scaler, gated by two new kwargs (piecewise_scale, defaulting to False, and subint_mad_numpieces, defaulting to the finest subint_numpieces value). Sparse segments (<16 unmasked points) and zero-MAD segments fall back to the global median/MAD. channel_scaler is already local in frequency and is unchanged. - surgical.py: register piecewise_scale and subint_mad_numpieces params and pass them through to comprehensive_stats. - default.cfg: add piecewise_scale=False, subint_mad_numpieces=None so the default behaviour is unchanged; wideband configs can opt in. Backward compatible: with the flag off, or when the MAD segmentation reduces to one piece, the output is bitwise-identical to the previous global MAD. Tests: Tsys-ramp unit test (global over-flags high-Tsys, piecewise balances; equal-excess spikes judged locally), exact backward-compat regression, sparse/all-masked/zero-MAD fallback edge cases, and a skipped psrchive integration placeholder. 183 passed, 1 skipped. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01ChmPorpBcsSQQeJ6Zaoutj --- coast_guard/clean_utils.py | 98 ++++++++++++++++- coast_guard/cleaners/surgical.py | 21 ++++ coast_guard/configurations/default.cfg | 2 +- tests/test_clean_utils.py | 142 +++++++++++++++++++++++++ 4 files changed, 259 insertions(+), 4 deletions(-) diff --git a/coast_guard/clean_utils.py b/coast_guard/clean_utils.py index ccd6a82..237fc67 100755 --- a/coast_guard/clean_utils.py +++ b/coast_guard/clean_utils.py @@ -110,13 +110,92 @@ def channel_scaler(array2d, **kwargs): return scaled +def _piecewise_mad_scale(detrended, numpieces, min_points=16): + """Scale a 1D detrended residual to robust sigma using a per-segment + median/MAD instead of a single global one. + + The residual is split into 'numpieces' roughly-equal segments using + the same np.linspace scheme as detrend(), and each segment is + centred on its own median and divided by its own MAD (scaled by + 1.4826). This keeps the definition of "sigma" local in frequency, so + a channel is judged against the noise floor of its own sub-band + rather than the whole band. On wideband receivers, where the + off-pulse RMS varies strongly with Tsys, this avoids inflating the + apparent significance of channels in high-Tsys sub-bands. + + Segments too sparse or degenerate to give a stable MAD fall back to + the global median/MAD of the whole residual: + * fewer than 'min_points' unmasked samples, or + * a zero MAD (all-equal or heavily-masked segment). + + Inputs: + detrended: A 1D (optionally masked) array of detrended residuals. + numpieces: The number of frequency segments to split into. + min_points: The minimum number of unmasked samples for a segment + to use its own median/MAD. (Default: 16) + + Output: + scaled: A 1D masked array of the residual in units of robust + sigma. Masked entries of the input remain masked. + """ + detrended = np.ma.masked_array(detrended, + mask=np.ma.getmaskarray(detrended)) + # Global fall-back statistics (used for sparse/degenerate segments). + global_median = np.ma.median(detrended) + global_mad = np.ma.median(np.abs(detrended - global_median)) + + scaled = np.ma.masked_array(np.empty(detrended.shape, dtype=float), + mask=np.ma.getmaskarray(detrended).copy()) + # Same segmentation scheme as detrend()'s numpieces splitting. + isplit = np.linspace(0, detrended.size, numpieces + 1, endpoint=1) + edges = np.round(isplit).astype(int) + for start, stop in zip(edges[:-1], edges[1:]): + seg = detrended[start:stop] + if np.ma.count(seg) < min_points: + # Too few unmasked points for a stable local MAD. + median, mad = global_median, global_mad + else: + median = np.ma.median(seg) + mad = np.ma.median(np.abs(seg - median)) + if mad == 0: + # Degenerate segment; local MAD would divide by zero. + median, mad = global_median, global_mad + scaled[start:stop] = (detrended[start:stop] - median)/(mad * 1.4826) + return scaled + + def subint_scaler(array2d, **kwargs): """For each sub-int detrend and scale it. + + By default each sub-int's detrended spectrum is scaled to robust + sigma using a single median/MAD taken over the whole band. On + wideband receivers the off-pulse RMS varies strongly with frequency + (Tsys), so a single global MAD inflates the apparent significance of + high-Tsys channels and systematically over-flags them. Passing + piecewise_scale=True instead computes the median/MAD per frequency + segment, using the same segmentation as the piecewise detrend, so the + scaling granularity matches the detrend and each channel is judged + against its own sub-band. + + Keyword arguments: + subint_order, subint_breakpoints, subint_numpieces: as for the + frequency-direction detrend. + piecewise_scale: If True, scale to sigma using a per-segment + (piecewise) median/MAD rather than a single global MAD. + (Default: False, i.e. original global behaviour.) + subint_mad_numpieces: Number of frequency segments to use for the + piecewise MAD. If None, the finest (last) subint_numpieces + value is used so the scaling lines up with the finest detrend + pass. Only used when piecewise_scale is True. (Default: None) """ # Grab key-word arguments. If not present use default configs. orders = kwargs.pop('subint_order', config.cfg.subint_order) breakpoints = kwargs.pop('subint_breakpoints', config.cfg.subint_breakpoints) numpieces = kwargs.pop('subint_numpieces', config.cfg.subint_numpieces) + # Piecewise-MAD scaling controls. Default to the backward-compatible + # global-MAD behaviour when these are absent. + piecewise_scale = kwargs.pop('piecewise_scale', False) + mad_numpieces = kwargs.pop('subint_mad_numpieces', None) if breakpoints is None: breakpoints = [[]]*len(orders) if numpieces is None: @@ -131,6 +210,14 @@ def subint_scaler(array2d, **kwargs): min(len(orders), len(breakpoints), len(numpieces))), errors.CoastGuardWarning) + # Resolve the segmentation for the piecewise MAD. When not given, default + # to the finest (last) detrend granularity so the scaling segments line up + # with the finest detrend pass. A value <= 1 means "one global segment", + # which reduces exactly to the original global-MAD behaviour. + if mad_numpieces is None: + last_numpieces = numpieces[-1] if len(numpieces) else None + mad_numpieces = last_numpieces if last_numpieces is not None else 1 + scaled = np.empty_like(array2d) nsubs = array2d.shape[0] for isub in np.arange(nsubs): @@ -138,9 +225,14 @@ def subint_scaler(array2d, **kwargs): for order, brkpnts, numpcs in zip(orders, breakpoints, numpieces): detrended = iterative_detrend(detrended, order=order, \ bp=brkpnts, numpieces=numpcs) - median = np.ma.median(detrended) - mad = np.ma.median(np.abs(detrended-median)) - scaled[isub,:] = (detrended-median)/(mad * 1.4826) # Scale MAD to be consistent with std for normal distribution + if piecewise_scale and mad_numpieces > 1: + # Local (per-segment) median/MAD, matching the piecewise detrend. + scaled[isub,:] = _piecewise_mad_scale(detrended, mad_numpieces) + else: + # Single global median/MAD over the whole band (original path). + median = np.ma.median(detrended) + mad = np.ma.median(np.abs(detrended-median)) + scaled[isub,:] = (detrended-median)/(mad * 1.4826) # Scale MAD to be consistent with std for normal distribution return scaled diff --git a/coast_guard/cleaners/surgical.py b/coast_guard/cleaners/surgical.py index a1301da..0cfc407 100755 --- a/coast_guard/cleaners/surgical.py +++ b/coast_guard/cleaners/surgical.py @@ -105,6 +105,25 @@ def _set_config_params(self): aliases=['iterations', 'iter'], nullable=True, help="Number of iterations to run the surgical cleaner [default = 1]") + self.configs.add_param('piecewise_scale', config_types.BoolVal, + aliases=['piecewise_scale'], + nullable=True, + help="If True, scale the frequency-direction " + "statistics to sigma using a per-segment " + "(piecewise) median/MAD that matches the " + "piecewise detrend, instead of a single " + "global MAD across the whole band. This " + "avoids over-flagging channels in high-Tsys " + "sub-bands on wideband receivers. " + "[default = False]") + self.configs.add_param('subint_mad_numpieces', config_types.IntVal, + aliases=['subint_mad_numpieces'], + nullable=True, + help="Number of frequency segments to use for " + "the piecewise MAD scaling (only used when " + "piecewise_scale is True). If None, the " + "finest subint_numpieces value is used. " + "[default = None]") self.parse_config_string(config.cfg.surgical_default_params) @@ -274,6 +293,8 @@ def _clean(self, ar): subint_order=self.configs.subint_order, \ subint_breakpoints=self.configs.subint_breakpoints, \ subint_numpieces=self.configs.subint_numpieces, \ + piecewise_scale=self.configs.piecewise_scale, \ + subint_mad_numpieces=self.configs.subint_mad_numpieces, \ aggressive=self.configs.aggressive \ ) if plot: diff --git a/coast_guard/configurations/default.cfg b/coast_guard/configurations/default.cfg index 5b3961f..f633807 100755 --- a/coast_guard/configurations/default.cfg +++ b/coast_guard/configurations/default.cfg @@ -10,7 +10,7 @@ combine_maxspan = 3600 # Max number of seconds a combined archive can span (psra combine_maxgap = 119 # Maximum gap between archives before starting a combined archive (psradd -G) # Cleaning -surgical_default_params = 'template=None,chan_breakpoints=None,chan_numpieces=1,chan_order=1,chanthresh=3,subint_breakpoints=1,subint_numpieces=1,subint_order=2;1,subintthresh=3,plot=None' +surgical_default_params = 'template=None,chan_breakpoints=None,chan_numpieces=1,chan_order=1,chanthresh=3,subint_breakpoints=1,subint_numpieces=1,subint_order=2;1,subintthresh=3,plot=None,piecewise_scale=False,subint_mad_numpieces=None' hotbins_default_params = 'calfrac=0.5,fscrunchfirst=False,iscal=False,onpulse=,threshold=5,tscrunchfirst=False' rcvrstd_default_params = 'badchans=None,badfreqs=None,badsubints=None,response=None,trimbw=0,trimfrac=0,trimnum=0' bandwagon_default_params = 'badchantol=0.8,badsubtol=0.8' diff --git a/tests/test_clean_utils.py b/tests/test_clean_utils.py index 258931b..4d2636d 100644 --- a/tests/test_clean_utils.py +++ b/tests/test_clean_utils.py @@ -277,3 +277,145 @@ def test_aggressive_uses_max_and_is_geq_average(self): aggressive=True) # max over diagnostics >= mean over diagnostics, elementwise. assert np.all(agg + 1e-9 >= avg) + + +# --------------------------------------------------------------------------- +# subint_scaler: piecewise-MAD scaling (frequency direction) +# --------------------------------------------------------------------------- +class TestSubintScalerBackwardCompat: + """With the flag off (or a single MAD segment) the piecewise path must + reduce EXACTLY to the original global-MAD behaviour.""" + + _KW = dict(subint_order=[1], subint_breakpoints=[[]], subint_numpieces=[1]) + + def test_flag_off_matches_global(self): + rng = np.random.default_rng(7) + arr = np.ma.masked_array(rng.normal(size=(4, 100))) + base = cu.subint_scaler(arr, piecewise_scale=False, **self._KW) + # Flag on but only one MAD segment -> identical to global. + same = cu.subint_scaler(arr, piecewise_scale=True, + subint_mad_numpieces=1, **self._KW) + npt.assert_array_equal(np.ma.getdata(base), np.ma.getdata(same)) + # Flag on, mad numpieces defaults to subint_numpieces[-1] == 1. + same2 = cu.subint_scaler(arr, piecewise_scale=True, **self._KW) + npt.assert_array_equal(np.ma.getdata(base), np.ma.getdata(same2)) + + def test_default_is_global(self): + # No piecewise kwargs at all -> original behaviour, reproduced by hand. + rng = np.random.default_rng(8) + arr = np.ma.masked_array(rng.normal(size=(2, 50))) + kw = dict(subint_order=[1], subint_breakpoints=[[]], subint_numpieces=[2]) + out = cu.subint_scaler(arr, **kw) + exp = np.empty_like(arr) + for i in range(arr.shape[0]): + d = cu.iterative_detrend(arr[i, :], order=1, bp=[], numpieces=2) + med = np.ma.median(d) + mad = np.ma.median(np.abs(d - med)) + exp[i, :] = (d - med) / (mad * 1.4826) + npt.assert_array_equal(np.ma.getdata(out), np.ma.getdata(exp)) + + +class TestSubintScalerPiecewise: + """Piecewise MAD keeps the 'sigma' definition local in frequency so a + Tsys gradient does not inflate high-Tsys channels' significance.""" + + _KW = dict(subint_order=[1], subint_breakpoints=[[]], subint_numpieces=[1]) + + @staticmethod + def _robust_spread(x): + x = np.asarray(x, dtype=float).ravel() + med = np.median(x) + return 1.4826 * np.median(np.abs(x - med)) + + def _tsys_ramp(self, nsub=4, nchan=512, seed=0): + # Per-channel off-pulse scatter ramps x5 across the band (1 -> 5). + rng = np.random.default_rng(seed) + sigma = 1.0 + 4.0 * np.arange(nchan) / (nchan - 1) + noise = rng.normal(size=(nsub, nchan)) * sigma[None, :] + return np.ma.masked_array(noise), sigma + + def test_global_overflags_high_tsys_piecewise_balances(self): + arr, _ = self._tsys_ramp() + g = np.ma.getdata(cu.subint_scaler(arr, piecewise_scale=False, **self._KW)) + p = np.ma.getdata(cu.subint_scaler(arr, piecewise_scale=True, + subint_mad_numpieces=8, **self._KW)) + low, high = slice(0, 64), slice(448, 512) + # Global MAD: high-Tsys channels have a much larger scaled spread -> + # they exceed any fixed threshold "just because" the noise is higher. + g_ratio = self._robust_spread(g[:, high]) / self._robust_spread(g[:, low]) + assert g_ratio > 3.0 + # Piecewise MAD: spread ~1 in both sub-bands (balanced, no over-flag). + p_low = self._robust_spread(p[:, low]) + p_high = self._robust_spread(p[:, high]) + assert 0.6 < (p_high / p_low) < 1.7 + npt.assert_allclose(p_low, 1.0, atol=0.35) + npt.assert_allclose(p_high, 1.0, atol=0.35) + + def test_equal_excess_outlier_judged_locally(self): + # Two spikes of EQUAL absolute excess, one in a low-Tsys sub-band + # (a real, locally-significant outlier) and one in a high-Tsys + # sub-band (within the local noise). + arr, _ = self._tsys_ramp(seed=1) + delta = 18.0 + lo_chan, hi_chan = 32, 480 + arr[0, lo_chan] += delta + arr[0, hi_chan] += delta + g = np.ma.getdata(cu.subint_scaler(arr, piecewise_scale=False, **self._KW)) + p = np.ma.getdata(cu.subint_scaler(arr, piecewise_scale=True, + subint_mad_numpieces=8, **self._KW)) + # Global MAD gives the two equal-excess spikes near-equal significance + # (it cannot tell the real RFI from the high-Tsys noise excursion). + assert abs(g[0, lo_chan]) > 5 and abs(g[0, hi_chan]) > 5 + assert 0.6 < abs(g[0, hi_chan]) / abs(g[0, lo_chan]) < 1.7 + # Piecewise MAD: the low-Tsys spike is a strong local outlier; the + # high-Tsys one sits within its sub-band's noise. + assert abs(p[0, lo_chan]) > 8 + assert abs(p[0, hi_chan]) < 5 + assert abs(p[0, lo_chan]) > 2 * abs(p[0, hi_chan]) + + +class TestPiecewiseMadEdgeCases: + """Sparse, all-masked and zero-MAD segments must fall back gracefully.""" + + def test_sparse_masked_and_degenerate_segments(self): + # 4 MAD segments of 16 channels each: + # seg0 fully masked, seg1 all-equal (MAD==0), seg2 <16 unmasked, + # seg3 normal noise. + rng = np.random.default_rng(4) + nchan = 64 + row = rng.normal(size=nchan) + row[16:32] = 5.0 # constant segment -> local MAD 0 + arr = np.ma.masked_array(row.reshape(1, nchan), + mask=np.zeros((1, nchan), dtype=bool)) + arr.mask[0, 0:16] = True # fully-masked segment + arr.mask[0, 32:34] = True # leaves 14 unmasked in seg2 (<16) + out = cu.subint_scaler(arr, subint_order=[0], subint_breakpoints=[[]], + subint_numpieces=[1], piecewise_scale=True, + subint_mad_numpieces=4) + # No NaN/Inf anywhere, even in the fallback segments. + assert np.all(np.isfinite(np.ma.filled(out, 0.0))) + # Masked input entries stay masked in the output. + assert np.all(out.mask[0, 0:16]) + assert np.all(out.mask[0, 32:34]) + + def test_helper_matches_global_for_one_piece(self): + # _piecewise_mad_scale with numpieces=1 == a single global MAD scale. + rng = np.random.default_rng(5) + d = np.ma.masked_array(rng.normal(size=200)) + med = np.ma.median(d) + mad = np.ma.median(np.abs(d - med)) + expected = (d - med) / (mad * 1.4826) + got = cu._piecewise_mad_scale(d, numpieces=1) + npt.assert_allclose(np.ma.getdata(got), np.ma.getdata(expected), atol=1e-12) + + +class TestPiecewiseMadIntegration: + """End-to-end check on a real wideband archive (needs psrchive + data).""" + + def test_uwl_tsys_gradient(self): + pytest.skip( + "Integration test: requires psrchive and a real UWL archive with a " + "strong Tsys gradient. Run the surgical cleaner with " + "piecewise_scale=True and assert the zapped-channel fraction in " + "high-Tsys sub-bands drops to that of low-Tsys sub-bands, while " + "injected narrowband RFI is still removed.") From 30c4e0947e3710be1cff3f44972df4a524925671 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 17 Jul 2026 00:09:06 +0000 Subject: [PATCH 2/4] Fix KeyError when a config omits the piecewise-scaling params add_param() registers a param's type but never stores its `default`, so a param only has a value if it appears in a parsed config string. The surgical cleaner reads self.configs.piecewise_scale and self.configs.subint_mad_numpieces unconditionally in _clean, but a receiver/custom config loaded via -C (set_override_config) replaces surgical_default_params with a string that does not list those keys -- bypassing the default.cfg copy -- so they were never set and the read raised KeyError. Every receiver config (UHF/UWL/L-band), i.e. the whole point of the wideband feature, hit this. Establish the safe (feature-off) defaults in _set_config_params by parsing 'piecewise_scale=False,subint_mad_numpieces=None' before parsing surgical_default_params. A config that lists the keys still overrides them; a config that omits them now falls back to feature-off instead of crashing. The redundant (and override-bypassed) copies are removed from default.cfg so the defaults have a single source of truth. Adds a regression test loading the surgical cleaner behind a receiver-style override both with and without the piecewise keys. 185 passed, 1 skipped. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01ChmPorpBcsSQQeJ6Zaoutj --- coast_guard/cleaners/surgical.py | 7 +++++ coast_guard/configurations/default.cfg | 2 +- tests/test_cleaners.py | 43 ++++++++++++++++++++++++++ 3 files changed, 51 insertions(+), 1 deletion(-) diff --git a/coast_guard/cleaners/surgical.py b/coast_guard/cleaners/surgical.py index 0cfc407..67e295f 100755 --- a/coast_guard/cleaners/surgical.py +++ b/coast_guard/cleaners/surgical.py @@ -124,6 +124,13 @@ def _set_config_params(self): "piecewise_scale is True). If None, the " "finest subint_numpieces value is used. " "[default = None]") + # Establish safe (feature-off) defaults for the optional piecewise + # params *before* parsing surgical_default_params. A receiver/custom + # config that overrides surgical_default_params need not list these + # keys; if it omits them they stay at these defaults (rather than + # raising KeyError when read in _clean), and if it sets them the + # override below wins. + self.parse_config_string('piecewise_scale=False,subint_mad_numpieces=None') self.parse_config_string(config.cfg.surgical_default_params) diff --git a/coast_guard/configurations/default.cfg b/coast_guard/configurations/default.cfg index ffa3e50..98200db 100755 --- a/coast_guard/configurations/default.cfg +++ b/coast_guard/configurations/default.cfg @@ -10,7 +10,7 @@ combine_maxspan = 3600 # Max number of seconds a combined archive can span (psra combine_maxgap = 119 # Maximum gap between archives before starting a combined archive (psradd -G) # Cleaning -surgical_default_params = 'template=None,chan_breakpoints=None,chan_numpieces=1,chan_order=1,chanthresh=3,subint_breakpoints=None,subint_numpieces=2;4,subint_order=2;1,subintthresh=3,plot=None,piecewise_scale=False,subint_mad_numpieces=None' +surgical_default_params = 'template=None,chan_breakpoints=None,chan_numpieces=1,chan_order=1,chanthresh=3,subint_breakpoints=None,subint_numpieces=2;4,subint_order=2;1,subintthresh=3,plot=None' hotbins_default_params = 'calfrac=0.5,fscrunchfirst=False,iscal=False,onpulse=,threshold=5,tscrunchfirst=False' rcvrstd_default_params = 'badchans=None,badfreqs=None,badsubints=None,response=None,trimbw=0,trimfrac=0,trimnum=0' bandwagon_default_params = 'badchantol=0.8,badsubtol=0.8' diff --git a/tests/test_cleaners.py b/tests/test_cleaners.py index a952ed6..e2d9079 100644 --- a/tests/test_cleaners.py +++ b/tests/test_cleaners.py @@ -115,3 +115,46 @@ def test_cleaner_get_config_string(self): cfgstr = cleaner.get_config_string() assert 'badchantol=' in cfgstr assert 'badsubtol=' in cfgstr + + +class TestSurgicalPiecewiseConfigResolution: + """Regression: the optional piecewise-scaling params must resolve even + when a receiver/custom config overrides ``surgical_default_params`` + without listing them. + + Reproduces the ``-C UHF/UWL/L`` path: ``set_override_config`` replaces + ``surgical_default_params`` with a receiver string that has no + ``piecewise_scale``/``subint_mad_numpieces`` keys. Previously the cleaner + never set them and reading them in ``_clean`` raised ``KeyError``; now they + default to feature-off but a config may still enable them. + """ + + # A receiver-style config string WITHOUT the new piecewise keys. + _RECEIVER = ("template=None,chan_breakpoints=None,chan_numpieces=1," + "chan_order=1,chanthresh=5,subint_breakpoints=None," + "subint_numpieces=8;16,subint_order=2;1,subintthresh=5," + "plot=None") + + def _load_with_override(self, override_string): + from coast_guard import config + cfg = config.cfg.get() + cfg.clear_overrides() + cfg.set_override_config('surgical_default_params', override_string) + try: + return cleaners.load_cleaner('surgical') + finally: + cfg.clear_overrides() + + def test_defaults_when_config_omits_them(self): + c = self._load_with_override(self._RECEIVER) + # No KeyError, and the safe (feature-off) defaults are in place. + assert c.configs.piecewise_scale is False + assert c.configs.subint_mad_numpieces is None + # The receiver's own params still took effect. + assert c.configs.subint_numpieces == [8, 16] + + def test_config_can_enable_piecewise(self): + c = self._load_with_override( + self._RECEIVER + ",piecewise_scale=True,subint_mad_numpieces=16") + assert c.configs.piecewise_scale is True + assert c.configs.subint_mad_numpieces == 16 From 9976ee537c0d0634234b982b6f7f93e510226e32 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 17 Jul 2026 01:24:34 +0000 Subject: [PATCH 3/4] Enable piecewise MAD in the UWL config and add a --piecewise CLI flag - UWL_3K_Murriyang.cfg: turn piecewise scaling on by default (piecewise_scale=True, subint_mad_numpieces=16 to match the finest subint_numpieces) since the 704-4032 MHz UWL band has a large Tsys gradient that a single global MAD over-flags. - clean_archive.py: add a -pw/--piecewise flag so pipelines that do not use a config file can enable it. apply_surgical_cleaner only appends piecewise_scale=True when the flag is given -- when absent it leaves the (default or -C) config's value untouched, so it never switches off a config that enables it (e.g. -C UWL). Adds a regression test that the shipped UWL config enables piecewise scaling. 186 passed, 1 skipped. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01ChmPorpBcsSQQeJ6Zaoutj --- clean_archive.py | 19 +++++++++++++++++-- .../receivers/UWL_3K_Murriyang.cfg | 6 +++++- tests/test_cleaners.py | 12 ++++++++++++ 3 files changed, 34 insertions(+), 3 deletions(-) diff --git a/clean_archive.py b/clean_archive.py index 004f1a7..24fbdb3 100644 --- a/clean_archive.py +++ b/clean_archive.py @@ -26,7 +26,7 @@ def apply_rcvrstd_cleaner(ar): rcvrstd_cleaner = cleaners.load_cleaner('rcvrstd') rcvrstd_cleaner.run(ar) -def apply_surgical_cleaner(ar, tmp, cthresh=None, sthresh=None, plot=False, aggressive=False, iterations=1): +def apply_surgical_cleaner(ar, tmp, cthresh=None, sthresh=None, plot=False, aggressive=False, iterations=1, piecewise_scale=False): """Apply the surgical cleaner to an archive in place. @@ -42,6 +42,10 @@ def apply_surgical_cleaner(ar, tmp, cthresh=None, sthresh=None, plot=False, aggr and algorithms. (Default: False) iterations: Number of times to run the surgical cleaner. (Default: 1) + piecewise_scale: If True, force per-segment (piecewise) MAD + scaling of the frequency-direction statistics on. If False, + the setting is left to the (default or -C) config, so this + does not switch off a config that enables it. (Default: False) Outputs: None - The archive is cleaned in place. @@ -62,6 +66,10 @@ def apply_surgical_cleaner(ar, tmp, cthresh=None, sthresh=None, plot=False, aggr param_parts.append("chanthresh={0}".format(cthresh)) if sthresh is not None: param_parts.append("subintthresh={0}".format(sthresh)) + if piecewise_scale: + # Only force it on; when the flag is absent leave the config's value + # untouched (so -C UWL, which enables it, is not overridden to off). + param_parts.append("piecewise_scale=True") surgical_parameters = ",".join(param_parts) surgical_cleaner.parse_config_string(surgical_parameters) surgical_cleaner.run(ar) @@ -121,6 +129,12 @@ def apply_bandwagon_cleaner(ar, badchantol=None, badsubtol=None): help="Whether to use more aggressive cleaning thresholds and algorithms") parser.add_argument("-i", "--iterations", type=int, dest="iterations", default=1, help="Number of iterations to run the surgical cleaner [default = 1]") + parser.add_argument("-pw", "--piecewise", dest='piecewise_scale', action='store_true', default=False, + help="Scale the frequency-direction statistics to sigma using a " + "per-segment (piecewise) MAD that matches the piecewise detrend, " + "instead of one global MAD across the whole band. Recommended for " + "wideband receivers with a strong Tsys gradient (already on in the " + "UWL config). [default = off]") parser.add_argument("-C", "--config", type=str, dest="config_path", default=None, help="Custom config file for misbehaving receivers. Inputting UHF or L will " "automatically load the MeerKAT config files for those receivers. Inputting " @@ -140,6 +154,7 @@ def apply_bandwagon_cleaner(ar, badchantol=None, badsubtol=None): output_path = args.output_path aggressive = args.aggressive iterations = args.iterations + piecewise_scale = args.piecewise_scale config_path = args.config_path #raise error if user tries to write multiple input files to one output name @@ -219,7 +234,7 @@ def apply_bandwagon_cleaner(ar, badchantol=None, badsubtol=None): subint_thresh if subint_thresh is not None else "cfg") apply_rcvrstd_cleaner(loaded_archive) - apply_surgical_cleaner(loaded_archive, template_path, cthresh=chan_thresh, sthresh=subint_thresh, plot=plot, aggressive=aggressive, iterations=iterations) + apply_surgical_cleaner(loaded_archive, template_path, cthresh=chan_thresh, sthresh=subint_thresh, plot=plot, aggressive=aggressive, iterations=iterations, piecewise_scale=piecewise_scale) apply_bandwagon_cleaner(loaded_archive, badchantol=badchantol, badsubtol=badsubtol) # Unload the Archive file diff --git a/coast_guard/configurations/receivers/UWL_3K_Murriyang.cfg b/coast_guard/configurations/receivers/UWL_3K_Murriyang.cfg index a3c6437..67eddd3 100644 --- a/coast_guard/configurations/receivers/UWL_3K_Murriyang.cfg +++ b/coast_guard/configurations/receivers/UWL_3K_Murriyang.cfg @@ -1,3 +1,7 @@ #crazy intense RFI mitigation for UWL fold-mode... -surgical_default_params = 'template=None,chan_breakpoints=None,chan_numpieces=1,chan_order=1,chanthresh=5,subint_breakpoints=None,subint_numpieces=8;16,subint_order=2;1,subintthresh=5,plot=None' +#piecewise_scale=True: the UWL band spans 704-4032 MHz with a large Tsys +#gradient, so scale each sub-int's stats to sigma per frequency segment +#(subint_mad_numpieces=16, matching the finest subint_numpieces) rather than +#with one global MAD, to avoid over-flagging high-Tsys sub-bands. +surgical_default_params = 'template=None,chan_breakpoints=None,chan_numpieces=1,chan_order=1,chanthresh=5,subint_breakpoints=None,subint_numpieces=8;16,subint_order=2;1,subintthresh=5,plot=None,piecewise_scale=True,subint_mad_numpieces=16' rcvrstd_default_params = 'badchans=None,badfreqs=962:1150;1857:1863;1725:1735;1080:1090;1022:1026;1918:1922;3070:3074,badsubints=None,response=704:4032,trimbw=0,trimfrac=0,trimnum=0' diff --git a/tests/test_cleaners.py b/tests/test_cleaners.py index e2d9079..bf21710 100644 --- a/tests/test_cleaners.py +++ b/tests/test_cleaners.py @@ -158,3 +158,15 @@ def test_config_can_enable_piecewise(self): self._RECEIVER + ",piecewise_scale=True,subint_mad_numpieces=16") assert c.configs.piecewise_scale is True assert c.configs.subint_mad_numpieces == 16 + + def test_shipped_uwl_config_enables_piecewise(self): + # The shipped UWL receiver config (loaded via -C UWL) should turn + # piecewise MAD scaling on by default. + import os + from coast_guard import config + uwl_path = os.path.join(config.base_config_dir, 'receivers', + 'UWL_3K_Murriyang.cfg') + overrides = config.read_file(uwl_path, required=True) + c = self._load_with_override(overrides['surgical_default_params']) + assert c.configs.piecewise_scale is True + assert c.configs.subint_mad_numpieces == 16 From ef7a726e368d96b996781dad9add5ce8c7825cd0 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 21 Jul 2026 22:53:03 +0000 Subject: [PATCH 4/4] Merge master (2D tools) and add template safety tests + nchan warning Brings the 2D-template subtraction and template-estimation safety changes (PR #7) into the piecewise-MAD branch. The merge was clean -- the 2D template loading (surgical._clean) and remove_profile1d changes sit in different regions from the piecewise subint_scaler work. On top of the merge: - clean_utils.check_template_nchan(): warn when a 2D (per-channel) template's channel count does not match the data. remove_profile_inplace indexes the template by data channel, so a mismatch mis-aligns the per-channel subtraction (or silently uses only the first data_nchan channels). Called from surgical._clean after the template shape is resolved. - Tests for the merged safety changes (all pure-Python, no psrchive): * remove_profile1d: all-zero template returns None and warns; the dot-product amplitude guess recovers the true scale; and the zero-median-template case (which broke the old median/median guess) now fits correctly. * check_template_nchan: warns on a 2D channel mismatch, stays quiet when the counts match and for a 1D template. 192 passed, 1 skipped. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01ChmPorpBcsSQQeJ6Zaoutj --- coast_guard/clean_utils.py | 29 ++++++++++++++ coast_guard/cleaners/surgical.py | 4 ++ tests/test_clean_utils.py | 67 ++++++++++++++++++++++++++++++++ 3 files changed, 100 insertions(+) diff --git a/coast_guard/clean_utils.py b/coast_guard/clean_utils.py index 2d0055b..3ac3d6f 100644 --- a/coast_guard/clean_utils.py +++ b/coast_guard/clean_utils.py @@ -507,6 +507,35 @@ def fft_rotate(data, bins): # original length (irfft otherwise assumes an even-length signal). return np.fft.irfft(phasor*np.fft.rfft(data), n=data.size) +def check_template_nchan(template, data_nchan): + """Warn if a 2D (per-channel) template's channel count does not match the + data being cleaned. + + remove_profile_inplace() indexes a 2D template by data channel + (template[ichan]), so a template with a different number of channels + than the data is mis-aligned -- or, if it has more channels, silently + uses only the first 'data_nchan' of them. A 1D template is broadcast to + every channel and is always fine, so no warning is issued for it. + + Inputs: + template: The template array (1D (nbin,) or 2D (nchan x nbin)). + data_nchan: The number of channels in the data being cleaned. + + Output: + None - Emits a CoastGuardWarning on a 2D channel-count mismatch. + """ + if np.ndim(template) <= 1: + return + template_nchan = np.shape(template)[0] + if template_nchan != data_nchan: + warnings.warn( + "2D template has %d channels but the data has %d. A 2D template is " + "indexed by data channel, so it must have the same number of " + "channels as the data; frequency-dependent subtraction may be " + "mis-aligned." % (template_nchan, data_nchan), + errors.CoastGuardWarning) + + def remove_profile1d(prof, isub, ichan, template, phs, return_params=False): """Fit and subtract a (rotated, scaled) template from a single profile. diff --git a/coast_guard/cleaners/surgical.py b/coast_guard/cleaners/surgical.py index 17fcc8f..2ea2de3 100644 --- a/coast_guard/cleaners/surgical.py +++ b/coast_guard/cleaners/surgical.py @@ -209,6 +209,10 @@ def _clean(self, ar): tuple(range(template_data.ndim - 1))).squeeze() template = template_data if template_data.ndim > 1 else template_phs + # A 2D template is indexed by data channel during subtraction, + # so warn if its channel count does not match the data. + clean_utils.check_template_nchan(template, patient.get_nchan()) + print('Estimating template and profile phase offset') if self.configs.template is None: diff --git a/tests/test_clean_utils.py b/tests/test_clean_utils.py index 4d2636d..6d452f0 100644 --- a/tests/test_clean_utils.py +++ b/tests/test_clean_utils.py @@ -9,7 +9,10 @@ import numpy.testing as npt import pytest +import warnings + from coast_guard import clean_utils as cu +from coast_guard import errors # --------------------------------------------------------------------------- @@ -419,3 +422,67 @@ def test_uwl_tsys_gradient(self): "piecewise_scale=True and assert the zapped-channel fraction in " "high-Tsys sub-bands drops to that of low-Tsys sub-bands, while " "injected narrowband RFI is still removed.") + + +# --------------------------------------------------------------------------- +# remove_profile1d: all-zero guard + dot-product amplitude estimate +# --------------------------------------------------------------------------- +class TestRemoveProfile1d: + @staticmethod + def _gaussian(nbin=64, centre=20, width=3.0): + x = np.arange(nbin, dtype=float) + return np.exp(-0.5 * ((x - centre) / width) ** 2) + + def test_all_zero_template_returns_none_and_warns(self): + prof = self._gaussian() + template = np.zeros(64) + with pytest.warns(errors.CoastGuardWarning): + (isub, ichan), resid = cu.remove_profile1d(prof, 2, 5, template, 0) + assert (isub, ichan) == (2, 5) + assert resid is None + # With return_params the params are None too. + with pytest.warns(errors.CoastGuardWarning): + out = cu.remove_profile1d(prof, 0, 0, template, 0, return_params=True) + assert out[1] is None and out[2] is None + + def test_dot_amplitude_recovers_scale(self): + template = self._gaussian() + prof = 3.0 * template + (_, _), resid, params = cu.remove_profile1d(prof, 0, 0, template, 0, + return_params=True) + npt.assert_allclose(params[0], 3.0, atol=1e-6) + npt.assert_allclose(resid, 0.0, atol=1e-6) + + def test_zero_median_template(self): + # Baseline-subtracted templates have (near) zero median, which made the + # old np.median(prof)/np.median(template) initial guess divide by ~0. + # The dot-product guess handles it. + template = self._gaussian() - np.median(self._gaussian()) + assert np.median(template) == 0.0 + prof = 2.5 * template + (_, _), resid, params = cu.remove_profile1d(prof, 0, 0, template, 0, + return_params=True) + npt.assert_allclose(params[0], 2.5, atol=1e-6) + npt.assert_allclose(resid, 0.0, atol=1e-6) + + +# --------------------------------------------------------------------------- +# check_template_nchan: warn on a 2D template / data channel-count mismatch +# --------------------------------------------------------------------------- +class TestCheckTemplateNchan: + def test_warns_on_2d_mismatch(self): + template2d = np.ones((8, 64)) + with pytest.warns(errors.CoastGuardWarning): + cu.check_template_nchan(template2d, 16) + + def test_no_warning_when_2d_matches(self): + template2d = np.ones((8, 64)) + with warnings.catch_warnings(): + warnings.simplefilter("error") # any warning -> test failure + cu.check_template_nchan(template2d, 8) + + def test_no_warning_for_1d_template(self): + template1d = np.ones(64) + with warnings.catch_warnings(): + warnings.simplefilter("error") + cu.check_template_nchan(template1d, 512)