Configuration File

TABASCAL makes use of a configuration file to fully define the model used, what satellites to include and what inference to do. The configuration file is in YAML format using only the most basic constructs and types. The configuration file is broken up into multiple sections. Each section will be described below.

Your file is merged onto the base configuration that ships with TABASCAL, key by key, so only the keys you want to change have to appear in it — everything else takes the documented default. A section header written with nothing under it is therefore inert: rfi: alone, like rfi: null, leaves the whole rfi section at its defaults rather than emptying it. Setting an individual key to null is a different thing entirely — it is a value, and the sections below say what each one means (rfi.min_elevation: null disables the elevation mask, data.noise: null reads the noise from the measurement set). Which of the two you get follows the default the key has: an empty value is ignored only where the default is a section of further keys, so a key whose default is a single value or a list — satellites.norad_ids:, written with nothing after it — is null like any other. Leaving a section empty, or setting it to null, never removes its defaults.

Model

The forward model used by TABASCAL is modular and configurable directly in the model section of the configuration file. An example of this section is shown below.

model:
  name: Custom
  components:
    # RFI Phase
    - trajectory:Orbit
    - trajectory:PhaseCalculationRFI
    # RFI Signal
    - rfi_signal:ComplexRFIVarAnt
    # RFI Visibilty
    - rfi_vis:RiemannVis
    # Astronomical Visibility
    - ast_vis:GPVisAst
    # Gains
    - gains:UnitaryGains

The components should be given in order of dependency. For example, trajectory:Orbit is specified before trajectory:PhaseCalculationRFI because the later depends on the output of the former. Each component is a class which defines the parameters (if any), their initialisation, their prior distribution, and its own forward model for the component. The component modules are located in tabascal/components/. The model component given in the configuration file should use the module name and then the class name. For example the component trajectory:Orbit is a class named Orbit that resides in the module file tabascal/components/trajectory.py

That order is checked when the model is assembled, before anything is computed. Each component declares the state keys it reads and the keys it writes, and a list that leaves out a component — or holds the right ones in the wrong order — is rejected by name, saying which key is missing, what produces it, and whether that producer is absent or merely listed too late.

Precision

TABASCAL runs in single precision (fp32) by default. It is set in the same model section:

model:
  precision: single   # default; or "double" for fp64

single (default). Halves device-memory use (~2×), which raises the largest problem size that fits on a GPU. On GPUs with first-class fp64 (e.g. Hopper/GH200) it is not faster in wall-clock — the win is memory capacity, not speed.
double. Required by some components, and recommended when fitting satellite trajectories, as the differentiable orbit models need fp64 accuracy.

The following components run in double precision only and raise a clear error under single, so set model.precision: double to use them:

trajectory:PhaseCalculationRFI
trajectory:NoDragOrbit
trajectory:Orbit

Both rfi_vis kernels (RiemannVis and the FFI RiemannVisFFI, see RFI-visibility kernels) run in either precision, as do the astronomical Gaussian process ast_vis:GPVisAst and both gain components, gains:ConstGains and gains:UnitaryGains.

Data

The data section of the configuration file includes only a few options to select the data to use. An exhaustive example is given below.

data:
  out_dir:
  truth_zarr:
  ms_path: path/to/ms_file.ms
  data_col: DATA
  freq:
  corr: xx
  noise:
  gain_table:
  save_rfi_per_sat: false

