Rewritten 2026-08-31. Everything before this was written around the tile pipeline and did not describe what we are actually building.
Find the RFI in a calibration observation using clustering, at native resolution. The clustering sorts pixels by kind, not by identity. Two narrowband transmitters at different frequencies are the same kind of thing and belong in the same bin. What must not happen is a bin holding interference and good data together, because that is the bin you cannot discard without paying for it.
Flagging is then discarding the bins that are interference, written into the FLAG column.
Then run bandpass and gain, apply the calibration, and show the corrected spectrum. If there is still a spike in it, it has not worked.
All three stages exist and run end to end on field 0. The corrected spectrum is NOT flat, so it has not worked.
Over-flagging is fine: the cluster mask discards 5.34% against tfcrop's 4.12%, a comparable budget. Under-flagging is the problem: the corrected spectrum keeps spikes at 1.17 to 1.30, 1.50 to 1.65 and 1.95 to 2.00 GHz, and at 1.17 to 1.30 the cluster mask is worse than tfcrop on nearly every antenna.
One cause, and it is the same one in stage 2 and stage 3. Interference that is both broad in frequency and always on has no detector. The spectral reference sits inside it and the temporal reference is blind to it. Stage 2 therefore leaves it in the sky bin, and stage 3 cannot remove what was never named. Solving that is the whole of the next piece of work.
The three figures to read, in order:
figures/named_clusters_field0.png for what the bins are,
figures/stage3_spectrum.png for the answer,
figures/stage3_spectrum_by_antenna.png for whether it is one antenna.
Stage 1, rebuilt 2026-08-31. See the progress log at the bottom for what changed and why.
scripts/fullres_pass.py reads the measurement set at full resolution and
gives four numbers per pixel, plus an exact zero mask.
Two are amplitude contrasts. ztime is how far the pixel sits from its own
channel measured over time. zspec is how far it sits from its
neighbouring channels, with the centre channel left out so a line cannot
lift the reference it is being measured against.
Two are the same references applied to phase, in radians. ptime is the
departure from that baseline's own average phase direction. pspec is the
departure from the circular mean of the neighbouring channels. Phase is
there because amplitude alone cannot separate an antenna that has slewed
off source from one with a weak receiver or a faint field. All three lose
amplitude. Only the slew loses phase coherence.
Nothing is thresholded. The values are kept so stage 2 can be re-run at any setting without touching the measurement set again.
Field 0 output is in /var/mnt/fast/vistok_fullres_v2/field0/, 16 spw,
37,200 rows by 2 polarisations by 128 channels each, 1.3 GB, two minutes to
produce. The older /var/mnt/fast/vistok_fullres/ is the previous
amplitude-only run off CORRECTED_DATA and should not be used.
scripts/fullres_check.py validates a stage 1 directory before anything is
built on it: shapes, non-finite values, that the amplitude contrasts are
centred on zero with robust width near one, and that the phase contrasts
stay inside plus or minus pi.
Alongside it there is a working chain to split a field to a scratch measurement set, write a mask into its FLAG column, calibrate, and measure what is left over.
Stage 2, the clustering, does not exist at native resolution. Nothing in the repo clusters the stage 1 output.
Stage 3, deciding how far to spread each flag once you know what kind of thing caused it, does not exist either.
Every clustering result this repo has ever quoted came from a different code path. That path cut the data into blocks of 8 integrations by 64 channels, computed statistics on each block, and clustered the blocks. The comparisons against tfcrop, the residual tables, the shuffle control, all of it was measured on those blocks.
That path never touches the full resolution output. They are two separate pipelines that do not meet. So none of those numbers say anything about whether clustering at native resolution works, and they should not be quoted as if they do.
Two further problems with the block features, in
src/visibility_tokenizer/features.py, which are worth remembering when
stage 2 gets designed:
Nothing in the feature list compares a channel against its neighbouring channels. A narrowband line raises the block average exactly the way a brighter baseline does, so brightness and interference arrive in the same number and cannot be separated.
The features are standardised against the median across all blocks. An interferer that is on in nearly every block at a given channel becomes the median, so it standardises to zero and disappears. Persistent transmitters are invisible to that method by construction.
Two numbers, both reported, neither assumed.
How much good data was discarded. A straight count of the flag fraction. Nothing else is needed to know how much was thrown away.
How much RFI was left behind. The corrected spectrum. Interference that survived the mask is still in the calibrated data as a spike, so the spectrum is the measurement. Flat is a pass. A spike is a failure and gets called one.
tfcrop gets run on the same field from the same starting point, and its flag fraction and its corrected spectrum sit next to ours.
No fake signals get injected into the data, and nothing is scored against a labelled truth set. The FLAG column in the measurement set is a previous flagger's output, not truth.
Broad and persistent interference has no detector, and it is the largest thing in this band. Spectral contrast cannot see it because its reference window sits inside the emitter. Temporal contrast cannot see it because it never turns off. Falloff with baseline length is right in principle, since a point source at the phase centre gives the same amplitude on every baseline and near-field interference does not, but as implemented it does not separate the bins: gap 0.57 against 3.08 for the spectral test.
The earlier worry, that clustering every pixel would spend all its bins splitting the noise, did not happen. With k = 120 the noise fragmented into about a hundred bins and left centroids free for everything else, which is what it was supposed to do.
Sixteen.
fullres_pass.pythe per pixel pass. Everything rests on this.fullres_check.pyvalidate a stage 1 directory before building on it.prepare_scratch.pysplit a field to a scratch measurement set, clear its flags.apply_mask.pywrite a mask into a scratch FLAG column.capture_flags.pyread a flag column back out lined up with the stage 1 rows, so a comparison is on the same pixels.calibrate_compare.pyflag, bandpass, apply, then measure the residual.spectrum_before_after.pythe corrected spectrum figure.grow_flags.pyspread a mask outward the way tfcrop spreads its own, so a comparison is fair.flags_per_integration.pywhere a mask's flags fall in time.
Stage 2 and stage 3, in run order:
stage2_features.pythirteen numbers per pixel, with context.stage2_cluster.pycluster every pixel, no pre-selection.stage2_name.pyname the bins from physics; prints the gap at each cut.plot_named_clusters.pythe stage 2 test figure.stage3_mask.pynamed bins to a flag mask.stage3_calibrate.pyflag, setjy, delay, bandpass, gain, applycal, measure.plot_stage3_spectra.pythe corrected spectrum, aggregate and per antenna.
Run them as:
pixi run -e ms python scripts/prepare_scratch.py --field 0 --reuse
pixi run -e ms python scripts/fullres_pass.py \
/var/mnt/fast/vistok_scratch/start_field0.ms \
--out-dir /var/mnt/fast/vistok_fullres_v2 --tag field0
pixi run -e ms python scripts/fullres_check.py \
/var/mnt/fast/vistok_fullres_v2/field0
pixi run -e ms python scripts/stage2_features.py \
/var/mnt/fast/vistok_fullres_v2/field0 \
--out /var/mnt/fast/vistok_stage2/field0
pixi run -e ms python scripts/stage2_cluster.py \
/var/mnt/fast/vistok_stage2/field0 --k 120
pixi run -e ms python scripts/stage2_name.py \
/var/mnt/fast/vistok_stage2/field0
pixi run -e ms python scripts/plot_named_clusters.py \
/var/mnt/fast/vistok_scratch/start_field0.ms \
/var/mnt/fast/vistok_stage2/field0 \
--fig figures/named_clusters_field0.png
pixi run -e ms python scripts/stage3_mask.py \
/var/mnt/fast/vistok_stage2/field0 \
--out /var/mnt/fast/vistok_stage3/field0
pixi run -e ms python scripts/stage3_calibrate.py \
/var/mnt/fast/vistok_scratch/start_field0.ms --field 0 \
--cluster-mask /var/mnt/fast/vistok_stage3/field0 \
--out out/stage3/field0
pixi run -e ms python scripts/plot_stage3_spectra.py out/stage3/field0
About twenty minutes end to end, most of it in the three calibrations.
Known problem in spectrum_before_after.py: each mask solves its own
bandpass on data that still has interference in it, then divides by it. In
a contaminated channel that solution is wrong and dividing makes the spike
larger. That panel is measuring the solve, not the mask. It needs one fixed
bandpass applied to every mask.
Science measurement set, never write to it:
/var/mnt/Data/Data/22A-256/start.ms
Scratch, one per field, flags cleared:
/var/mnt/fast/vistok_scratch/start_field{0..4}.ms
Stage 1 output, current: /var/mnt/fast/vistok_fullres_v2/field0/
Stage 1 output, superseded (amplitude only, off CORRECTED_DATA), do not
use: /var/mnt/fast/vistok_fullres/field{0..4}/
Fields: 0 is 3C138, 1 is 3C48, 2 is J2355+4950, 3 is J0102+5824, 4 is pa30, which is the target and is 8.4 times fainter. Sixteen spectral windows, 128 channels each, 0.994 to 2.006 GHz.
Three faults were found in stage 1 by reading the source, all of which would have made stage 2 impossible as specified.
fullres_pass.py:150 defaulted to reading CORRECTED_DATA. That default was
set for a comparison against tfcrop, which runs on corrected data. It is
wrong for this work: we find the RFI first and calibrate afterwards, so
reading CORRECTED_DATA means hunting for interference in data a previous
solve has already reshaped. Now defaults to DATA.
fullres_pass.py:188 took the absolute value and discarded phase entirely.
Fixed: the complex values are carried through and two phase contrasts are
written.
Correlations. Only RR and LL are kept, correlations [0, ncorr-1] of four.
RL and LR are dropped. This is left as it was, deliberately. Highly
polarised interference would show in the cross hands, but adding them
doubles the problem before the parallel hands have been shown to work.
Revisit once stage 2 works.
The first attempt referenced phase per (polarisation, channel) across the
whole array, mirroring how the amplitude reference works. The check script
caught it. ptime came out with a robust width of 2.234 radians, and a
uniformly random phase gives 2.33, so the statistic was carrying nothing.
The reason is that amplitude and phase do not behave the same way under averaging across baselines. Every baseline sees a comparable amplitude, so a per-channel amplitude median across the array is meaningful. Every baseline has its own delay and fringe, so averaging their unit phase vectors cancels to nothing and the reference direction is arbitrary.
Referenced per baseline instead, the width drops to 1.274 radians on spw 0 and to 0.03 to 0.05 radians on most spw.
A per (baseline, scan) reference would be tighter still, and was rejected.
A scan that is entirely slewing would set its own reference and the slew
would disappear, which is the same structural blindness that ztime has
against an always-on emitter. The drift-immune version is a local-in-time
phase contrast, and stage 2 can compute that from the stored arrays plus
the time and antenna columns without re-reading the measurement set.
Produced by scripts/fullres_pass.py, checked by
scripts/fullres_check.py, written to
/var/mnt/fast/vistok_fullres_v2/field0/check.json.
Amplitude contrasts behave as they should everywhere: ztime median
between -0.05 and +0.01 with robust width 0.94 to 1.01, zspec median
between -0.03 and 0.00 with width 0.63 to 1.00. A deviation expressed in
robust sigma should sit at zero with width one, and it does.
Phase coherence varies a lot by spw, which is itself informative:
| spw | ptime width, rad | pspec width, rad |
|---|---|---|
| 0 | 1.274 | 1.749 |
| 1 | 0.485 | 0.553 |
| 2 | 0.266 | 0.296 |
| 3 | 0.042 | 0.038 |
| 4 | 0.040 | 0.036 |
| 5 | 0.041 | 0.038 |
| 6 | 0.052 | 0.048 |
| 7 | 0.226 | 0.252 |
| 8 | 0.044 | 0.036 |
| 9 | 0.059 | 0.052 |
| 10 | 0.101 | 0.098 |
| 11 | 0.288 | 0.295 |
| 12 | 0.201 | 0.219 |
| 13 | 0.043 | 0.037 |
| 14 | 0.037 | 0.030 |
| 15 | 0.035 | 0.031 |
Most of the band is phase coherent to better than a twentieth of a radian, which means a pixel that departs from its baseline reference by even a tenth of a radian is a real outlier there. That is a sharp discriminator and it is the reason phase was worth adding.
Spw 0 is an order of magnitude worse than the rest and spw 1, 2, 7, 11 and 12 are intermediate. Not yet explained. It could be genuine interference, or the low end of the band being less sensitive. It is not being treated as either until stage 2 has something to say about it.
Zeros run at 2.02% in every spw except spw 1 at 3.33% and spw 4 at 3.00%. The consistency across spw suggests a structural cause rather than an interference one, which is what stage 2 should find and name.
Nothing from stage 1. Stage 2 and stage 3 are logged below.
Read figures/named_clusters_field0.png before anything below. It is the
test and it is generated by scripts/plot_named_clusters.py.
Three scripts, run in order:
stage2_features.py gives every pixel thirteen numbers. The four stage-1
contrasts, the zero mask, and eight more that place the pixel in context:
what its channel neighbourhood is doing, what its own baseline was doing a
few integrations either side, what its two antennas were doing across the
array at that moment, what the whole array was doing at that moment in both
amplitude and phase, the baseline length, and how amplitude falls off with
baseline length. 0.5 GB per spw, three seconds per spw.
stage2_cluster.py clusters every pixel. All 152,371,200 of them for field
0, no pre-selection, k = 120, three minutes. Two metric choices, both
documented in the file: robust standardisation so no feature dominates by
its units, and arcsinh so a handful of extreme pixels do not take every
centroid. Neither cuts anywhere. k is large on purpose: the noise is
allowed to fragment into as many bins as it likes, because fragmenting the
noise costs nothing and absorbing the interference would cost everything.
stage2_name.py names the bins from their statistics, using rules written
from physics with one shared number, sigma = 3.
| name | bins | pixels | share | keep |
|---|---|---|---|---|
| sky | 105 | 145,494,425 | 95.487% | keep |
| zeros | 8 | 2,676,224 | 1.756% | discard |
| low gain | 2 | 929,584 | 0.610% | discard |
| slewing | 1 | 1,032,297 | 0.677% | discard |
| shadowing | 0 | 0 | 0.000% | discard |
| RFI broad continuous | 2 | 907,947 | 0.596% | discard |
| RFI broad bursty | 2 | 1,330,723 | 0.873% | discard |
| RFI narrow continuous | 0 | 0 | 0.000% | discard |
| RFI narrow bursty | 0 | 0 | 0.000% | discard |
Discarded 4.513%. Kept 95.487%.
Zeros come out clean. Eight bins, 1.756%, exactly zero everywhere in the band, and nothing else lands in them.
There is a genuine interference bin. RFI broad continuous sits well above
the band median at 1.15 to 1.35, 1.5 to 1.7 and 1.95 to 2.0 GHz, and drops
away between them. That is what an interference bin is supposed to look
like.
low gain sits below the band median across the whole band, which is the
right shape for an instrument fault rather than a frequency-specific one.
The zspec rule separates cleanly. The bins sit 2.42 below the cut and
5.50 above it, a gap of 3.08, so where exactly the cut was put does not
matter.
The largest emitter in the band is still being kept. sky carries the full
1.53 to 1.56 GHz spike and the 1.62 GHz spike at the same height as the
band median. That is under-flagging, and it is the single most important
thing wrong.
The reason is structural and worth stating precisely. zspec references a
pixel against its neighbours plus and minus eight channels, about 4 MHz.
The 1.5 to 1.6 GHz feature is roughly 100 MHz wide, so its own reference
sits entirely inside it and the contrast is near zero. ztime references a
pixel against its own channel over time and is blind to anything always on.
An emitter that is both broad and persistent escapes both, and on this data
that describes the strongest thing in the band.
No narrowband bins at all, and no shadowing bin.
The four-way interference split does not separate. The narrow-against-broad cut has a gap of 0.31 and the continuous-against-bursty cut a gap of 0.24, against 3.08 for the interference test itself. Those two splits are landing inside a continuum, so the labels they produce are not meaningful yet even where they fire.
The context features were first computed as raw levels. Phase coherence varies by a factor of thirty across the band, from 0.03 radians in spw 14 to 1.27 in spw 0, so one scale over the whole band put every pixel in the least coherent windows past a cut that was reasonable in the coherent ones. The run named 35% of the field slewing and found no interference at all. Fixed by referencing each context feature against its own channel over time, so it asks whether this moment is unusual for this channel rather than whether the number is large.
Then uv_slope was added to catch the broad persistent case, and passed
through that same per-channel normalisation, which destroyed it. A
persistent emitter has the same slope at every timestamp in its channel, so
normalising over time subtracted exactly the signal the feature existed to
find. With the normalisation in place it named 22% of the field
interference in bins whose spectra tracked the band median across the whole
band, which is ordinary data being thrown away.
uv_slope is computed and the clustering uses it, but it is NOT a naming
rule by default. Its gap either side of the cut is 0.57 against 3.08 for
zspec, so it does not separate the bins and where the cut lands is
arbitrary. --use-uv-slope reproduces that run: it discards 32.377%,
mostly into a 22% bin that looks like sky. Turning it off is not hiding the
result, it is applying a criterion fixed in advance, that a rule is only
used when the bins fall either side of it with a gap.
The slewing rule originally tested phase coherence alone. Its bin then sat ABOVE the band median at exactly the interference frequencies, because interference scrambles phase too. Slewing must lose amplitude as well as coherence, since coming off source takes the signal away while interference adds power. With the amplitude condition the slewing bin sits below the band median across the whole band with no spikes, which is the right shape.
Broad and persistent interference has no detector. Spectral contrast cannot see it because its reference is inside it. Temporal contrast cannot see it because it is always on. Baseline-length falloff is the right idea physically, a point source at the phase centre gives the same amplitude on every baseline and near-field interference does not, but as implemented it does not separate the bins.
Until that is solved the corrected spectrum will keep its largest spike, and stage 2 has not met the bar.
Read figures/stage3_spectrum.png and
figures/stage3_spectrum_by_antenna.png.
Three scripts:
stage3_mask.py turns the named bins into a per-spw boolean mask, keeping
the pixels named sky and discarding the rest. No threshold, no growth.
stage3_calibrate.py runs the full chain per mask from a completely
unflagged scratch copy: unflag, apply the mask, setjy, delay, bandpass,
gain, applycal, then measure the median spectrum for the array and for each
antenna. Each mask is calibrated by itself, which is the right comparison
when the question is what calibration a mask can produce.
plot_stage3_spectra.py draws it. Three curves and no more: uncorrected,
tfcrop, clusters.
Over-flagging, a straight count from flagdata(mode='summary'):
| mask | flagged |
|---|---|
| none | 0.000% |
| tfcrop | 4.123% |
| clusters | 5.338% |
Comparable budget. tfcrop discards 4.12%, the cluster mask 5.34%. This is the one place the work has clearly improved: the old tile-based mask discarded about nine times what tfcrop did, and at native resolution the two are within thirty percent of each other.
Under-flagging, read off the corrected spectrum: FAILED, for both masks.
Neither corrected spectrum is flat. Both keep large spikes at 1.17 to 1.30, 1.50 to 1.65 and 1.95 to 2.00 GHz, reaching 100 or more against a calibrated continuum of about 8 to 10.
The spikes are not an artefact of the spectral window edges. That was checked: the outer eight channels of every window are blanked in the figure and the spikes survive it.
At 1.17 to 1.30 GHz the cluster mask leaves a spike reaching about 100 where tfcrop leaves about 20. The per-antenna figure shows this on nearly every antenna, not on a few. At 1.50 to 1.65 they are comparable. At 1.95 to 2.00 the cluster mask is again higher.
So on the measurement that matters, the cluster mask is not better than tfcrop, and over a wide part of the band it is worse. That is the state of it. It costs slightly more data and it leaves more interference behind.
The stage 2 failure and the stage 3 failure are the same failure. The
emitters that survive into the corrected spectrum are at 1.17 to 1.30 and
1.50 to 1.65 GHz, and those are precisely the regions where the stage 2
figure shows the sky bin still carrying the emitter at full height. The
mask cannot remove what the clustering did not name.
Broad and persistent interference has no detector. Until it does, the rest of the chain is sound but is being fed an incomplete mask.
In the per-antenna figure, ea10 and ea21 are much cleaner than the other twenty three, with almost no spikes outside 1.5 to 1.65 GHz. Not investigated. It could be their position in the array, or it could be that those antennas simply have more of their data flagged. Not measured either way, so nothing is being claimed from it.