"""Configurable frequency- and spatial-mode flagging for data cubes.
Flaggers inspect a data cube and a noise proxy, usually a Stokes-I cube paired
with Stokes V, then either remove outlying channels / spatial modes
(``action='filter'``) or retain them with zero weight
(``action='zero_weight'``). :class:`FlaggerRunner` applies an ordered pipeline
and records the combined selection so it can be replayed on another cube::
runner = FlaggerRunner.load('flagger.ini')
filtered_i, filtered_v = runner.run(i_cube.copy(), v_cube.copy())
filtered_noise = runner.apply_last(noise_cube)
Flaggers may modify their input cubes; pass copies when the originals must be
preserved. Pipelines can be assembled in Python or loaded from an INI file.
"""
import configparser
import inspect
import matplotlib.pyplot as plt
import numpy as np
import tables
from matplotlib.colors import LogNorm
from . import datacube, fgfit, psutil, sphcube
[docs]
class BaseFlagger(psutil.SimpleConfig):
"""Base configuration shared by all flaggers.
Every flagger has a display ``name`` and an ``action``: ``'filter'``
removes selected modes, while ``'zero_weight'`` keeps their data and sets
their weights to zero.
"""
def __init__(self, **kargs):
"""Configure the flagger from keyword or string-valued settings."""
psutil.SimpleConfig.__init__(self)
self.add('name', 'Flagger', str)
self.add('action', 'filter', str)
self.parse_dict(kargs)
[docs]
def do_flag(self, i_cube, v_cube, verbose=True):
"""Apply this flagger to a data cube and its noise proxy.
Args:
i_cube: data cube used for Stokes-I statistics and geometry.
v_cube: matching noise-proxy cube, commonly Stokes V.
verbose (bool): report the number of selected modes.
Returns:
tuple: the flagged ``(i_cube, v_cube)``.
"""
raise NotImplementedError()
[docs]
def is_applicable(self, i_cube):
"""Whether this flagger supports the type of ``i_cube``."""
return True
[docs]
class BaseUVFlagger(BaseFlagger):
"""Base class for flaggers that select spatial modes."""
def __init__(self, **kargs):
BaseFlagger.__init__(self, **kargs)
[docs]
def do_flag(self, i_cube, v_cube, verbose=True):
"""Apply the spatial-mode mask returned by ``get_outliers``.
``get_outliers`` returns ``True`` for modes to flag. Filtering removes
those columns; zero-weighting retains the data and changes both cubes'
weight arrays in place.
Returns:
tuple: the flagged ``(i_cube, v_cube)``.
"""
uv_idx = self.get_outliers(i_cube, v_cube)
if uv_idx is not None:
n_filter = uv_idx.sum()
n_tot = len(uv_idx)
if verbose:
print(f'{self.name}: {self.action} {n_filter} / {n_tot} '
f'({n_filter / float(n_tot) * 100:.1f} %)')
if self.action == 'zero_weight':
i_cube.weights.data[:, uv_idx] = 0
v_cube.weights.data[:, uv_idx] = 0
elif self.action == 'filter':
i_cube.filter_uv_from_index(~uv_idx)
if v_cube != i_cube:
v_cube.filter_uv_from_index(~uv_idx)
return i_cube, v_cube
[docs]
class BaseFreqsFlagger(BaseFlagger):
"""Base class for flaggers that remove or zero-weight channels."""
def __init__(self, **kargs):
BaseFlagger.__init__(self, **kargs)
[docs]
def do_flag(self, i_cube, v_cube, verbose=True):
"""Apply the channel mask returned by ``get_outliers``.
``get_outliers`` returns ``True`` for channels to flag. Filtering
returns frequency-sliced cubes; zero-weighting retains the channels and
changes both cubes' weight arrays in place.
Returns:
tuple: the flagged ``(i_cube, v_cube)``.
"""
freqs_idx = self.get_outliers(i_cube, v_cube)
if freqs_idx is not None:
n_filter = freqs_idx.sum()
n_tot = len(freqs_idx)
if verbose:
print(f'{self.name}: {self.action} {n_filter} / {n_tot} '
f'({n_filter / float(n_tot) * 100:.1f} %)')
if self.action == 'zero_weight':
i_cube.weights.data[freqs_idx] = 0
v_cube.weights.data[freqs_idx] = 0
elif self.action == 'filter':
i_cube = i_cube.get_slice_from_idx(~freqs_idx)
if v_cube != i_cube:
v_cube = v_cube.get_slice_from_idx(~freqs_idx)
return i_cube, v_cube
[docs]
class FixedFreqsFlagger(BaseFreqsFlagger):
"""Flag configured frequencies or frequency ranges.
Frequencies are specified in MHz as comma-separated values or half-open
ranges, for example ``'121.3,123.1-124.0'``. By default every configured
channel is flagged. ``ratio_min_v_var`` can instead retain a configured
block unless its foreground-subtracted noise-proxy variance exceeds the
variance outside all configured blocks by that factor.
"""
def __init__(self, **kargs):
"""Configure fixed-frequency flagging.
Args:
**kargs: ``freqs`` is the MHz selection string;
``ratio_min_v_var`` enables variance-gated selection when
non-zero; ``ratio_min_v_var_poly_sub`` is the polynomial degree
used before comparing variances. ``name`` and ``action`` come
from :class:`BaseFlagger`.
"""
BaseFreqsFlagger.__init__(self, **kargs)
self.add('freqs', '', str)
self.add('ratio_min_v_var', 0, float)
self.add('ratio_min_v_var_poly_sub', 4, float)
self.parse_dict(kargs)
[docs]
def get_outliers(self, i_cube, v_cube):
"""Channels matching the configured ranges and variance criterion.
Returns:
ndarray: boolean mask of length ``n_freqs``; ``True`` means flag.
"""
fmhz = i_cube.freqs * 1e-6
fwidth = psutil.robust_freq_width(i_cube.freqs) * 1e-6
def fct_norm(f):
return np.round((f - fmhz[0]) / fwidth).astype(int)
norm_fmhz = fct_norm(fmhz)
idxs = []
for s in self.freqs.split(','):
if '-' in s:
a, b = s.split('-')
idxs.append(np.isin(norm_fmhz, np.arange(fct_norm(float(a)), fct_norm(float(b)))))
else:
idxs.append(np.isin(norm_fmhz, fct_norm(float(s))))
idx_full = np.sum(idxs, axis=0).astype(bool)
if self.ratio_min_v_var == 0:
return idx_full
idx_outside = ~idx_full
fitter = fgfit.PolyForegroundFit(self.ratio_min_v_var_poly_sub, 'power_poly')
fit_res = fitter.run(v_cube, v_cube)
v_var_outside = fit_res.sub.data[idx_outside].var()
outliers = np.zeros_like(fmhz)
for idx in idxs:
if fit_res.sub.data[idx].var() / v_var_outside > self.ratio_min_v_var:
outliers += idx
return outliers.astype(bool)
[docs]
class UVDirectionFlagger(BaseUVFlagger):
"""Flag Cartesian UV cells along a configured angular direction.
The direction is interpreted modulo 180 degrees. ``extend`` controls the
half-width of the selected UV line relative to the median U-cell spacing.
"""
def __init__(self, **kargs):
"""Configure the UV direction.
Args:
**kargs: ``direction_deg`` gives the direction in degrees and
``extend`` its width in UV-cell units; see :class:`BaseFlagger`
for ``name`` and ``action``.
"""
BaseUVFlagger.__init__(self, **kargs)
self.add('direction_deg', 0, float)
self.add('extend', 0.8, float)
self.parse_dict(kargs)
[docs]
def get_outliers(self, i_cube, v_cube):
"""UV samples intersecting the configured direction.
Returns:
ndarray: boolean mask of length ``n_vis``; ``True`` means flag.
"""
du = self.extend * np.median(np.diff(i_cube.uu))
d_rad = np.radians(self.direction_deg) % np.pi
theta = np.arctan2(i_cube.uu, i_cube.vv) % np.pi
d = (theta - d_rad - np.pi / 2) % np.pi
uv_idx = ((d > (np.pi / 2 - np.arctan(du / i_cube.ru)))
& (d < (np.pi / 2 + np.arctan(du / i_cube.ru))))
return uv_idx
[docs]
def is_applicable(self, i_cube):
"""Only Cartesian cubes provide the required UV coordinates."""
return isinstance(i_cube, datacube.CartDataCube)
[docs]
class SigmaClipper(psutil.SimpleConfig):
"""Shared robust detrending and one-sided sigma clipping.
The statistic is optionally polynomial-detrended after excluding an
initial set of extreme samples. Only positive residual outliers are
selected.
"""
def __init__(self, **kargs):
"""Configure robust clipping.
Args:
**kargs: ``nsigma`` sets the final MAD threshold;
``detrend_poly_deg`` sets the polynomial degree (0 disables
detrending); ``detrend_nsigma_clip`` sets the preliminary mask.
Concrete flaggers also use ``stokes`` to select I, V, dI, or dV
and ``sefd`` to choose SEFD rather than sample standard
deviation.
"""
psutil.SimpleConfig.__init__(self)
self.add('nsigma', 5, float)
self.add('detrend_poly_deg', 2, int)
self.add('detrend_nsigma_clip', 10, float)
self.add('stokes', 'V', str)
self.add('sefd', True, bool)
self.parse_dict(kargs)
[docs]
def get_sigma_clip_mask(self, x, y):
"""Detrend ``y`` and select high positive residuals.
Args:
x: coordinate used for polynomial detrending.
y: statistic to clip.
Returns:
ndarray: boolean outlier mask with the shape of ``y``.
"""
med = np.nanmedian(y)
rms = psutil.mad(y)
if self.detrend_poly_deg > 0:
mask = abs(y - med) < self.detrend_nsigma_clip * rms
mw_fct = np.poly1d(np.polyfit(x[mask], y[mask], self.detrend_poly_deg))
y = y - mw_fct(x)
med = np.nanmedian(y[mask])
rms = psutil.mad(y[mask])
return (y - med) > self.nsigma * rms
[docs]
class FreqsSigmaClipFlagger(SigmaClipper, BaseFreqsFlagger):
"""Sigma-clip channel noise or variance.
For each frequency, clip either the estimated SEFD (``sefd=True``) or the
standard deviation across spatial modes. ``stokes='I'`` uses ``i_cube``;
other values use the noise-proxy ``v_cube``.
"""
def __init__(self, **kargs):
"""Configure channel clipping; see :class:`SigmaClipper`."""
BaseFreqsFlagger.__init__(self, **kargs)
SigmaClipper.__init__(self, **kargs)
[docs]
def get_outliers(self, i_cube, v_cube):
"""Channels with anomalously high SEFD or variance.
Returns:
ndarray: boolean mask of length ``n_freqs``; ``True`` means flag.
"""
cube = i_cube if self.stokes == 'I' else v_cube
y = cube.estimate_freqs_sefd() if self.sefd else cube.data.std(axis=1)
return self.get_sigma_clip_mask(i_cube.freqs, y)
[docs]
def is_applicable(self, i_cube):
"""SEFD estimation requires a Cartesian cube; variance does not."""
if self.sefd:
return isinstance(i_cube, datacube.CartDataCube)
else:
return True
[docs]
class FreqsWeightsFlagger(BaseFreqsFlagger):
"""Flag channels whose median weight falls below a smooth trend.
A polynomial is fitted after excluding zero and five-MAD weight outliers.
Channels whose median effective weight divided by that trend is below
``ratio`` are selected.
"""
def __init__(self, **kargs):
"""Configure weight-trend flagging.
Args:
**kargs: ``ratio`` is the minimum trend-relative weight and
``trend_poly_deg`` the polynomial degree; see
:class:`BaseFlagger` for ``name`` and ``action``.
"""
BaseFreqsFlagger.__init__(self, **kargs)
self.add('ratio', 0.6, float)
self.add('trend_poly_deg', 2, int)
self.parse_dict(kargs)
[docs]
def get_outliers(self, i_cube, v_cube):
"""Channels whose median effective weight is below its fitted trend.
Returns:
ndarray: boolean mask of length ``n_freqs``; ``True`` means flag.
"""
weights = np.median(v_cube.weights.get(), axis=1)
mask = weights > 0
med = np.median(weights[mask])
rms = psutil.mad(weights[mask])
mask = abs(weights - med) < 5 * rms
mw_fct = np.poly1d(np.polyfit(i_cube.freqs[mask], weights[mask], self.trend_poly_deg))
return weights / mw_fct(i_cube.freqs) < self.ratio
[docs]
def is_applicable(self, i_cube):
"""Only Cartesian cubes are supported."""
return isinstance(i_cube, datacube.CartDataCube)
[docs]
class UVWeightsFlagger(BaseUVFlagger):
"""Flag UV samples with insufficient weight across the band.
For each UV sample, take its minimum effective weight over frequency and
compare it with the largest such minimum. Values below ``threshold`` times
that reference are selected.
"""
def __init__(self, **kargs):
"""Configure the relative minimum-weight ``threshold``."""
BaseUVFlagger.__init__(self, **kargs)
self.add('threshold', 0.01, float)
self.parse_dict(kargs)
[docs]
def get_outliers(self, i_cube, v_cube):
"""UV samples below the relative minimum-weight threshold.
Returns:
ndarray: boolean mask of length ``n_vis``; ``True`` means flag.
"""
min_weights = v_cube.weights.get().min(axis=0)
return min_weights < self.threshold * min_weights.max()
[docs]
def is_applicable(self, i_cube):
"""Only Cartesian cubes are supported."""
return isinstance(i_cube, datacube.CartDataCube)
[docs]
class UVSigmaClipFlagger(SigmaClipper, BaseUVFlagger):
"""Sigma-clip SEFD or variance across Cartesian UV samples.
The clipped statistic is evaluated per UV sample and detrended against
baseline radius. ``stokes`` selects ``'I'``, ``'V'``, ``'dI'``, or ``'dV'``;
the ``d`` variants first difference adjacent frequency channels.
"""
def __init__(self, **kargs):
"""Configure UV clipping; see :class:`SigmaClipper`."""
BaseUVFlagger.__init__(self, **kargs)
SigmaClipper.__init__(self, **kargs)
[docs]
def get_outliers(self, i_cube, v_cube):
"""UV samples with anomalously high SEFD or variance.
Returns:
ndarray: boolean mask of length ``n_vis``; ``True`` means flag.
"""
if self.stokes == 'I':
cube = i_cube
elif self.stokes == 'dI':
cube = i_cube.make_diff_cube()
elif self.stokes == 'dV':
cube = v_cube.make_diff_cube()
else:
cube = v_cube
y = cube.estimate_uv_sefd().data[0] if self.sefd else cube.data.std(axis=0)
return self.get_sigma_clip_mask(i_cube.ru, y)
[docs]
def is_applicable(self, i_cube):
"""Only Cartesian cubes provide the required UV statistics."""
return isinstance(i_cube, datacube.CartDataCube)
[docs]
class LMThetaMaxFlagger(BaseUVFlagger):
"""Flag spherical modes beyond a variance-derived angular boundary.
The boundary is inferred from frequency variance as a function of
``m / ell``. ``th_in`` defines the inner region used to set the reference
variance; ``relative_threshold`` selects the low-variance modes from which
the limiting angle is estimated.
"""
def __init__(self, **kargs):
"""Configure the inner-angle and relative-variance thresholds."""
BaseUVFlagger.__init__(self, **kargs)
self.add('th_in', 0.5, float)
self.add('relative_threshold', 0.001, float)
self.parse_dict(kargs)
[docs]
def get_outliers(self, i_cube, v_cube):
"""Spherical modes beyond the inferred angular boundary.
Returns:
ndarray: boolean mask of length ``n_modes``; ``True`` means flag.
"""
ml_r = i_cube.mm / (1. * i_cube.ll)
var_df_i = i_cube.data.var(axis=0)
th = self.relative_threshold * var_df_i[ml_r < self.th_in].mean()
th_max = np.arcsin(np.median(sorted(ml_r[var_df_i < th])[:20]))
idx_uv = i_cube.mm > np.clip(np.sin(th_max) * i_cube.ll, 1, max(i_cube.ll))
return idx_uv
[docs]
def is_applicable(self, i_cube):
"""Only spherical-harmonic cubes are supported."""
return isinstance(i_cube, sphcube.SphDataCube)
[docs]
class Flag:
"""Serializable frequency and spatial masks produced by a flagging run.
``idx_freqs`` and ``idx_uv`` select the channels and spatial modes kept
from the original cube. The corresponding ``*_zero_weights`` masks refer
to the retained grid and select samples whose weights are set to zero.
"""
def __init__(self, freqs, uu, vv, idx_uv, idx_uv_zero_weights, idx_freqs,
idx_freqs_zero_weights):
"""Create a flag selection for a specific cube geometry.
Args:
freqs: Original frequency coordinates in Hz.
uu: First original spatial coordinate. For Cartesian cubes this is
the U coordinate in wavelengths.
vv: Second original spatial coordinate, using the same convention
as ``uu``.
idx_uv: Boolean keep mask on the original spatial axis.
idx_uv_zero_weights: Boolean zero-weight mask on the retained
spatial axis.
idx_freqs: Boolean keep mask on the original frequency axis.
idx_freqs_zero_weights: Boolean zero-weight mask on the retained
frequency axis.
"""
self.freqs = freqs
self.uu = uu
self.vv = vv
self.idx_uv = idx_uv
self.idx_freqs = idx_freqs
self.idx_uv_zero_weights = idx_uv_zero_weights
self.idx_freqs_zero_weights = idx_freqs_zero_weights
[docs]
def save(self, filename):
"""Save the masks and their coordinate grids to HDF5.
Args:
filename: Destination HDF5 filename.
"""
with tables.open_file(filename, 'w') as h5_file:
group = h5_file.create_group("/", 'flag', 'Flag metadata')
h5_file.create_array(group, 'idx_uv', self.idx_uv, "Idx flag UV")
h5_file.create_array(group, 'idx_uv_zero_weights', self.idx_uv_zero_weights,
"Idx zero weights UV")
h5_file.create_array(group, 'idx_freqs', self.idx_freqs, "Idx flag freqs")
h5_file.create_array(group, 'idx_freqs_zero_weights', self.idx_freqs_zero_weights,
"Idx zero weights freqs")
h5_file.create_array(group, 'freqs', self.freqs, "Frequencies (Hz)")
h5_file.create_array(group, 'uu', self.uu, "U (lambda)")
h5_file.create_array(group, 'vv', self.vv, "V (lambda)")
[docs]
@staticmethod
def load(filename):
"""Load a :class:`Flag` from HDF5.
Args:
filename: HDF5 file created by :meth:`save`.
Returns:
The stored :class:`Flag`.
"""
with tables.open_file(filename, 'r') as h5_file:
idx_uv = h5_file.root.flag.idx_uv.read()
idx_uv_zero_weights = h5_file.root.flag.idx_uv_zero_weights.read()
idx_freqs = h5_file.root.flag.idx_freqs.read()
idx_freqs_zero_weights = h5_file.root.flag.idx_freqs_zero_weights.read()
freqs = h5_file.root.flag.freqs.read()
uu = h5_file.root.flag.uu.read()
vv = h5_file.root.flag.vv.read()
return Flag(freqs, uu, vv, idx_uv, idx_uv_zero_weights, idx_freqs, idx_freqs_zero_weights)
[docs]
def apply(self, cube):
"""Apply these masks to a compatible cube.
The input is copied. Rejected channels and spatial modes are removed,
then the zero-weight masks are applied to the retained grid.
Args:
cube: Cube with the same original frequency and spatial layout as
the cube used to create this flag.
Returns:
A flagged copy of ``cube``.
"""
assert np.allclose(cube.freqs, self.freqs)
assert len(cube.ru) == len(self.uu)
cube = cube.copy()
cube.filter_uv_from_index(self.idx_uv)
cube = cube.get_slice_from_idx(self.idx_freqs)
if cube.weights is not None:
cube.weights.data[self.idx_freqs_zero_weights] = 0
cube.weights.data[:, self.idx_uv_zero_weights] = 0
return cube
[docs]
class FlaggerRunner:
"""Run an ordered flagger pipeline on two matched cubes.
Every flagger derives a mask from the Stokes-I and noise cubes and applies
it to both, keeping their sampling identical. The runner retains the
inputs, outputs, and combined :class:`Flag` from the most recent run for
plotting or reuse.
The pipeline can modify the supplied cubes. Pass copies to :meth:`run`
when their original state must be retained.
"""
def __init__(self, verbose=True, min_weights=2):
"""Create an empty pipeline.
Args:
verbose: Report how many samples each operation affects.
min_weights: Remove samples below this weight before running the
configured flaggers. Set to ``None`` to disable this step.
"""
self.flaggers = []
self.original_i_cube = None
self.original_v_cube = None
self.filtered_i_cube = None
self.filtered_v_cube = None
self.flag = None
self.verbose = verbose
self.min_weights = min_weights
[docs]
def add(self, flagger):
"""Append a flagger to the pipeline.
Args:
flagger: A :class:`BaseFlagger` instance.
"""
self.flaggers.append(flagger)
[docs]
def run(self, i_cube, v_cube):
"""Run the pipeline on two cubes with matching sampling.
Flaggers that do not support the supplied cube type or metadata are
skipped. Filtering and zero-weight selections are applied identically
to both cubes.
Args:
i_cube: Stokes-I data cube.
v_cube: Noise or Stokes-V cube with the same geometry.
Returns:
Tuple containing the flagged Stokes-I and noise cubes.
Note:
The input cubes can be modified. Use ``cube.copy()`` when their
original state is needed later.
"""
i_cube_before_flagging = i_cube.copy()
if self.min_weights is not None:
i_cube.filter_min_weight(self.min_weights, verbose=self.verbose)
v_cube.filter_min_weight(self.min_weights, verbose=self.verbose)
self.original_i_cube = i_cube.copy()
self.original_v_cube = v_cube.copy()
for flagger in self.flaggers:
if flagger.is_applicable(i_cube):
i_cube, v_cube = flagger.do_flag(i_cube, v_cube, verbose=self.verbose)
self.filtered_i_cube = i_cube
self.filtered_v_cube = v_cube
idx_freqs, idx_uv, _, _ = datacube.get_common_idx(i_cube_before_flagging,
self.filtered_i_cube)
if self.filtered_i_cube.weights is not None:
idx_uv_zero_weights = (self.filtered_i_cube.weights.get().sum(axis=0) == 0)
idx_freqs_zero_weights = (self.filtered_i_cube.weights.get().sum(axis=1) == 0)
else:
idx_uv_zero_weights = None
idx_freqs_zero_weights = None
if isinstance(i_cube_before_flagging, datacube.CartDataCube):
self.flag = Flag(
i_cube_before_flagging.freqs, i_cube_before_flagging.uu,
i_cube_before_flagging.vv, idx_uv, idx_uv_zero_weights, idx_freqs,
idx_freqs_zero_weights)
elif isinstance(i_cube_before_flagging, sphcube.SphDataCube):
self.flag = Flag(
i_cube_before_flagging.freqs, i_cube_before_flagging.ll,
i_cube_before_flagging.mm, idx_uv, idx_uv_zero_weights, idx_freqs,
idx_freqs_zero_weights)
else:
raise ValueError('Cube is not of a supported format')
return self.filtered_i_cube, self.filtered_v_cube
[docs]
def apply_last(self, cube):
"""Apply the most recently generated flag to a compatible cube.
Args:
cube: Cube with the same original geometry as the last inputs.
Returns:
A flagged copy of ``cube``.
"""
return self.flag.apply(cube)
[docs]
def plot(self, figsize=(10, 12), **fig_kargs):
"""Plot diagnostics from the most recent run.
Removed samples are marked with crosses and retained samples whose
weights were set to zero with plus signs.
Args:
figsize: Matplotlib figure size.
**fig_kargs: Additional arguments passed to
:func:`matplotlib.pyplot.subplots`.
Returns:
The Matplotlib figure containing frequency and spatial
diagnostics.
"""
fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(ncols=2, nrows=3, figsize=figsize,
**fig_kargs)
idx1, idx1_uv, _, _ = datacube.get_common_idx(self.original_i_cube, self.filtered_i_cube)
freqs_idx_filtered = ~idx1
uv_idx_filtered = ~idx1_uv
if self.filtered_i_cube.weights is not None:
freqs_idx_zero_w = self.filtered_i_cube.weights.data.sum(axis=1) == 0
uv_idx_zero_w = self.filtered_i_cube.weights.data.sum(axis=0) == 0
ax1.plot(self.original_i_cube.freqs * 1e-6, self.original_i_cube.data.var(axis=1),
c=psutil.lblue, label='I before')
ax1.plot(self.original_i_cube.freqs * 1e-6, self.original_v_cube.data.var(axis=1),
c=psutil.lorange, label='V before')
ax1.plot(self.filtered_i_cube.freqs * 1e-6, self.filtered_i_cube.data.var(axis=1),
c=psutil.dblue, label='I after')
ax1.plot(self.filtered_v_cube.freqs * 1e-6, self.filtered_v_cube.data.var(axis=1),
c=psutil.dorange, label='V after')
ax1.plot(self.original_i_cube.freqs[freqs_idx_filtered] * 1e-6,
self.original_v_cube.data.var(axis=1)[freqs_idx_filtered],
marker='x', ls='', c=psutil.black)
if self.filtered_i_cube.weights is not None:
ax1.plot(self.filtered_i_cube.freqs[freqs_idx_zero_w] * 1e-6,
self.filtered_i_cube.data.var(axis=1)[freqs_idx_zero_w],
marker='+', ls='', c=psutil.black)
ax1.set_yscale('log')
ax1.set_xlabel('Freqs [MHz]')
ax1.set_ylabel('Variance [Unormalized]')
ax1.legend()
ax2.plot(
self.original_i_cube.freqs * 1e-6,
self.original_v_cube.make_diff_cube_interp().estimate_freqs_sefd(),
c=psutil.lorange, label='dV before')
ax2.plot(
self.original_i_cube.freqs * 1e-6,
self.original_i_cube.make_diff_cube_interp().estimate_freqs_sefd(), c=psutil.lblue,
label='dI before')
ax2.plot(
self.filtered_v_cube.freqs * 1e-6,
self.filtered_v_cube.make_diff_cube_interp().estimate_freqs_sefd(),
c=psutil.lorange, label='dV after')
ax2.plot(
self.filtered_i_cube.freqs * 1e-6,
self.filtered_i_cube.make_diff_cube_interp().estimate_freqs_sefd(), c=psutil.dblue,
label='dI after')
ax2.plot(self.original_i_cube.freqs[freqs_idx_filtered] * 1e-6,
self.original_v_cube.make_diff_cube_interp().estimate_freqs_sefd()[freqs_idx_filtered],
marker='x', ls='', c=psutil.black)
ax2.set_yscale('log')
ax2.set_xlabel('Freqs [MHz]')
ax2.set_ylabel('SEFD [Unormalized]')
ax2.legend()
if isinstance(self.original_i_cube, datacube.CartDataCube):
self.original_v_cube.weights.plot_uv(ax=ax3)
ax3.set_title('Weights')
self.original_v_cube.estimate_uv_sefd().plot_uv(ax=ax4)
ax4.set_title('SEFD V')
self.original_i_cube.make_diff_cube().estimate_uv_sefd().plot_uv(ax=ax5)
ax5.set_title('SEFD d_nu I')
self.original_v_cube.make_diff_cube().estimate_uv_sefd().plot_uv(ax=ax6)
ax6.set_title('SEFD d_nu V')
for ax in [ax3, ax4, ax5, ax6]:
ax.scatter(
self.original_i_cube.uu[uv_idx_filtered],
self.original_i_cube.vv[uv_idx_filtered], c=psutil.black, s=40, marker='x')
if self.filtered_i_cube.weights is not None:
ax.scatter(
self.filtered_i_cube.uu[uv_idx_zero_w],
self.filtered_i_cube.vv[uv_idx_zero_w], c=psutil.black, s=40,
marker='+')
elif isinstance(self.original_i_cube, sphcube.SphDataCube):
self.original_i_cube.plot_lm(ax=ax3, action_fct=np.mean)
ax3.set_title('Mean Stokes I')
self.original_v_cube.plot_lm(ax=ax4, action_fct=np.var, norm=LogNorm())
ax4.set_title('Var V')
self.original_i_cube.make_diff_cube().plot_lm(ax=ax5, action_fct=np.var, norm=LogNorm())
ax5.set_title('Var d_nu I')
self.original_v_cube.make_diff_cube().plot_lm(ax=ax6, action_fct=np.var, norm=LogNorm())
ax6.set_title('Var d_nu I')
for ax in [ax3, ax4, ax5, ax6]:
ax.scatter(
self.original_i_cube.ll[uv_idx_filtered],
self.original_i_cube.mm[uv_idx_filtered], c=psutil.black, s=40, marker='x')
else:
raise ValueError('Cube is not of a supported format')
fig.tight_layout()
return fig
[docs]
@staticmethod
def load(filename):
"""Build a flagger pipeline from an INI configuration file.
The ``[flagger]`` section defines an ordered, comma-separated
``pipeline``. Each named section must provide a ``type`` matching a
flagger class from this module. Its remaining values are passed to the
class as keyword arguments.
For example::
[flagger]
pipeline = bad_channels, noisy_baselines
[bad_channels]
type = FixedFreqsFlagger
freqs = 115.0-115.2, 121.5
action = filter
[noisy_baselines]
type = UVSigmaClipFlagger
nsigma = 5
action = zero_weight
Args:
filename: Path to the configuration file.
Returns:
A :class:`FlaggerRunner` containing the configured pipeline.
"""
config = configparser.RawConfigParser()
config.read(filename)
flagger_runner = FlaggerRunner()
assert config.has_option('flagger', 'pipeline'), f"Secton 'flagger' missing from {filename}"
for flagger_name in config.get('flagger', 'pipeline').split(','):
flagger_name = flagger_name.strip()
if len(flagger_name) > 0:
items = dict(config.items(flagger_name))
assert 'type' in items, f"No type found for '{flagger_name}'"
assert items['type'] in globals(), "Type '{}' incorrect".format(items['type'])
klass = globals()[items['type']]
assert inspect.isclass(klass), "Type '{}' incorrect".format(items['type'])
del items['type']
items['name'] = flagger_name
flagger = klass(**items)
flagger_runner.add(flagger)
return flagger_runner