out_dir: Where the run writes: plots/ and results/, the latter holding the used_orbits file. Also available as the -od flag, which takes precedence. Defaults to the directory ms_path names, so a run on real data needs neither — products land beside the visibilities. Given without an ms_path it is read as a tab-sim simulation directory, and the MS and truth zarr are looked for inside it at <out_dir>/<out_dir name>.ms and .zarr, which is the layout sim-vis writes. Naming neither an MS nor a directory is the one combination that cannot be resolved. Note that -od moves the outputs but not the MS when ms_path names one; -ms is how a run is moved onto other visibilities. (Renamed from sim_dir, which meant the inputs and the outputs at once — see #207.)
truth_zarr: The tab-sim simulation truth, read by ast.init: truth, rfi.init: truth and plots.truth. Defaults to <MS directory>/<MS name>.zarr — beside the visibilities, which is where sim-vis writes the two together — so a simulation needs it only when its truth has been moved away from its MS. A real observation has none: with any of those three options set, a missing truth stops the run before it fits, naming which option asked for it; with none of them set it is a printed “No tab-sim truth available” and the run proceeds.
ms_path: Path to the Measurement Set (MS) to run on. This can also be given at runtime of tabascal with the -ms flag, which takes precedence over the config. Leave both unset and the MS is looked for at <out_dir>/<out_dir name>.ms, the layout sim-vis writes — so a simulation can name the directory instead. Real data wants this key, since an MS taken from a telescope has no relationship to the name of the directory it happens to sit in.
data_col: The data column within the MS file to use as the observed data. Default is DATA but can be any column that exists in the MS file.
corr: This is the correlation product to run on, default xx. It is matched against the MS’s POLARIZATION::CORR_TYPE by identity, not by position, so it names the correlation you want rather than an axis index: yy selects YY whether the MS holds all four correlations or only that one. Linear (xx, xy, yx, yy), circular (rr, rl, lr, ll) and Stokes (i, q, u, v) names are accepted. Requesting a correlation the MS does not hold is an error naming what it does hold.
freq: A single frequency in Hz to read instead of the whole band; the nearest channel is used. The request must fall inside the band — a frequency more than half a channel beyond the nearest centre is an error rather than a silent read of the edge channel, since argmin always returns a channel and a units slip would otherwise pass unnoticed. null (the default) reads every channel. SIGMA_SPECTRUM is narrowed to the same channel, so the noise cannot come back on a channel the visibilities did not.
skip_ms_write: Leave the run’s results in the zarr and skip copying them into the Measurement Set. Default false. A fit writes its products twice — once into results/, and again into the MS’s data columns — and on a large observation the second copy is a substantial part of a run that is otherwise finished. With this set the run stops after the zarr, and tab2MS copies the same products across later, so nothing is lost and the choice is deferred rather than made in advance. Useful when the MS is on slow or shared storage, when several fits are being compared and only one of them will be kept, or when the MS is read-only. Also available as the --skip-ms-write flag, which turns it on whatever the config says; the flag only enables, so a config with it set true is not overridden by leaving the flag off.
row_chunk: The number of MS rows read and written per chunk, or null (the default) to leave dask-ms’s own chunking alone. It is a throughput and memory knob with no effect on results: a larger chunk means fewer, bigger reads, which is usually faster until the chunk stops fitting comfortably in host memory. The right value depends on the machine and the file rather than on the observation, which is why there is no default worth shipping — leave it unset unless MS I/O is measurably costing you. Must be a positive whole number if given.

noise: The per-visibility noise in Jy. Leave it null to use the MS’s own noise columns, which are read per baseline, and per channel where the MS resolves them that far, rather than averaged to one number: the antennas of a real array differ in sensitivity and a bandpass is not flat, so a single value mis-weights every visibility. On a real low-frequency array the per-baseline SIGMA has been measured to span a factor of ~30, so a scalar under-weights the quietest baselines by up to ~200x. It matters most when fitting gains, because the per-antenna noise correlates with the per-antenna gain (on the same data, sigma_a ~ amplitude_a^0.76, R = 0.96) — a uniform-noise likelihood cannot tell a loud antenna from a noisy one, so the fitted gain absorbs the noise structure.

With noise: null the noise is resolved from the MS, most specific column first:

column	shape read	constant in time	varying in time
`SIGMA_SPECTRUM`	`(row, chan, corr)`	`(n_bl, n_freq)` — per baseline and channel	`(n_bl, n_freq, n_time)`
`SIGMA`	`(row, corr)`	`(n_bl,)` — per baseline	`(n_bl, 1, n_time)` — per baseline and timestep

The column is checked against the MS’s whole channel axis before any single-channel selection is applied: a SIGMA_SPECTRUM that covers a different number of channels from DATA is a disagreement about the observation, and narrowing both to one channel would let it pass without settling which is right.

The time axis is kept only when it carries something. Each column is read cell by cell — per baseline for SIGMA, per (baseline, channel) for SIGMA_SPECTRUM — and a cell whose rows are all bit-identical is one measurement written into every row: it collapses to that value, so an MS with a uniform SIGMA weights exactly as it always did. If any cell’s rows differ, the column is saying the noise changed over the observation, and the whole column is kept time-resolved. There is no tolerance and no threshold: an MS that writes a constant noise writes the identical value in every row, and a column that varies is taken at face value, because nothing here can tell a corrupted row from a timestep on which the noise really was different — both are a positive, finite number the MS wrote down. One odd-but-finite row is therefore enough to make a column vary, and it is kept rather than median-ed away. Validity masking applies either way: a non-positive or non-finite entry is never a measurement. The read prints one line when it finds time variation, saying the noise is being kept time-resolved. A time-resolved SIGMA is (n_bl, 1, n_time) and never (n_bl, n_time) — whenever the observation has as many channels as timesteps, nothing in that shape says which axis is which, and every consumer would weight the visibilities by the wrong one in silence.

Where a column does collapse, it collapses with a median per cell, over the rows that carry a value — non-positive and non-finite rows are dropped before the median, not after it, so a few invalid rows cannot cost a baseline the rows that measured it properly. Cells left with no valid value at all — dead baselines, flagged channels — take the median of the valid cells rather than a zero that would divide the likelihood by nothing, with a warning saying how many. In the time-resolved case the same fill runs one axis further out: an entry that measured nothing takes its own cell’s median over the timesteps that did, and only a cell that measured nothing at all falls back on the median over the cells that did. Such cells are normally flagged, and a flagged cell is dropped from the likelihood whatever noise it was given; the fill itself does not read the flags, so an unflagged cell with no valid noise is included at the median. Take the warning as saying which cells to check rather than as a guarantee they carry no weight. SIGMA_SPECTRUM is optional in the MS format: if it is absent, unreadable, covers a different set of channels from the ones being read, or holds no positive finite value anywhere — a column that was never filled in — the read falls through to SIGMA with a warning. If neither column is usable nothing is invented — a noise scale invented on the spot would silently re-weight the entire fit — so the read warns and leaves the noise unset. Setting noise is then what supplies it, and is how such an MS is run; left null, the run stops with an error naming this option.

A scalar may still be given, applying to every visibility. 0 is rejected rather than treated as “no override” — zero noise divides the likelihood by nothing, so whichever was meant, say it. true/false are rejected too: YAML booleans reach float() as 1.0 and 0.0, and a run on a uniform 1 Jy noise nobody asked for is plausible enough to go unnoticed. A path to an .npz is also accepted, for a noise measured out of band; it must carry one of, most specific first:

key	shape	meaning
`sigma_bl_freq`	`(n_bl, n_freq)`	per-baseline, per-channel noise, used as given
`sigma_bl`	`(n_bl,)`	per-baseline noise, used as given
`s_ant`	`(n_ant,)`	per-antenna noise, combined as `sqrt(s_p^2 + s_q^2) / sqrt(2)`

sigma_bl and s_ant must be one-dimensional, s_ant must cover every antenna the observation’s baselines are formed from, and the values read must be real, positive and finite: every entry of sigma_bl_freq or sigma_bl, and, for s_ant, every antenna this observation actually correlates. A file of per-antenna noise may legitimately cover a whole array, and an entry for an antenna these baselines never use cannot mis-weight anything, so those entries are deliberately not policed. A complex array is rejected rather than read as its real part, and so are boolean and string arrays — only an integer or floating-point array is read as a noise, since astype(float) would otherwise turn a file of flags into a uniform 1 Jy and parse a file of text into whatever the strings spell. An override is used exactly as given, so it is not repaired the way an MS column is: a file carrying an entry that is not a noise is rejected, naming the key and how many entries offend, rather than having part of the user’s own answer filled in for them. The median fill above applies to the MS columns only. No override is time-resolved — there is no key for a time axis, and a scalar, per-baseline or per-(baseline, channel) override applies to every timestep — so an MS whose noise genuinely varies over the observation is best read from its own columns.

The values are checked again after conversion to the precision the run works in (model.precision). Under single precision a value like 1e-50 underflows to zero — which would divide the likelihood by nothing — and 1e40 overflows to infinity, so a file that is valid as written but not at the run’s precision is an error naming the dtype, not a silent re-weighting.

noise is in the frame of data_col, i.e. before gain_table is applied. With a gain table, leave it null and let the MS’s SIGMA carry the noise rather than pre-scaling it by hand — the table is divided out of the noise as well as the data, so a value given here would be scaled a second time.

gain_table: Path to a CASA calibration table — from gaincal, flux-calibrate, or anything else that writes one — or an ordered list of them. The gains are divided out of the visibilities and the noise once, when the MS is read:
```
vis_obs   = DATA  / (g_p conj(g_q))
sigma_cal = SIGMA / |g_p conj(g_q)|
```
so everything downstream — priors, the RFI and astronomical models, chi², and the results written back to the MS — lives in one frame, the calibrated one, and the gains are applied once instead of on every forward pass. Carrying the noise with the data is the point of using a table rather than scaling the data by hand: get it wrong and chi² is off by |g|². Scale ast.gp_cov.std by k alongside it if you are setting a width by hand rather than measuring it.

The division happens in memory, so data_col on disk is untouched and still raw. The results writer therefore removes the same layer again when it writes the TAB_* columns and the weights; the run hands it this list automatically, and tab2MS takes it as -gt — see Output.

A table is solved on whatever (frequency, time) grid the calibrator chose, so it is placed on the observation’s grid first. Matching is by value, not by index — a time within 1 ms and a channel within 1e-6 of the band centre frequency count as the same sample — which is what lets a table solved on a master MS apply to a subset carved out of it, in whatever channel order the subset was written. The times matched are the MS’s own TIME column on the scale it declares (times_mjd), not the UTC-normalised times_jd, because a caltable’s TIME is a copy of the MS’s: declared frame to declared frame is exact, where the UTC coordinate of a TAI-declared MS is 37 s away.

Where the grids do not coincide the gains are interpolated linearly in amplitude and unwrapped phase, never in real and imaginary parts: two unit gains 60° apart average to |g| = 0.87 in real/imag, so the data would be divided by a gain no antenna ever had and the flux scale would move by 13 %. The phase is unwrapped in two dimensions — along frequency within each timestep, and then the timesteps onto one another by whole turns — and stays a real surface through both interpolation stages, exp(i·phase) being applied once at the end. A B table winds across the band through a residual delay and around the ±π branch cut in time, so interpolating the stored angles across such a step averages the two sides of the cut into a gain pointing the wrong way; rebuilding a complex gain between the two stages loses the branch just as badly, and tears the band by a whole turn at whichever channel crossed the cut. Beyond the solved range the edge value is held: a table that does not reach the start of the observation calibrates it with the earliest solution it does have.

The table has to sample the phase below half a turn, and nothing here can check that it did. Unwrapping recovers a phase only where the solutions sample it: a genuine change of more than π between two adjacent solved samples — between two solved channels, or between two solution intervals of a phase slewing faster than the calibration cadence follows — is simply not in the table, because complex solutions carry the phase modulo 2π. Such a step aliases to the shorter branch and is taken as such: a true 0 → 1.5π evolution between two solutions interpolates to −0.25π half way, not +0.75π. This is the same assumption np.unwrap makes, applied coherently across the band rather than channel by channel, and it is undetectable from the data — it is a requirement on the calibration that produced the table (solve finely enough in both axes), not something the interpolation can warn about. A step of exactly π is the boundary case, and is left as written rather than turned into a half turn of the opposite sign.

A flagged or zero solution is an absence, not a value: the interpolation bridges across it from the solutions either side, exactly as it bridges a coordinate the table never sampled. Only an antenna with no valid solution anywhere has nothing to interpolate from — its gain is 1 and its visibilities are flagged, even when flags: false, since a visibility nobody calibrated is not data.

With a list, each table is interpolated onto the observation’s grid and only then are they composed, g_total = Π gᵢ in the order given. The two orders disagree: two amplitudes ramping 1 → 3 give 2 × 2 = 4 half way when each is interpolated first, and (1 + 9) / 2 = 5 when the product is interpolated, which is an artefact of fitting a quadratic with a straight line.

Each table prints one coverage line as it is read, giving the fractions of the observation’s samples whose gain was taken exactly from a solution, interpolated between solutions, held from an edge, or left unsolved — a table that turns out to cover the observation mostly by extrapolation says so rather than being applied in silence. Each sample is classified by the support it was actually built from, so a table whose solutions run along a diagonal reports the edge-holds it really performed rather than the rectangle its two axes span. The coverage line is also what catches a table whose TIME was written on some other unit or scale, which is assumed rather than read: it would report no exact cover at all.
save_rfi_per_sat: Also store the fitted RFI visibility split per satellite in the results .zarr, as rfi_vis_src (sample, src, bl, freq, time) with a norad_id coordinate naming the satellite behind each src. Default false.

It is a diagnostic for astronomical signal leaking into the RFI model: a genuine satellite is a clean streak in exactly one per-source image, while a feature that appears in several is sky flux the model has split across satellites — which reduced chi² is blind to, since the split costs it nothing. Write the sources into MS columns with tabascal rfi-per-sat and image one at a time.

Off by default because it is not free. The stored array is n_rfi times the size of rfi_vis, and filling it costs n_rfi extra evaluations of the run’s own RFI-visibility op — one per satellite, each over the whole source axis with the other satellites’ amplitudes held at zero, which is what keeps the evaluation inside the RFI-axis sharding rather than gathering the fine grids onto one device. That is a few forward passes’ worth of work at the end of a fit, not a second fit, but it is n_rfi × the storage forever.

The multiplier is on disk only. Each (sample, satellite) block is written into the results zarr as it is evaluated — the store is created before the first evaluation and filled one chunk at a time — so the decomposition is never assembled in memory. What it adds to the writing process is a small multiple of one block, sizeof(rfi_vis) / n_sample, whatever n_rfi is: the block on the device, its copy on the host, and whatever the chunk write buffers. On disk it is chunked one satellite per chunk, so reading or imaging a single source does not pull the rest of them in either.

Padded sources are not stored: under sharding the satellite list is padded to a multiple of the device count with dark dummies, and only the real satellites get a src slice. The sources sum back to rfi_vis exactly in exact arithmetic; in floating point the split re-associates the op’s single reduction over (source, integration sample), which is round-off on fitted grids (~2e-16 relative in double, ~6e-8 in single) but is bounded by the fine-grid terms rather than by the coarse visibilities — see the per-satellite columns for the case where those two differ by everything.

Plots

The plots section defines which plots of predictions to create. These produced plots are saved in the plots/ directory alongside the MS file TABASCAL was run on. These are defined as follows

plots:
  init: True
  truth: False
  prior: True
  prior_samples: 100

The above configuration will plot the predictions for the initial parameters and the prior distributions. To plot the prediction from the prior, 100 samples will be drawn from the prior and then pushed through the forward model to get a predicition. The truth plots can only be included if TABASCAL is being run on a simulation created with sim-vis.

Inference

The inference section defines the type of inference that will be done. An example is given below.

inference:
  opt: True

opt: Optimisation will be done to find the maximum a posteriori (MAP) point. It is the only inference this section offers.

Optimisation

The opt section gives the parameters for the optimiser. An example is given below.

opt:
  epsilon: 1e-1
  max_iter: 10000
  dual_run: True
  trace_path: null

epsilon: This is the optimiser step size
max_iter: This is the number of iterations to run the optimiser for in each optimisation run.
dual_run: This determines whether the optimiser will be run a second time starting from where the previous run left off but with an epsilon that is 10x smaller.
trace_path: Path of an .npz file to record the optimiser’s convergence to, or null (the default) to record nothing. See Optimiser trace.

Optimiser trace

Two models cannot be compared on their loss curves. The loss is a negative log joint, so its prior term scales with the latent dimension of whichever parameterisation is running — gains:UnitaryGains, which fits no gain at all, and gains:ConstGains, which fits \(2 n_\text{ant} - 2\) parameters, do not put their losses on a common scale. And loss per iteration hides the cost of an iteration, so a model that converges in fewer but more expensive steps looks better than it is.

Setting trace_path records, once per optimiser iteration, the quantities that can be compared: the wall-clock time the iteration finished, and the metrics fixed by the data rather than by the parameterisation. The metrics are read out of the same forward pass as the loss, so they cost a few elementwise reductions rather than a second evaluation of the model, and the run is otherwise unchanged. With trace_path left null nothing is recorded and the optimiser takes the same compiled path it always did.

Set the TAB_LOSS_TRACE environment variable to override trace_path for a single run, for tracing a config that should not be edited. Under multiple processes every process traces — they all evaluate the same model and must run the same program — and one file is written, by the first process.

The file is written once, at the end of the run, and holds one array per key, each of length max_iter (or 2 * max_iter with dual_run) with one entry per iteration:

Key	Present	Meaning
`loss`	always	The optimiser’s own loss, the negative log joint divided by `vis_obs.size`
`time_s`	always	Seconds from the start of the first iteration to the end of this one, measured after the device sync that reads the loss — so it bounds work completed, not work dispatched
`chi2`	always	Reduced chi-squared of the observed visibilities, flag-masked and weighted by the resolved noise, exactly as the value printed at init and opt
`vis_ast_nrmse`	with truth	RMSE of the recovered astronomical visibilities against the simulation truth, over the representative noise
`vis_rfi_nrmse`	with truth	RMSE of the recovered RFI visibilities against the simulation truth, over the representative noise

The nrmse keys need a truth to score against, so they appear only on a dataset simulated with sim-vis. Every entry is recorded at the parameters that produced that iteration’s gradient, i.e. before its update — so the first entry is the value at the initialisation, and the last is one update behind the reported optimum.

Astronomical Signal

The ast section defines the prior distrbution, intialisation and forward model parameters for the astronomical signal. An example is given below.

ast:
  init: prior
  mean: 0
  freq_pad_factor: 2.0
  time_pad_factor: 2.0
  baseline_block_size: auto
  cutoff: 1e-6
  gp_cov:
    std: data
    corr_freq: null
    fov_deg: 5
    gammas: [5, 5]

init: Where the parameters start. sample, the default, draws one realisation of the prior’s shape — note that it is narrower than the prior itself by a factor of √2, which #217 tracks; prior starts at the prior mean, i.e. at whatever mean is set to; data starts at the observed visibilities themselves, RFI and all, attributed to the sky; zeros starts at an identically zero sky; and truth starts at the simulated astronomical visibilities. zeros and prior are the same starting point at the default mean: 0, and part company once the prior mean is the data. truth reads the tab-sim .zarr beside the measurement set, so it is only possible on a dataset simulated with sim-vis — as is plots.truth, which reads the same truth; a config asking for either without a readable zarr stops before the run starts rather than deep inside component setup.
mean: The mean of the prior distribution: 0 (equivalently zeros), the default, or data, the observed visibilities.
freq_pad_factor: This defines the size of the padding used when modelling the signal in the Fourier domain. The signal is modelled in the Fourier domain where periodicity is assumed on some interval. If freq_pad_factor: 1.0 is given then the interval is the interval of the data itself and will lead to periodic solutions.
time_pad_factor: This defines the padding used in the time axis of the signal. It is the time axis equivalent to freq_pad_factor.

baseline_block_size: The number of baselines GPVisAst transforms per step of its scan over the baseline axis: auto, the default, sizes the block so that one step’s padded grid stays inside a fixed budget; a whole number sets it outright; null puts every baseline in a single step. Turning the latent modes back into visibilities means padding them up to the padded Fourier grid, transforming, and cropping the padding away again, so doing every baseline at once holds an (n_bl, n_freq_pad, n_time_pad) array several times over — most of it discarded by the crop. Each padded axis is n + 2 * floor(n * (pad_factor - 1) / 2), so at the default factor of 2.0 it is about twice the data axis. The scan replaces n_bl in that shape with the block. It is purely a memory strategy: baselines are independent, so the result does not depend on it, and unlike the RFI scans there is no checkpoint on the body — the transform is affine in the parameters, so its derivative is a linear map with no primal intermediates to store.

ast_vis:GPVisAstDFT is the same component with the padded grid never formed: the pad, the transform and the crop are one linear, separable map, so the same visibilities come out of one small matrix per axis, sized by the modes that survive cutoff rather than by the padded grid. Its only transient is the intermediate of two matrix products, so auto there is a single step over every baseline – the scan has nothing large left to bound, and its stack of per-step outputs is a cost of its own. On one GH200, 2016 baselines at 150 integrations, value and gradient together, each case in its own process (the peak is a process high-water mark):

channels	`GPVisAst`, best block	`GPVisAst`, unblocked	`GPVisAstDFT`
8	273.9 MB, 1.0 ms	441.2 MB, 1.0 ms	218.6 MB, 0.2 ms
32	828.6 MB, 4.0 ms	1711.2 MB, 3.3 ms	683.8 MB, 0.4 ms
128	2386.4 MB, 14.7 ms	5926.7 MB, 12.7 ms	1799.1 MB, 1.1 ms

The two agree to round-off in value and in gradient, which is what its tests check. An explicit baseline_block_size is still honoured on both.

auto exists because a block that does not bind costs scan steps for nothing. On a single-channel observation, where the padded grid is 1 by 180, a fixed block of 128 split 4560 baselines into 36 steps and cost 13 % of the optimiser’s time with no memory saved; sized from the grid, the same observation runs in one step, and a wide band still blocks.

cutoff: This is the relative cutoff for Fourier components. The power spectrum is calculated and a mode is kept when its value is above cutoff times the largest on each axis, reducing the number of parameters to fit. It does not change the width of the prior: the normalisation to std happens after the cut, so changing it changes how many modes are fitted, and with them the correlation structure the prior can express, but not how much sky flux it expects. That is why it sits out here with freq_pad_factor and baseline_block_size — the other knobs on the size of the problem — rather than inside gp_cov, which is about its statistics. It must be a positive number below 1, since it is relative to the largest mode on each axis and at 1 every mode is cut. null is 1e-6.
gp_cov: The prior covariance of the astronomical visibilities: a width, the correlation scales that set the knee, and the roll-off past it. It is diagonal in the Fourier domain, which is how it is applied, and its diagonal is a power spectral density — but what you are setting here is the covariance of a Gaussian process, not the power spectrum an observation is trying to measure.

rfi.gp_cov takes the same keys at the same level, and the same validator checks both: each of std, corr_freq, corr_time and fov_deg must be a finite positive number, and gammas an ordered pair of them — one for the frequency axis and one for the time axis, in that order. A key the section does not read is refused by name rather than ignored, so a gamma written for gammas is caught rather than silently doing nothing.

The two sections differ in exactly two ways, and both are stated rather than left to be inferred from where a key sits. This one takes its time scale as a field of view, because a sky source’s coherence time is its fringe rate — a consequence of how far off axis it is, rather than a number to pick — where an emitter’s is a property of the emitter and is given directly as corr_time. And the two read corr_freq: null differently: no roll-off at all here, half the observed band for the RFI. Both are defaults for a signal that is not there to be characterised, and they decline to commit in different ways: the astronomical null names no delay scale at all, while the RFI null names half the observed band.

The parameters of the covariance are defined as

std: The width of the prior on the astronomical visibilities, in Jy — the standard deviation of the complex visibility about its prior mean, \(\sqrt{E|V-\mu|^2}\). At the default mean: 0 that is exactly rms|V|. Under mean: data the prior is centred on the observed visibilities and std is the scatter allowed around them rather than the total amplitude, which is then \(\sqrt{\mathtt{std}^2 + |\mu|^2}\). The power spectrum is normalised to it, so this is the prior width: not its square, not its square root, and not a constant times it.

data, the default, measures it for you: rms|V| per baseline over the samples nothing has flagged. At the default mean: 0 that is exactly what std is defined as, so there is no conversion in between (under mean: data the two differ, since the prior is then centred on the visibilities rather than on zero), and it gives one width per baseline rather than one for all of them.

It measures whatever is in the unflagged samples, RFI included. Where the contamination has been flagged, the estimate is the sky. Where it has not, the estimate is the contamination: on the shipped 8A simulation, whose RFI is unflagged because modelling it is the job, data returns about 11 Jy against a true sky under 3. It prints what it measured and warns when nothing – neither the MS nor the gain table – excludes a single sample. Note that it reads the MS’s flags directly, whatever data.flags is set to, along with the samples no gain table could calibrate — those answer different questions. data.flags decides whether tabascal declines to fit a flagged sample; this decides where it is willing to measure a clean sky from. A strong emitter flagged by some other task is exactly the data tabascal is here to recover, so you want it in the likelihood and out of the amplitude estimate, and that is what you get. An uncalibratable sample is out of both: it kept a unity gain where its neighbours were divided by a real one, so it is not on the same flux scale as the data around it. #220 tracks estimating it from the RFI model’s own view of which samples are clean, which is what would make the default estimate the sky rather than the sky plus whatever sits on top of it.

To set it by hand, take the visibility amplitude in a channel with no visible RFI. The shipped configs are set from the true sky of the simulation each runs on, at half the widest baseline’s amplitude, so every baseline’s sky sits inside 2 σ of the prior. Covering the widest at 1 σ instead was measurably worse on both counts: on the FixedOrbit case it cost 17 % of the astronomical accuracy (NRMSE 0.2089 against 0.1784) and raised the coherent bias from 1.2 σ to 1.8 σ. Where no such channel exists, the amplitude of the shortest baselines is the usual stand-in, since they see the most sky flux and the least resolved structure. It is a prior width, not a measurement, so being within a factor of a few is enough — what matters is that it is the right order of magnitude and that it means the same thing from one config to the next.

Why it is a width and not a power. cutoff exists to drop modes that carry no power, and gammas and fov_deg shape the roll-off; none of them is a statement about how bright the sky is, so none of them moves this number. The spectrum is normalised to std after the cut, which is what makes that true.

Every baseline gets this width. The time-axis knee is that baseline’s own maximum fringe rate, so each baseline’s spectrum is normalised separately; unnormalised, the shipped configuration spans a factor of 1.23 across its baselines.

It is a deliberately weak, generic prior rather than a model of the sky: a real sky puts more correlated flux on short baselines than on long ones, which resolve it out, so giving every baseline the width you read off the shortest ones is permissive on the longest. That permissiveness is what lets an astronomical model absorb RFI, so it is worth knowing about; #219 tracks giving the width a baseline-length envelope. The alternative — letting the width fall out of each baseline’s fringe-rate knee — is not better: that varies with temporal bandwidth, which is not sky amplitude either.

rfi.gp_cov.std means the same thing — the typical rms|V| in Jy — so one measurement sets either prior. What differs is internal: vis_ast is the modelled quantity, while rfi_A is a per-antenna amplitude the visibility is quadratic in, so the RFI prior takes a square root that this one does not. In both sections the number you read off the data is the number you write down.
corr_freq: The bandwidth, in Hz, over which the astronomical signal stays correlated. It sets the knee along the frequency axis, and that axis is a delay: the modes there are fftfreq(n_freq, chan_width), in inverse Hz, so the knee is 1 / (2 * pi * corr_freq) seconds — the same conversion rfi.gp_cov.corr_freq makes. null, the default, is no roll-off along that axis at all: every delay mode is kept and none is preferred.

It has no effect on a single-channel observation, where the only delay mode is zero, which is why every shipped config leaves it unset. On a wide band it is worth setting: the sky is smooth in frequency, so its power belongs at low delay, and the value to give is a bandwidth comparable to the one observed rather than a narrow one. A useful sanity check: the delay axis only runs out to 1 / (2 * chan_width) — 2.4 µs for 209 kHz channels, whatever the channel count — and the knee sits at k0 = 1 / (2 * pi * corr_freq). Comparing the two says how much of the axis the roll-off touches. At corr_freq = chan_width / pi (67 kHz for those channels) the knee lands exactly on the far end, where the outermost mode is still suppressed to 0.32 of the peak: that is where the roll-off bites hardest across the band, not where it stops mattering. It takes about ten times smaller — under 7 kHz — before the outermost mode is within 3 % of the peak and the axis is effectively flat. A correlation bandwidth of one observed band puts the knee at 1 / (pi * n_freq) of the axis, keeping the lowest delays and cutting the rest, which is what a sky smooth in frequency should look like.
fov_deg: The field of view in degrees used to set the maximum astronomical fringe rate (the knee k0 of the time-axis power spectrum). It is the full field of view, i.e. the angular diameter out to the first null of the primary beam; the maximum source offset from the phase centre is fov_deg / 2. When omitted, it defaults to the primary-beam field of view of the telescope, 2 * 1.22 * lambda / D (null-to-null), from the dish diameter D and frequency read from the MS file.
gammas: The rate of drop off in the power spectrum. As \(\gamma \rightarrow \infty\), the power spectrum tends to a Gaussian with width set by corr_freq in the frequency axis and inferred from fov_deg in the time axis.

A fixed sky of discrete sources

The GP sky above is flexible: it has per-baseline freedom, so a free per-antenna gain can be absorbed into it and is a flat direction of the likelihood. A source with a known position, flux and shape is rigid, and so anchors the gain. That sky is configured with ast.point_sources and needs two components, FixedDiscreteSky to read the catalogue and DiscreteSkyVis to turn it into visibilities:

model:
  components:
    - ast_signal:FixedDiscreteSky   # must come before DiscreteSkyVis
    - ast_vis:GPVisAst              # optional
    - ast_vis:DiscreteSkyVis
    - gains:UnitaryGains

ast:
  point_sources:
    - {name: Fornax A, ra: 50.6738, dec: -37.2083, I: 750.0, ref_freq_mhz: 154.0, alpha: -0.77}
  source_block_size: 128

FixedDiscreteSky writes the sky into the model state and DiscreteSkyVis reads it, so FixedDiscreteSky must be listed first. DiscreteSkyVis accumulates into vis_ast rather than assigning it, so it composes with ast_vis:GPVisAst — with both listed, vis_ast is the GP plus the fixed sources, and since vis_ast is zeroed before the components run the two may be given in either order. Neither component has any free parameters.

Sources are either the inline list above or a path to an OSKAR sky model file:

ast:
  point_sources: /path/to/sky.osm

The inline form requires ra and dec in degrees and I in Jy, and optionally takes name, ref_freq_mhz, alpha, Q, U, V, rm, fwhm_major_arcsec, fwhm_minor_arcsec and position_angle_deg. That list is exhaustive and case-sensitive: any other field is an error rather than being ignored, since a typo would otherwise change the source in silence — fwhm_maj for fwhm_major_arcsec would leave a Gaussian modelled as a point. An optional field may be given as null to mean “unset”, but any other unreadable value (alpha: "", ref_freq_mhz: false) is an error rather than a fallback to the default. The file form is the OSKAR sky model, which is what Karabo emits: one source per line, fields separated by whitespace and/or commas, # starts a comment, and blank lines are skipped.

I is the integrated flux of the source, and that is what appears on a zero-length baseline in any direction — there is no 1/n applied to it. (The measurement equation carries the sky brightness as B/n, because the solid-angle element is dl dm / n; for a discrete source of integrated flux S the brightness is S times a delta function of solid angle, and the two factors of n cancel exactly.)

#	Column	Units	Meaning
1	RA	deg	Right ascension
2	Dec	deg	Declination
3	I	Jy	Stokes I at the reference frequency
4	Q	Jy	Stokes Q
5	U	Jy	Stokes U
6	V	Jy	Stokes V
7	Reference frequency	Hz	Frequency at which the Stokes fluxes are quoted
8	Spectral index	—	`alpha` in `I(nu) = I (nu / nu_ref)**alpha`
9	Rotation measure	rad/m²
10	FWHM major	arcsec	Gaussian major axis; `0` for a point source
11	FWHM minor	arcsec	Gaussian minor axis; `0` for a point source
12	Position angle	deg	Major-axis position angle, north through east

A row need not carry all twelve fields, but only three lengths are meaningful, and they follow OSKAR’s own fixed-format reader:

3 to 9 columns — the leading columns above, with the rest defaulting to zero. ra dec I is a valid flat-spectrum point source.
11 columns — the legacy layout: columns 1-8 as above, then FWHM major, FWHM minor and position angle, with no rotation measure (it defaults to zero). This is not the 12-column layout with one field missing; read as though it were, the major axis lands in the rotation-measure column and a perfectly ordinary Gaussian is rejected as polarised.
12 columns — the full modern layout.

10 columns, or 13 and more, are rejected: 10 is a half-specified shape (a major axis with no minor axis or position angle) and 13+ is not the format.

Two points about the columns themselves:

Polarisation is parsed but not modelled. A non-zero Stokes Q, U or V, or a non-zero rotation measure, is rejected with an error naming the source, rather than being silently dropped. The whole format is accepted so that modelling polarisation later (issue #151) widens what these values mean rather than changing what the file may contain.
A spectral index needs a reference frequency. A non-zero spectral index with no positive reference frequency is an error, since falling back to a flat spectrum would put the source at the wrong flux in every channel with nothing in the output to say so. A zero spectral index is a flat spectrum and needs no reference frequency.

A source with a non-zero FWHM is an elliptical Gaussian, which multiplies the point-source visibility by the uv-plane envelope

\[G(u,v) = \exp\left(-\frac{\pi^2}{4\ln 2}\left(a^2 u'^2 + b^2 v'^2\right)\right)\]

for FWHM \(a\) (major) and \(b\) (minor) in radians, where \((u', v')\) is the baseline in wavelengths rotated into the source frame, \(u' = u\sin\phi + v\cos\phi\) along the major axis and \(v' = u\cos\phi - v\sin\phi\) along the minor axis. The position angle \(\phi\) follows the radio convention, measured from north (the \(m\) axis) through east (the \(l\) axis), so position_angle_deg: 0 puts the major axis north-south and a north-south baseline is the one that resolves the source out. A zero FWHM gives \(G = 1\) exactly, so points and Gaussians are the same code path and the same component.

A source more than 90 degrees from the phase centre is modelled as given, not rejected — the w term is computed exactly over the whole sphere — but it raises a warning naming the sources, since in practice it means a swapped or mis-signed coordinate rather than a real field.

source_block_size: The number of sources DiscreteSkyVis handles per step of its scan over the catalogue, a whole number defaulting to 128. The geometric delay array is (n_bl, n_time, n_src), which for a large catalogue is the biggest array in the model; the scan replaces n_src in that shape with source_block_size, at the cost of recomputing each block in the backward pass. It is purely a memory strategy — the result does not depend on it.
uvw_sign: The sign applied to each of the \(u\), \(v\) and \(w\) axes of the measurement set’s UVW column before DiscreteSkyVis uses it — three values, each exactly +1 or -1, defaulting to [-1, -1, -1].

The visibility equation above is written for the baseline \(b = \mathrm{ANTENNA2} - \mathrm{ANTENNA1}\), but a UVW column may hold either that baseline or its negative, depending on the software that wrote it, and nothing in the data says which. The default negates because that is what tab-sim writes — bl_uvw = ants_uvw[a1] - ants_uvw[a2], the same convention tabascal uses to form baselines internally — so on simulated data the fixed sky lands on top of the simulated sources. A measurement set from another toolchain may carry the opposite convention, and then the right value is [1, 1, 1]; the per-axis form is there for the rarer case of a column that differs on only some axes.

Getting it wrong is not obvious from the fit. Negating all three axes conjugates every visibility, which is the sky mirrored through the phase centre: the fixed sources sit in the wrong place, and since that corruption is smooth it is largely what a gain solved against a fixed sky will absorb, leaving the optimisation looking converged. The correct value is a property of the dataset rather than of the model, so it is worth establishing once per instrument or pipeline — for instance by checking the UVW column against the antenna positions and ANTENNA1/ANTENNA2 of the same row, or by imaging a bright known source and confirming it is not reflected through the phase centre.

Note that the catalogue fluxes are in the same scale as the data the model is fit to. With data calibrated to Jy these are physical Jy; without that, the data are in raw correlator units and a Jy catalogue flux is meaningless.

RFI signal

The rfi section defines the prior distribution over the RFI signal. An example of this section is given below.

The prior is a Gaussian process in the Fourier domain rather than over a dense covariance matrix. That is a numerical requirement and not a preference: a Cholesky factorisation’s jitter is absolute, so at the variance this prior carries it regularises far too weakly and returns NaN in single precision.

rfi:
  init: sample
  mean: 0
  min_elevation: 0
  freq_pad_factor: 2.0
  time_pad_factor: 2.0
  n_int_freq: 1
  time_int_factor: 1
  baseline_block_size: 128
  poly_interp_stencil: 1
  path_order: 3
  cutoff: 1e-9
  gp_cov:
    std: null
    corr_freq: 1e6
    corr_time: 24
    gammas: [3, 3]

All parameters in this section that overlap with those of the ast section have the same definition, except that init and mean accept one more value:

gp_cov.std: The width of the prior on the RFI signal, in Jy: the RFI’s typical rms|V|. The same quantity as ast.gp_cov.std — read an amplitude off the data and write it here, for either prior.

The two get there by different arithmetic, and this is where what the number means differs. vis_ast is the modelled quantity, so the astronomical prior’s width is the latent’s width. rfi_A is a per-antenna amplitude and the visibility is quadratic in it (\(V^\text{RFI}_{pq} = A_p A_q^*e^{i\Delta\phi}\)), so the per-antenna width is \(\sqrt{\texttt{std}}\) and carries units of \(\sqrt{\text{Jy}}\). Both translations are internal: the number you write is rms|V| in Jy in both sections.

The spectrum is normalised to std after the cut and after the roll-off, so — exactly as on the astronomical side — gammas and cutoff change which modes are fitted and how they correlate, not how much RFI the prior expects. (sum(pk) is std / 2; the factor is the complex latent, see _LATENT_POWER.)

It is one source’s width, and it is stated for rfi_signal:ComplexRFIVarAnt. The normalisation sum(pk) = std / 2 is exact and survives the transform: latent_to_signal pads the coefficients, inverts with norm="forward" and crops, so each retained coefficient reaches every output point with unit magnitude and \(E\lvert A\rvert^2 = 2\sum \texttt{pk} = \texttt{std}\). For independent antennas that gives \(E\lvert V_{pq}\rvert^2 = \texttt{std}^2\). Padding and cropping change which modes exist and how they correlate, not the amplitude.

It is an expectation, so any one observation scatters around it — by more where the correlation structure leaves fewer independent samples in the grid, which is a property of the spectrum rather than of the normalisation.

It is also the instantaneous visibility of a zero-mean, unmasked source: the width the prior is on, not a prediction of what a run will see. Three model steps sit in between, all of them deliberate:
- A non-zero mean — data, est, matched-filter — adds its own power, \(E\lvert V_{pq}\rvert^2 = (\texttt{std} + \lvert m_p\rvert^2)(\texttt{std} + \lvert m_q\rvert^2)\) for mean amplitudes \(m\). Sources then carry non-zero mean visibilities as well, so their total no longer generally scales as \(\sqrt{N}\).
- min_elevation zeroes a source while it is below the cut, so one visible for a fraction \(f\) of the observation shows \(\texttt{std}\sqrt{f}\) across the whole of it — and one that never rises, exactly zero.
- The visibility kernels average the fine grid, and fringes that turn within an integration cancel there.
Two further factors sit between it and the RFI a run realises even at zero mean, and neither is corrected for. They are properties of the model, and correcting them would make the same number mean different widths in different configurations:
- N satellites realise \(\sqrt{N}\) times it, for ComplexRFIVarAnt. The width applies to each source and their visibilities add — in quadrature, because each source has zero mean visibility. Three satellites at std: 10 expect a total RFI near 17 Jy. ComplexRFIConstAnt’s sources do not add in quadrature: each has a non-zero mean visibility \(E[V_s] = \texttt{std}\,e^{i\phi_s}\), so they add coherently by however much the geometric phases align.
- rfi_signal:ComplexRFIConstAnt realises another \(\sqrt{2}\). It broadcasts one amplitude to every antenna, so its visibility is \(\lvert A\rvert^2\) where ComplexRFIVarAnt’s is \(A_p A_q^*\) with independent draws, and \(E\lvert A\rvert^4 = 2(E\lvert A\rvert^2)^2\).
Both are pinned by tests that sample each component through its own antenna structure.

data measures it from the observed visibilities, the same measurement ast.gp_cov.std: data makes — the same function, in fact. Two things differ, and only one of them is a choice:
- It is a scalar, since this prior normalises a single spectrum where the astronomical model carries a width per baseline.
- The MS’s flags are kept. The astronomical estimate excludes them because whatever a flag means, the sample is not clean sky. The reverse does not follow: a flag says something is wrong here, not RFI is here, and a dead antenna is neither RFI nor a scale to set an RFI prior from. Measuring the flagged samples is also biased high even when they are RFI, wherever the flagger thresholds on amplitude: the flagged half is then the bright half — on the shipped 8A simulation that returns 1.54x the true RFI where measuring everything returns 1.015x. Samples no gain table could calibrate are dropped, which is the one exclusion both priors share.
It measures the total RFI and hands it to each source, so it does not correct for the \(\sqrt{N}\) above: on an N-satellite run the prior is that much wider than what was measured. Deliberate — a number and a measurement that produced different priors would be worse — and stated here because it is the one place the two std: data options differ in effect rather than only in mask.

It also tends to sit above the RFI itself, since the sky and the noise are in the visibilities too — tends to, not always: the sky and the RFI can cancel coherently on a given sample. A good one exactly where it matters — RFI that dominates the sky dominates the measurement — and far too wide where the RFI is faint, which is the mirror of the astronomical estimate’s own limitation. #220 is what fixes both.

On the shipped 8A simulation, whose three satellites give a true RFI rms|V| of 11.0 Jy: data measures 11.2 per source and the prior realises about 19. The null default sets 46.4 per source, which carries the same \(\sqrt{3}\) and realises about 80. Roughly 1.8x the true RFI against 7.3x — both figures on the same footing, which is the comparison that matters.

When null it is twice the largest visibility in the observation. That is a maximum where this key says typical — 46.4 Jy against a true RFI rms|V| of 11.0 on the shipped 8A simulation, where data gives 11.2. #227 tracks making it a statistic this key’s name describes.
init / mean: matched-filter (alias mf) estimates the per-satellite light curves directly from the visibilities the run has already loaded, by matched-filtering them against the known satellite trajectory phase, and seeds the RFI amplitude with them. It is the same seed as est without the file: no imaging step, no rfi.est, and no matching of light curves to satellites by name, since the estimator is handed satellites.norad_ids and returns the curves in that order. See Estimating the light curves from the data.

The only additional parameters are

gp_cov: The prior covariance of the RFI signal, taking the same keys as ast.gp_cov at the same level and checked by the same validator: std, corr_freq, corr_time and gammas. Every one of them is optional, and each null has its own meaning, given below.
- corr_freq: the bandwidth, in Hz, over which the RFI signal stays correlated, setting the knee along the frequency axis. That axis is a delay, so the knee is 1 / (2 * pi * corr_freq) seconds — the same conversion ast.gp_cov.corr_freq makes. null is half the observed band, which is not what the astronomical corr_freq: null means: that one is no roll-off at all. The difference is which default declines to commit. Note that a flat delay spectrum is not a smooth signal: equal power at every delay leaves the channels uncorrelated, its normalised channel covariance being [1, 0, 0, ...]. So the astronomical null states no frequency structure rather than gentle structure — which costs nothing on the single-channel observations every shipped config runs, and is exactly why ast.gp_cov.corr_freq tells you to set a value on a wide band. An emitter is coherent over some band by construction, so here half the observation is the least-committal guess that is still a coherence scale.
- corr_time: the time, in seconds, over which the RFI signal stays correlated, setting the knee along the time axis. null is half the observation. It is given as a time where the astronomical block asks for fov_deg, because an emitter’s coherence time is a property of the emitter, while a sky source’s is its fringe rate — a consequence of how far off axis it sits rather than a number to pick.
- gammas: the roll-off exponent on the frequency and time axes, in that order. null takes the component’s own value: [3, 3] for rfi_signal:ComplexRFIVarAnt, [100, 100] for rfi_signal:ComplexRFIConstAnt. The two have never agreed, and the difference is preserved rather than unified, since making them agree would change one of the two models rather than fix a bug.
cutoff: the relative power below which a k-mode is dropped from the latent grid. It therefore sets the number of fitted RFI parameters, which the run prints as (n_k_fq, n_k_tm) beside the resolved values — and nothing about what the prior believes, since the spectrum is normalised to gp_cov.std after the cut. That is why it sits beside freq_pad_factor and baseline_block_size rather than inside gp_cov, exactly as ast.cutoff does. It is relative to the largest mode on each axis, so it must be below 1 — at 1 every mode is cut and nothing is left to fit. Below 1 is necessary rather than sufficient: a value near enough to 1 to round to it in the working precision cuts everything too, which the run refuses with a message naming the cutoff rather than a shape error from inside the transform. null takes the component’s own value: 1e-9 for rfi_signal:ComplexRFIVarAnt, 1e-6 for rfi_signal:ComplexRFIConstAnt.
min_elevation: Elevation in degrees below which a satellite’s RFI signal is held at zero, so it is only modelled while it is up. The default is 0, which masks a satellite exactly while it is below the geometric horizon. Set it to null to disable masking entirely and model every satellite over the whole observation.

The default is the one the base configuration ships, so omitting the key means 0, not no mask. Running with no mask at all takes an explicit min_elevation: null.

While a satellite is below the horizon it contributes no signal, but an unmasked model still carries a full set of free parameters for it over those times. Those parameters have no signal of their own to constrain them, so they are free to absorb signal that belongs elsewhere — the astronomical sky, or another RFI source — to the extent that the RFI signal prior admits it and the fringe rates overlap. Masking removes the parameters rather than relying on the fit to leave them alone. This is why 0 rather than null is the default: a satellite below the horizon is not a modelling choice, it is simply not there. It shows in the image domain: an unmasked run can reconstruct a static feature near the horizon out of those parameters, which the mask removes.

Each satellite gets its own in-view window, evaluated on the observation time grid and expanded over each integration, so an integration is never partially masked. Setup fails if a satellite is never above the cut, since it would then be modelled nowhere. tabascal light-curve is the exception: it is measuring rather than fitting, so a satellite that never rose comes back as a zero curve, named in a warning, and the satellites that were up are still measured.

Raising the cut above 0 additionally excludes the low-elevation part of each pass, where the fringe rate is lowest and the overlap with other components is therefore greatest. How far to raise it is observation-dependent and is not currently calibrated, so no value above 0 is recommended here. Note that masking is about which parameters exist, not about subtraction quality, and reduced \(\chi^2\) is largely insensitive to it — judge the effect on the recovered sky model.

The RFI signal is modelled on a grid finer than the data, then averaged back down onto it. The fine grid is n_freq * n_int_freq by n_time * n_int_time, where each count is the number of fine samples per data cell on that axis. The two axes are configured differently, because only one of them can be estimated: there is no observable that fixes the frequency count, so n_int_freq is given directly, while the time count follows from how fast the RFI fringe winds and is therefore derived rather than written down. time_int_factor scales that derivation. There is no rfi.n_int_time key. The time count is derived, not chosen, so there is nothing to set: rfi.time_int_factor scales the derivation instead.

n_int_freq: This is the amount of over-sampling in the frequency domain that is used and then averaged back down to the data sampling rate. It therefore determines the number of samples per frequency channel that are used in the averaging to correctly calculate the fringe-winding loss (band-smearing). Band-smearing can be caused by both the phase variation over the channel width due to the geometric phase as well as the intrinsic signal of the RFI sources. The default is 1, i.e. no over-sampling.
time_int_factor: In the time axis the number of integration samples needed to accurately model fringe-winding loss (time-smearing) is calculated based solely on the fringe rate due to the movement of the RFI source as well as the signal to noise ratio with

\[N^T_\text{int} \geq \pi \nu_F \Delta t \sqrt{\frac{\lvert V^\text{RFI}_\text{inst} \rvert}{6 \sigma_n}}\]

where \(N^T_\text{int}\) is the number of integration samples used per time step, \(\Delta t\) is the integration time for a single sample, \(\nu_F\) is the fringe frequency of the source due to its movement, \(\lvert V^\text{RFI}_\text{inst} \rvert\) is the instantaneous RFI visibility amplitude, and \(\sigma_n\) is the visibility noise of a single data point. This parameter (time_int_factor) determines the factor by which to increase this oversampling.

For the legacy fine-grid route the estimate is per baseline, since \(\nu_F\) is, and TabConfig.estimate_rfi_sampling then reduces it to the single fine-grid count: the count is the largest per-baseline rate, rounded up to a size with enough divisors for the stride binning that groups baselines by how finely each needs to be sampled. That rounding applies to every run, not only to one using a RiemannVisVariable component: the number of divisors required is at least min_time_bins + 1, and only the extra divisors a RiemannVisVariable component needs are conditional on selecting one. min_time_bins and max_time_bins bound the number of quantile levels the grouping places, which caps the number of stride groups rather than guaranteeing that many — distinct levels can round onto the same stride. With no satellites configured there is no fringe rate to estimate from at all, every baseline falls back to a single required sample, and the count comes out at 2. So the effective count is not exactly time_int_factor times the formula above: the factor scales the per-baseline rates going in, and the binning decides what comes out. The value actually used is printed during setup.

baseline_block_size: The number of baselines RiemannVis calculates per step of its scan over the baseline axis, a whole number defaulting to 128. The Riemann sum is formed on a (n_bl, n_rfi, n_freq_fine, n_time_fine) fine grid before anything is reduced, which is n_rfi * n_int_freq * n_int_time times the size of the visibilities it reduces to and, under reverse-mode automatic differentiation, is what the tape holds; the scan replaces n_bl in that shape with baseline_block_size, or with n_bl itself where that is smaller, and recomputes each block in the backward pass rather than keeping it. It is purely a memory strategy — baselines are independent, so the result does not depend on it — at the cost of recomputation and of one scan step per block. null is the setting for a single block over every baseline: the checkpoint stays, so the backward pass still recomputes the fine grid instead of storing it, but the grid is formed whole. That bounds the tape and not the peak, and it is the setting that trades the memory back for the scan’s step overhead. Only RiemannVis reads it: the FFI kernels bound the same term inside the compiled kernel, and the Variable components carry their own baseline grouping. See RFI-visibility kernels.
poly_interp_stencil: Read by the PolyInterpVis family, the data-grid route (see The data-grid RFI route and its kernel boundary). The number of data cells on either side of a cell that its interpolant runs through, on each axis, so that 2 * poly_interp_stencil + 1 coarse samples make each cell’s fine samples: 1, the default, is the quadratic through a cell and its two neighbours; 0 holds each cell’s value across it; 2 is the quartic through five. An axis with fewer cells than the stencil needs – a single channel, say – takes the widest stencil it can hold, which for one cell is its own value. A whole number of at least 0.
poly_time_sampling: Read by PolyInterpVisVariable and PolyInterpVisVariableFFI. max_groups is 1 or 2 (default 2); split_at is a positive whole sampling count or null (default). The polynomial routes round each fringe-rate requirement up to a positive odd count: the even-count fine_offsets grid is shifted half a sample left, and its phase error is not controlled by the midpoint-quadrature estimate above. The non-variable polynomial route uses the maximum of these counts, with no legacy divisor rounding. Each variable group builds independent interpolation tables at its own count and compacts its antenna inputs, so one group reproduces that non-variable route. With two groups allowed, setup searches every threshold between distinct requirements and minimises sum(n_g * number_of_antennas_in_group). Shared antennas count in both groups. This score is materialisation work, not predicted runtime. split_at restricts the search to requirements at or below that count versus those above it; empty splits and splits that do not strictly beat one group fall back to one group. Equal scores favour fewer groups, then the lowest threshold. The legacy stride keys do not affect polynomial sampling; with no satellites its count is 1.
poly_analytic: Read by rfi_vis:PolyInterpVisHybrid and rfi_vis:PolyInterpVisHybridFFI. Baselines whose rounded requirements exceed quadrature_limit use an analytic time integral; the rest retain quadrature. quadrature_limit: null uses the measured GH200 VJP crossover: 166 samples in single precision, 56 in double, because the optimiser runs the transpose every iteration. Re-measure on very different hardware; quadrature’s memory traffic makes operation counts a poor predictor. An explicit non-negative count overrides the cut; 0 sends every baseline through the analytic route. Empty groups create no operator. segments splits the analytic integral into equal pieces, reducing local curvature by its square; terms sets the curvature series length. Defaults are segments: 2, terms: 6, cubic_terms: 3: measured accuracy was 1.21e-5 at a fifth of the 4/16 cost. This measurement is not a bound for arbitrary geometry. The broader stress tests through 3 curvature turns at the cell edge use an explicit 4/16 expansion; increase segments or terms for larger curvature or tighter accuracy. cubic_terms sets the residual cubic-phase expansion length; 0 drops it for comparison. Segments and curvature terms must be positive whole numbers, and cubic terms non-negative. The analytic route integrates the full amplitude polynomial against quadratic phase and a perturbative cubic correction; delay derivatives above order three are omitted. Frequency integration and its offsets are unchanged. The FFI route requires ri_kernels from interp-analytic and differentiates amplitudes through both groups; analytic phase and delay inputs are fixed. The older poly_time_sampling options are read only by the two quadrature-only variable components.
path_order: Read by trajectory:FixedOrbitCoarse alone. The degree of the polynomial in time fitted through each cell’s fine-grid path, whose coefficients – the path and its first path_order time derivatives at the cell centre – rfi_vis:PolyInterpVis rebuilds the fine phase from. 3, the default, keeps the phase difference a baseline sees to a few thousandths of a degree over an 8 s integration on the test geometry; 1 is a few tenths of a degree there and worse on longer baselines. A cell with fewer fine samples than the degree needs is fitted to the degree it can hold. A whole number of at least 0.

RFI light curve estimates

rfi.est points at a measured light curve file, used by init: est and mean: est to seed the RFI signal. This is the interchange format between tabascal and whatever measures the light curves, so it is deliberately strict.

The file is either a .zarr store (read with xarray.open_zarr) or a .npz, and must contain all four of

name	shape	contents
`light_curves`	`(n_src, n_time, n_freq)`, real	apparent flux per source, in Jy
`norad_ids`	`(n_src,)`	NORAD id of each row of `light_curves`
`times`	`(n_time,)`	UTC Modified Julian Date, in days, strictly increasing
`freqs`	`(n_freq,)`	frequency in Hz, strictly increasing

and should also carry

name	shape	contents
`time_scale`	scalar	`"utc"` — the scale `times` is on. A store attribute in the zarr form, an array in the npz

In the zarr form the four required names are coordinates or variables of light_curves, whose dimensions must be exactly norad_ids, times and freqs — declared in any order, since they are identified by name and transposed on read. time_scale is a store attribute. A minimal writer:

import numpy as np, xarray as xr

xr.Dataset(
    {"light_curves": (("norad_ids", "times", "freqs"), curves)},
    coords={"norad_ids": np.array([25544, 27386]), "times": times_mjd_utc, "freqs": freqs_hz},
    attrs={"time_scale": "utc"},
).to_zarr("light_curves.zarr")

time_scale is checked, not merely recorded: a file declaring anything other than utc is refused rather than converted, since the reader cannot know what another writer meant by it and converting there would make the reader a second place the format’s scale is decided. Fix the file instead — rewrite times as UTC MJD (tabascal.time.to_utc_mjd) and stamp utc.

A file that carries no time_scale is read as UTC with a warning. Light-curve files written by tabascal before the stamp existed took their times from the measurement set’s TIME column as declared, so one written from a UTC-declared MS — the overwhelmingly common case — is already correct, while one written from a TAI- or TT-declared MS is offset by the leap seconds and is indistinguishable from a correct file. Regenerate any estimate measured on a non-UTC MS with tabascal light-curve; the rest can be left alone or re-stamped.

Rows are matched to satellites by NORAD id, never by position, so the order of sources in the file does not have to match satellites.norad_ids. Samples are interpolated onto the observation’s own time and frequency grid, so the file’s sampling does not have to match the observation either.

Both are strict because their failure modes are silent. A light curve attached to the wrong satellite still has the right shape and still optimises — it just seeds the prior from another satellite. A file whose sampling is assumed rather than declared is resampled wrongly by an unknown amount. Neither surfaces as an error, only as a worse fit, so a file that cannot state which satellite and which sample times it describes is rejected rather than guessed at.

Times are absolute (MJD on a stated scale) rather than seconds from the start of a particular observation, so a light curve is interpretable on its own and can be reused across measurement sets covering the same pass. Both halves of that matter, and the stated scale is UTC. A Julian day number is a number until a scale says what it counts: a Measurement Set declares the scale of its TIME column in a MEASINFO record and is free to declare TAI, whose numbers name instants 37 s from the ones the same numbers name on UTC. An axis written on whatever the measuring MS happened to declare would be reusable only against another MS that happened to declare the same thing — and would be resampled by the difference without anything raising. tabascal reads the observation’s own times onto UTC before sampling an estimate, and tabascal light-curve writes them the same way, so both ends of the format are on one scale.

light_curves is a flux in Jy, not the modelled amplitude rfi_A. The RFI visibility is quadratic in rfi_A (\(V^\text{RFI}_{pq} = A_p A_q^* e^{i\Delta\phi}\)), so rfi_A carries units of \(\sqrt{\text{Jy}}\) and the estimate is seeded with \(\sqrt{\lvert \text{light\_curves} \rvert}\). Supplying an amplitude where a flux is expected is squared away silently, so the value is wrong rather than the shape — give the flux the source would show in the visibilities, on the same scale as rfi.gp_cov.std.

Some further details:

light_curves must be real — the magnitude \(\lvert \hat{S} \rvert\). A complex array is rejected rather than truncated: the cast to float64 would keep \(\text{Re}(\hat{S})\) and drop \(\text{Im}(\hat{S})\) behind nothing but a numpy warning, which on an uncalibrated column discards most of the signal without changing the shape of the result. tabascal light-curve writes the magnitude under this name and keeps the native complex estimate alongside it as light_curves_complex, which the reader ignores.
Each satellite must appear exactly once. A repeated NORAD id has no single answer to which row belongs to it, and resolving that by file order is the thing id-matching exists to avoid, so it is rejected — merge the passes or drop one before using the file as an estimate.
Labels that are not integer NORAD ids never match a satellite and are dropped, so a file may carry named sources (e.g. Fornax A) alongside the satellites without filtering beforehand, and those may repeat freely.
Samples outside the file’s coverage are zero, on either axis — the file says nothing there, which is the same “no signal known” convention the elevation mask uses. An axis of length 1 is held constant instead, since a single sample carries no gradient to interpolate along; a single-frequency light curve therefore applies across the whole band rather than being zeroed outside it.
The file does not have to cover every satellite in the fit. Satellites with no light curve are initialised at zero and named in a warning, so light curves can be measured for a subset — the bright or well-characterised sources — while the rest are still modelled and fitted, just without an informative starting point. It is an error only if no configured satellite is found, which would otherwise silently reduce the whole estimate to zeros.
A file written by tabascal light-curve --fit-offset carries the along-track offset fit beside the curves: tau_best (n_src,), tau_grid (n_tau,), z2_tau (n_src, n_tau, n_freq), z2_best, best_chan, significance, null_mean, null_std, detected (n_src,) bool, r_best (n_src, n_time, n_freq) and offset_threshold_sigma. Readers of the format ignore all of them. tau_best is not a diagnostic, though: it is the along-track offset each curve was extracted at, so a run seeded from such a file has to model the same trajectory — through the epoch-shifted records of the next bullet — or it fits a light curve measured on one trajectory against a model of another. See Fitting the along-track time offset.

Estimating the light curves from the data

rfi.init: matched-filter (alias mf) needs no file at all: it measures the light curves from the visibilities the run has already loaded. For a satellite on a known trajectory the RFI contribution to baseline \((p, q)\) is \(A_p A_q^* e^{i(\phi_p - \phi_q)}\) with \(\phi\) the geometric phase, so the unit-modulus template \(T_{pq} = e^{i(\phi_p - \phi_q)}\) de-rotates it. The maximum-likelihood estimate of the source visibility at each channel and timestep is the inverse-variance-weighted, de-rotated baseline average

\[\hat{S}[f, t] = \frac{\sum_{pq} w_{pq} T_{pq}^{*} V_{pq}}{\sum_{pq} w_{pq}}, \qquad w_{pq} = \frac{1}{\sigma_{pq}^2},\]

with standard error \(1/\sqrt{\sum w}\). The satellite’s fringe adds coherently after de-rotation while the sky and the noise do not, so \(\hat{S}\) isolates the satellite, and \(\sqrt{\lvert \hat{S} \rvert}\) is the per-antenna amplitude the model is seeded with — the same quantity est reads out of a file.

The weights are the run’s own noise, resolved per baseline and per channel as far as the MS resolves it (see data.noise); the template carries no gain. Uniform weights on an array whose antennas differ in sensitivity would over-weight the loud baselines. With no usable noise anywhere the curves are still measured, but they are then both unweighted and unscaled: the error and the z statistic come back as NaN rather than as a floor derived from a noise nobody stated. Flagged samples and autocorrelations are excluded, a satellite is not filtered for while it is below rfi.min_elevation (those times seed at zero, exactly as the elevation mask holds them there), and a channel and timestep where every baseline is flagged seeds at zero rather than at a measured value.

The same estimate is available as a standalone tool, tabascal light-curve, which writes it in the interchange format above so it can seed a later run through rfi.est. See Usage.

The z statistic reported alongside the curves reads Re(S_hat) and so assumes the column is phase calibrated. The magnitude statistic reported beside it is robust to a phase common to every baseline, but not to an uncalibrated antenna gain, which decorrelates the coherent sum and shrinks the estimate itself — see Usage.

Two limitations are worth knowing. The estimator assumes the satellite is exactly where its orbit record says it is: a position error scatters the per-baseline phases and costs coherence, which shows as an under-estimated flux rather than as an error. And it does not filter the astronomical signal out first, so a bright source in the field contributes to \(\hat{S}\) wherever its fringe rate overlaps the satellite’s.

Satellites

The satellites section determines which satelites to include in the model and the prior distribution to use for there trajectories. An example is given below.

satellites:
  norad_ids: [20452, 38833, 45854]
  norad_ids_path: null
  extra_orbit_dir: null
  extra_orbit_max_age_days: null
  remote_max_age_days: 3
  cache_reuse_max_age_days: 1

norad_ids: List of the NORAD IDs of the satellites to include. TABASCAL requests the record whose epoch is closest to the observation from the IAU CPS SatChecker service (via the satchecker-client package) — its get-nearest-omm endpoint for observations from 2026-07-12 onwards and get-nearest-tle before that, falling back to the other archive if the first has nothing acceptable. Cache misses run concurrently with a bounded five-worker pool; no account or credentials are required. Every ID listed here must resolve to an acceptable record: otherwise preflight stops before reading the visibilities and names each failure. TABASCAL never silently drops a configured satellite from the RFI model.

An empty list (or null) is valid only for a model that does not use a satellite trajectory component — a stationary-RFI or astronomical-only run. If the model.components list includes one that consumes orbital records (FixedOrbit, Orbit, NoDragOrbit), configuring no IDs is a configuration error rather than a run that models nothing. Either separator form of a component reference is recognised, so trajectory:FixedOrbit and trajectory.FixedOrbit behave identically here.
norad_ids_path: Optional path to a text file of NORAD IDs, one per line; blank lines and # comments are ignored and malformed lines are reported with their line number. When set it takes precedence over norad_ids, and the -np/--norad-path CLI flag takes precedence over both.
extra_orbit_dir: Optional path to an additional directory of local orbit files, searched per NORAD ID before the managed cache and SatChecker. Every *.json file in the directory is considered; files must be pandas-oriented JSON tables carrying either NORAD_CAT_ID, TLE_LINE1 and TLE_LINE2, or NORAD_CAT_ID, EPOCH and the seven OMM element columns. The kind is inferred, so a Space-Track gp/gp_history export drops in unconverted. For each requested satellite the valid record whose epoch is closest to the observation is chosen — by epoch distance, regardless of format — and, if it is accepted (see extra_orbit_max_age_days), it wins outright and no service call is made for that satellite. Files that cannot be read or lack either required column set are skipped. Records that fail validation are rejected, allowing that satellite to fall through to the managed cache and SatChecker. Legacy date-named files and the bundled Space-Track fixtures remain supported. This directory can be given at runtime with the --extra-orbit-dir flag. tabascal can also write one itself: tabascal light-curve --fit-offset --write-shifted-tle DIR saves the detected satellites’ records with their epochs moved by the fitted along-track offset, so a run pointed at DIR models the corrected trajectory rather than the one the elements arrived with, and tabascal search writes the same records for every satellite it detects — see Records with a fitted time offset. The search also sets this key for you: the satellites fragment it emits names the directory it wrote, with extra_orbit_max_age_days: null, so merging the fragment is enough. (The ORBIT_CACHE_DIR environment variable is a different thing: it relocates where the managed cache is stored, and is not an additional source of records.)
extra_orbit_max_age_days: Maximum allowed absolute difference, in days, between an extra_orbit_dir record’s epoch and the observation epoch. null (default) applies no age limit, preserving exact replay of used_orbits_*.json; 0 accepts an epoch match within TLE precision. A rejected local record falls through to the managed cache and SatChecker. The age comes from the record itself — line 1 for a TLE, the EPOCH field for an OMM — not from the filename or modification time.
remote_max_age_days: Hard ceiling, in days, on how far a SatChecker or managed-cache record may be from the observation. A TLE’s epoch is re-derived locally from line 1; an OMM has no lines, so its EPOCH field is used after being range-checked. Every accepted remote record’s source, provider, endpoint, signed offset and absolute age is logged. null explicitly removes the ceiling.

This ceiling is also what makes the endpoint fallback work. Neither SatChecker endpoint reports that it has nothing near the epoch requested — get-nearest-omm answers a pre-2026 request with its earliest record — so an over-age response is the signal that the record belongs to the other archive, and the other endpoint is then asked.

The default of 3 is provisional. It is a hard backstop against obviously unsuitable remote records — for one observation, SatChecker’s per-satellite fallback silently returned records ~31 days old, worth ~9,663 km of ISS position error — and not a claim that a three-day-old element set gives adequate positional accuracy. The calibrated, observation-specific suitability policy that should replace it is tracked in issue #101; it may end up rejecting records younger than three days for some orbits and baselines, or accepting older ones where independently justified.
cache_reuse_max_age_days: Request-avoidance threshold for the per-NORAD cache (default 1). A cached record this close to the observation avoids a request. An older cached record triggers an exact-epoch nearest lookup — including against the fallback archive, since holding a stale record is not the same as the archive having answered — but remains an offline fallback if it is within remote_max_age_days. A response replaces it only when strictly closer to the observation. null always reuses the nearest acceptable cached record. When both limits are set, this value must not exceed the hard ceiling.

Gains

The gains section defines the prior over the antenna gains and how they are initialised. An example is given below.

gains:
  init: prior
  amp_mean: 1.0
  amp_std: 10
  phase_mean: 0.0
  phase_std: 30
  ref_ant: null
  fix_flux_scale: true

init: How the gain parameters are initialised. prior (the default) starts at the prior mean. gains:ConstGains additionally accepts a path to a previously measured gain — see A constant gain per antenna.
amp_mean: The centre of the prior over the gain amplitude. It must be positive and finite. gains:ConstGains fits the log amplitude, so it is the median of a lognormal — the centre in log space — and it is read only when fix_flux_scale is false. See below.
amp_std: The standard deviation of the prior over the gain amplitude, as a percentage of amp_mean. amp_std: 10 with amp_mean: 1.0 is a 10 % spread. null defaults to 20, a 20 % spread.
phase_mean: The mean of the prior over the gain phase, in radians.
phase_std: The standard deviation of the prior over the gain phase, in degrees. null defaults to 180 — half a turn, which is effectively uniform over the circle. See below.
ref_ant, fix_flux_scale: Read by gains:ConstGains only; see below.
r_seed: The random seed the gain component draws with.

null means “unset”; 0 means zero. Every key in this section is defaulted when, and only when, it is null or absent — a written-down value is taken at its word. A zero seed is the seed it says. A zero width is an error naming the key, because it is a degenerate distribution rather than an absent one: it pins every gain to its mean and leaves the fit nothing to move. Negative and non-finite values are errors for the same reason.

The default prior widths

amp_std: null resolves to 20 % and phase_std: null to 180°. Both are wide on purpose, because in these components the prior width is not only a prior. The fitted parameter is always a standard normal \(z\), and the width is what carries it to the gain — directly in ConstGains (\(\texttt{phase} = \texttt{mean} + \texttt{phase\_std}\cdot z\), and \(|g| = e^{\texttt{amp\_std}\, z}\)). A narrow prior is therefore a short lever. Under a per-coordinate optimiser the step in \(z\) is set by opt.epsilon whatever the gradient is, so the phase moves by \(\texttt{phase\_std} \cdot \texttt{epsilon}\) per iteration. At the old 1° default and the default epsilon of 1e-2 that is 0.01° an iteration: the whole 500-iteration budget could not cross a radian, and a run whose antennas genuinely differed by tens of degrees ended where it started, looking converged. phase_std divides the phase on the way in as well (ConstGains reads a measured gain through it), so a narrow width also starts such a fit tens of \(\sigma\) from the prior mean. (UnitaryGains fits no gain and reads neither width.)

Why 180° is “effectively uniform”. The prior is Gaussian but the phase is an angle, so what the model sees is the wrapped normal,

\[p(\theta) = \frac{1}{2\pi}\left(1 + 2\sum_{k \ge 1} e^{-k^2\sigma^2/2}\cos k\theta\right),\]

which is uniform to within \(2e^{-\sigma^2/2}\): 1.4 % at 180°, and 5e-9 at 360°. The extra six decades of flatness buy nothing — no fit is sensitive to a 1.4 % tilt in the prior around the circle — and a full turn costs something, since the likelihood is \(2\pi\)-periodic in the phase and therefore periodic in \(z\) with period \(2\pi/\sigma\): doubling \(\sigma\) halves that period and puts twice as many whole-turn copies of every optimum inside the prior’s bulk. Half a turn is also the width at which \(|z| \le 1\) covers the whole circle, so every phase — including one read from a calibration table — starts inside the prior rather than outside it.

Set them explicitly when you know better. These are the widths for an array you have no prior information about. A well-behaved instrument with a recent calibration justifies a much tighter prior, and tightening it is what makes the prior do work; leaving it at the default only says that the data should decide.

A constant gain per antenna

ConstGains fits one complex gain per antenna, constant over time and frequency — the static direction-independent gain the array is known to have:

\[V^\text{OBS}_{pq} = g_p g_q^* \left( V^\text{AST}_{pq} + V^\text{RFI}_{pq} \right)\]

It adds only \(2 n_\text{ant} - 2\) parameters, and unlike a fixed gain it is constrained by the data.

model:
  components:
    - trajectory:FixedOrbit
    - rfi_signal:ComplexRFIConstAnt   # not ComplexRFIVarAnt — see below
    - rfi_vis:RiemannVis
    - ast_signal:FixedDiscreteSky
    - ast_vis:DiscreteSkyVis
    - gains:ConstGains

gains:
  ref_ant: null
  fix_flux_scale: true

A gain is only identifiable against a model term it cannot deform, which is what the component list above is for (issue #124):

Pair it with rfi_signal:ComplexRFIConstAnt. With the per-antenna RFI model ComplexRFIVarAnt the RFI amplitude \(A_p\) is already free per antenna, so \(g_p A_p (g_q A_q)^*\) is unchanged by \(g_p \rightarrow c_p g_p\) together with \(A_p \rightarrow A_p / c_p\): the gain is an exact flat direction of the RFI term and only the astronomical model constrains it. Setup warns when the two are combined rather than refusing — the pairing is a modelling rule, not a hard error — and names the flat direction.
The astronomical GP absorbs a gain the same way. ast_vis:GPVisAst has per-baseline freedom, so a gain solved against it alone is a reparametrisation of an already-free vis_ast. A rigid sky — ast_signal:FixedDiscreteSky with ast_vis:DiscreteSkyVis — is what anchors the gain.

The gauge. The gain is purely relative: it carries no absolute flux scale and no absolute phase, and both are removed by construction rather than fitted.

The overall phase is unobservable, so ref_ant’s phase is pinned to exactly 0 and the other \(n_\text{ant} - 1\) phases are free.
The overall amplitude is degenerate with the RFI source amplitude and the astronomical amplitude, so the log amplitudes are carried by \(n_\text{ant} - 1\) parameters on an orthonormal basis of the zero-sum subspace, giving \(\sum_p \log |g_p| = 0\), i.e. a geometric mean \(|g|\) of exactly 1. Left free it simply drifts — in one earlier run it settled at a median \(|g|\) of 0.70, with the sky model absorbing the reciprocal — which is a nuisance direction that buys nothing and slows convergence.

Both directions are removed from the parameters, not just from the value they map to. Writing \(n_\text{ant}\) amplitude parameters and subtracting their mean would give the same gains and the same prior, but would leave the all-ones direction of that latent space invisible to every visibility: flat in the likelihood however much data there is, curved only by the prior. Such a coordinate wrecks the conditioning of the optimisation and makes a likelihood-only Fisher matrix singular, so there is no such coordinate.

The prior on \(|g_p|\) is lognormal: amp_std is used as the standard deviation of \(\log|g|\), which agrees with a fractional spread to first order and keeps the gain positive by construction. amp_mean is the median of that prior — its centre in log space, and the value the fit starts at — rather than its arithmetic mean, which is the slightly larger \(\texttt{amp\_mean} \cdot e^{\sigma^2/2}\) for \(\sigma\) the log-space spread. And it is that only when the flux scale is free: under the zero-sum gauge the geometric mean of \(|g|\) is 1 by construction and amp_mean merely sets the scale that amp_std’s percentage is taken of, so a non-unit amp_mean with fix_flux_scale: true raises a warning saying so. A value that is not positive and finite is an error rather than a default.

ref_ant: The antenna whose phase is pinned to 0. null (the default) selects the first antenna with any unflagged data. An antenna every one of whose baselines is flagged everywhere is not constrained by any visibility, so it cannot be the reference the others are measured against, and naming one explicitly is an error rather than a silently unpinned fit.

One reference pins one connected group. A phase is measured relative to another antenna’s along a chain of baselines that carry data, so if the unflagged baselines split the array into groups that share no baseline, pinning ref_ant in one group leaves every other group’s overall phase unconstrained. Setup computes the connected components of the unflagged baseline graph and stops, naming the groups, rather than fitting a model with a flat direction in it. Fit the groups separately, or flag the smaller ones out of the run. Heavy flagging is fine as long as what survives still connects the array.
fix_flux_scale: Whether to keep the zero-sum log-amplitude constraint above. true is the default. false frees the overall amplitude and is accepted only with ast_signal:FixedDiscreteSky in model.components, since a fixed-flux sky is the one thing in the model that can set the scale; without it the run stops at setup with the degeneracy spelled out. Nothing else changes: the phase reference is still pinned either way.

The absolute flux scale is assumed to be set by the data — e.g. by a REAL_DATA_FLUXCAL column, or by the fixed sky above — and if it ever needs fitting it belongs in a separate scalar component rather than in the gain.

Initialising at a measured gain. gains.init optionally starts the fit at a gain measured elsewhere, which is a much better starting point than the prior mean:

gains:
  init: /path/to/gains.npz     # or /path/to/caltable.B

An .npz carrying the per-antenna gain under the key gain, shape (n_ant,) complex.
A calibration table, read with read_caltable() — one tabascal wrote with tab2MS, or one from CASA. A caltable is resolved over frequency and time and ConstGains is not, so it is reduced to the median \(|g|\) and the mean phase direction over each antenna’s valid samples, with a warning naming the largest deviation when the solutions actually do vary. Flagged solutions are dropped rather than counted as zero, and an antenna the table has no solution for at all falls back to unit gain, with a warning.

Either way the gain is projected into the gauge rather than taken as given: the mean log amplitude is subtracted (unless fix_flux_scale: false) and the phases are referenced to ref_ant. Projecting an already-projected gain changes nothing, so the same file can be handed back to a second run. A zero or non-finite gain is an error — the fit is in log amplitude, which such a value has no value at.