Skip to content

virgil.imaging

Regularizers, priors and helpers for image reconstruction with Image. Ways to choose a regularization weight are the L-curve's corner, the discrepancy principle, classic MaxEnt (LCurve.classic_maxent) and, for Gaussian-field images, the Laplace evidence (log_evidence).

Methods

CLEAN

Högbom's CLEAN (Högbom 1974, A&AS 15, 417) finds the peak of the residual dirty image and adds a fraction (the loop gain) of a point source there. For data linear in the image, that residual is proportional to \(-\partial\chi^2/\partial c_p\), the gradient of χ² with respect to the flux of a point at each pixel. clean uses that gradient directly, so it works for closure, kernel and DISCO phases and V² as well. With \(J e_p\) the change in the whitened residuals per unit flux at pixel \(p\), each iteration picks the pixel whose Gauss–Newton step lowers χ² most, the largest \(g_p^2/|J e_p|^2\) with \(g_p<0\), and adds gain times that step, \(-g_p/(2|J e_p|^2)\). This is matching pursuit; for linear data, where \(|J e_p|\) is the same everywhere, it is exactly Högbom's CLEAN. The norms \(|J e_p|\) are computed once, at the start, with one Jacobian–vector product per pixel. Normalizing by them matters next to an analytic star, where flux in a pixel is nearly the same as the star's: the gradient there is small, but so is \(|J e_p|\), and an unnormalized search would pile flux beside the star.

Components are added to a fixed base scene, with fluxes relative to its (for a System base, to its total, as siblings of its components). They are grey unless a spectrum is given, as the image is in SPARCO. Without a base, the components alone make the image, starting from one at the centre of the grid (closure phases do not fix the position, so this is also the anchor).

Each iteration only adds flux, so an early step that overshoots cannot be undone by later ones. Every refit_every iterations a major cycle (Clark 1980, A&A 89, 377; Cotton–Schwab, Schwab 1984, AJ 89, 1076) refits the fluxes of all the components at once by non-negative least squares on the linearized residuals: flux can move between components, and components whose flux falls to zero are removed. A last major cycle runs when CLEAN stops at the target or at max_iterations. The base is fixed unless base_priors names some of its parameters; each major cycle then fits those together with the fluxes with fit, since fitting them in turn would let the components absorb an error in the base, or the base absorb flux it does not model. A companion much closer to the star than λ/B has its flux and separation nearly degenerate, whatever fits it.

Multi-scale CLEAN (Cornwell 2008, IEEE JSTSP 2, 793, arXiv:0806.2228; Rau & Cornwell 2011, arXiv:1106.2745): with several scales_mas, each component is a pixel convolved with a circular Gaussian of FWHM \(s\) (\(s=0\) is a point), and the image is \(\sum_s G_s * F_s\). The search runs over every (scale, pixel) with the same score \(g^2/|J e|^2\), and the major cycles refit all scales together. Each shape is cut to the support and renormalized there, so a component's flux is the flux it puts in the image. Cornwell multiplies the peak residual at each scale by a bias \(1 - 0.6\,s/s_{\max}\) to favour small scales; the score here is the χ² decrease of each component's step, which does not depend on its size, so no bias is needed. scale_bias applies one, \(1-b\,s/s_{\max}\), if the broad components grab flux that belongs to points.

Iteration stops when χ² per data point reaches target_chi2_red (the discrepancy principle; it relies on correct error bars), after max_iterations, or when χ² has stalled: no pixel lowers it, or it has fallen by less than a fraction stall_tolerance over the last stall_window iterations, and a major cycle does not help. That happens when noise puts the truth's own χ² per point above the target; chi2_red shows how close it came.

Evidence

With a Gaussian-field image, whitened latents \(z\sim N(0,I)\) and Jacobian \(J\) of the whitened residuals with respect to them, the Laplace evidence at the MAP is

\[\log Z \approx -\tfrac12\chi^2 - L - \tfrac12|z|^2 - \tfrac12\log\det(I+J^\top J),\]

and the Laplace posterior of the latents is \(N(z_0,(I+J^\top J)^{-1})\). The log determinant is \(\sum\log(1+s^2)\) over the singular values \(s\) of \(J\), computed from the singular values of \(J\) (QR-reduced when far from square): the Gram matrix's rounding error can exceed the identity on high signal-to-noise data and leave \(I+J^\top J\) numerically indefinite. With the thin SVD \(J=USV^\top\) a posterior draw is \(z_0 + V\,\mathrm{diag}((1+s^2)^{-1/2})\,\varepsilon_1 + (I-VV^\top)\varepsilon_2\).

\(-\tfrac12\chi^2 - L\) is model_loglike at the MAP less its data-only normalization, \(-\sum\log\sigma - \tfrac n2\log 2\pi\) (in its von Mises form for unprojected phases, and with the correlated closure phases' determinant), which is the same for every model and hyperparameter on a given dataset and cancels in differences of \(\log Z\). For plain data \(L=0\). Otherwise \(L\) is the rest of the likelihood's normalizer: for data with calibration gains or closure offsets (which the likelihood marginalizes) it is \(\tfrac12\log\det\) of their covariance factor, and for extra observables the log ratio of effective to quoted errors. It is evaluated at the MAP as the likelihood evaluates it, so that evidences for data with different gain or offset widths compare; its curvature with respect to \(z\) is left out of \(J\), as Gauss–Newton leaves out the residuals' own second derivatives.

Error scale

In MacKay's evidence framework (MacKay 1992, eqs. 4.9–4.10; Bishop 2006, eqs. 3.91–3.95) the noise level is a hyperparameter. Write the noise precision as \(\beta=1/s^2\), so that the likelihood is \(\exp(-\beta\chi^2/2)\) with χ² on the quoted errors. The Laplace evidence is maximized when

\[\frac1\beta = s^2 = \frac{\chi^2}{N-\gamma}, \qquad \gamma=\sum_i\frac{\beta\lambda_i}{1+\beta\lambda_i},\]

where \(N\) is the number of independent data and the \(\lambda_i\) are the eigenvalues of the Gauss–Newton curvature \(J^\top J\) in the field's whitened latents, in which the prior's curvature is the identity. \(\gamma\) is the effective number of parameters the data measure: a direction with \(\beta\lambda\gg1\) is fixed by the data and counts as one, one with \(\beta\lambda\ll1\) is fixed by the prior and counts as none. It is solved for \(\beta\) by Newton's method, which converges monotonically from \(\beta=0\) because \(\gamma\) is concave in \(\beta\). Each measured parameter uses up one datum's worth of scatter, so χ² ≈ N − γ for honest error bars, and the ordinary \(s^2=\chi^2/N\) is biased low for the same reason as the 1/N estimate of a sample variance. Only the field's latents are counted in \(\gamma\); each other measured parameter raises \(s\) by about 1/(2N), negligible in SPARCO fits.

One scale per observable. A single \(s\) assumes every error bar is off by the same factor. The V² and closure phases of one instrument are calibrated differently: in a MATISSE N-band contest file the V² gave χ² per point 0.005 and the closure phases 0.49. With by_observable=True each kind of observable \(b\) ("vis", "phi" and each kind in extras, pooled over the datasets) gets its own precision \(\beta_b\) on its rows of the residual vector, and the evidence is maximized where, for every block,

\[\frac1{\beta_b}=s_b^2=\frac{\chi_b^2}{N_b-\gamma_b}, \qquad \gamma_b=\beta_b\operatorname{tr}(A^{-1}G_b), \qquad A=I+\sum_b\beta_bG_b,\]

with \(G_b=J_b^\top J_b\) from block \(b\)'s rows. \(\gamma_b\) is the sum over the block's rows of the diagonal of the hat matrix \(B^{1/2}JA^{-1}J^\top B^{1/2}\) (\(B\) the diagonal matrix of each row's \(\beta_b\)): the share of the measured parameters that block pays for. The blocks are coupled, so the equations are solved by iterating \(\beta_b\leftarrow N_b/(\chi_b^2+\gamma_b/\beta_b)\) (the same fixed point as MacKay's, but monotone, so it cannot oscillate) from \(\beta_b=N_b/\chi_b^2\) to a relative change of \(10^{-10}\). With one block it reproduces the single scale. \(N_b\) counts independent data: closure phases from four or more telescopes count their independent combinations, and their periodic penalty residuals belong to the block's χ² and \(\gamma_b\) but are not counted as data. The blocks' covariances must be \(s_b^2D\): calibration gains, closure-phase offsets, marginalized flux scales or differential phases with a finite prior_width add nuisance covariance that does not scale with the quoted errors, and raise a ValueError. Prefer the single scale when a block has few data: \(s_b\) scatters by about \(1/\sqrt{2(N_b-\gamma_b)}\).

References: D. J. C. MacKay (1992), "Bayesian interpolation", Neural Computation 4, 415 (doi); C. M. Bishop (2006), Pattern Recognition and Machine Learning, §3.5; S. F. Gull (1989), "Developments in maximum entropy data analysis", in Maximum Entropy and Bayesian Methods, Kluwer, 53–71, for the same \(N-\gamma\) argument applied to maximum entropy (LCurve.classic_maxent).

Regularizers and helpers for image reconstruction.

An image fit is a call to fit whose model contains an Image, with a prior on its log_brightness (see :func:image_priors) and usually a regularizer, which adds a penalty on the pixel fluxes b to the loss:

Regularizer Penalty Least squares (LM)? A prior density?
TSV w Σ (Δx b)² + (Δy b)² yes no
TV w Σ √((Δx b)² + (Δy b)² + ε²) no no
MaxEntropy w Σ b log(b / q) no no
Laplacian w Σ (∇²b)² yes no
StarletL1 w Σ √(s² + ε²) over starlet details s no no
LogSum w Σ log(1 + b / (ε b̄)) no no
Centroid ½ |centroid / σ|² yes yes

Differences Δ are between neighbouring pixels, with zeros beyond the edges. TSV and the Laplacian favour smooth images, TV piecewise-flat ones, maximum entropy images close to a default q; StarletL1 (an L1 norm of the starlet coefficients; an L1 norm of the pixels is constant, since they sum to one) and LogSum favour sparse ones, and LogSum is not convex, so start it from a clean model. l_curve sweeps the weight w.

When fit's model function returns one model per dataset, every regularizer acts on the first model only: right for transformed copies of one scene (e.g. a Rotated epoch), but an Image that appears only in a later model is not regularized.

Closure, kernel and DISCO phases do not fix an image's position: something must (an analytic star at the origin, a Centroid prior, or a centred prior mean). diagnose checks a fit for this and other common pitfalls.

TSV

Bases: _ImageRegulariser

Total squared variation: weight * Σ (Δx b)² + (Δy b)².

A quadratic smoothness penalty, so it can be fitted by Levenberg–Marquardt (Kuramochi et al. 2018, ApJ 858, 56).

Parameters:

Name Type Description Default
weight float

Strength of the penalty.

 required
path str

Path of the Image in the model (e.g. "env"); None if the model is the Image itself.

 None
smoothed bool

Penalize the Image's smoothed() pixels (see blur_fwhm_mas) instead of the raw ones. Default False.

 False

weight = np.asarray(weight, dtype=float) instance-attribute

path = path instance-attribute

smoothed = bool(smoothed) instance-attribute

__init__(weight, path=None, smoothed=False)

residuals(model)

value(model)

TV

Bases: _ImageRegulariser

Total variation: weight * Σ √((Δx b)² + (Δy b)² + ε²).

Favours piecewise-flat images with sharp edges (Rudin, Osher & Fatemi 1992, Physica D 60, 259). ε makes it differentiable where the image is flat. It is averaged over the four flips of the image, so it does not depend on which neighbour each forward difference pairs a pixel with. It is not isotropic: a sharp edge along a diagonal costs about a fifth more per unit length than one along a pixel axis.

Parameters:

Name Type Description Default
weight float

Strength of the penalty.

 required
epsilon float

Smoothing scale, as a fraction of the mean pixel flux (default 1e-2).

 0.01
path str

Path of the Image in the model.

 None
smoothed bool

Penalize the Image's smoothed() pixels instead of the raw ones (default False).

 False

weight = np.asarray(weight, dtype=float) instance-attribute

epsilon = np.asarray(epsilon, dtype=float) instance-attribute

path = path instance-attribute

smoothed = bool(smoothed) instance-attribute

__init__(weight, epsilon=0.01, path=None, smoothed=False)

value(model)

MaxEntropy

Bases: _ImageRegulariser

Maximum entropy: weight * Σ b log(b / q), with b and q unit-sum.

The relative entropy of the image with respect to a default image q (Skilling & Bryan 1984, MNRAS 211, 111; BSMEM, Baron & Young 2008).

Parameters:

Name Type Description Default
weight float

Strength of the penalty.

 required
prior array - like

The default image q, with the image's shape and positive where the image's support is; normalized here. Default: flat over the support.

 None
path str

Path of the Image in the model.

 None
smoothed bool

Penalize the Image's smoothed() pixels instead of the raw ones (default False).

 False

prefers_lbfgs = True class-attribute instance-attribute

weight = np.asarray(weight, dtype=float) instance-attribute

prior = None if prior is None else np.asarray(prior) instance-attribute

path = path instance-attribute

smoothed = bool(smoothed) instance-attribute

__init__(weight, prior=None, path=None, smoothed=False)

value(model)

residuals(model)

Residuals r with value = 0.5 Σ r² (b and q sum to 1): r = t √(2 w q φ(t)) for t = (b - q) / q and φ(t) = ((1 + t) ln(1 + t) - t) / t², so that the entropy can be fitted by Levenberg–Marquardt. Zero where b is.

metric_rows(model)

Rows whose JᵀJ is the curvature w b of the entropy in ln b (the Skilling–Bryan metric): a positive Gauss–Newton form. In the Image's softmax logits its diagonal is w b (1 - b).

Laplacian

Bases: _ImageRegulariser

Squared Laplacian: weight * Σ (∇²b)².

∇²b is the five-point Laplacian, with zeros beyond the edges. A quadratic penalty on curvature, smoother than TSV, so it can be fitted by Levenberg–Marquardt (Thiébaut & Young 2017).

Parameters:

Name Type Description Default
weight float

Strength of the penalty.

 required
path str

Path of the Image in the model.

 None
smoothed bool

Penalize the Image's smoothed() pixels instead of the raw ones (default False).

 False

weight = np.asarray(weight, dtype=float) instance-attribute

path = path instance-attribute

smoothed = bool(smoothed) instance-attribute

__init__(weight, path=None, smoothed=False)

residuals(model)

value(model)

StarletL1

Bases: _ImageRegulariser

L1 norm of the starlet details: weight * Σ √(s² + ε²).

s are the detail coefficients of starlet at every scale. An L1 norm favours few non-zero coefficients, so the image is built from few compact structures of any size (Honma et al. 2014, PASJ 66, 95). The coarse plane, which carries the flux, is not penalized.

Parameters:

Name Type Description Default
weight float

Strength of the penalty.

 required
scales int

Number of detail planes (default 4); the largest holds structure about 2**(scales - 1) pixels across.

 4
epsilon float

Smoothing scale, as a fraction of the mean pixel flux (default 1e-2).

 0.01
path str

Path of the Image in the model.

 None
smoothed bool

Penalize the Image's smoothed() pixels instead of the raw ones (default False).

 False

weight = np.asarray(weight, dtype=float) instance-attribute

scales = int(scales) class-attribute instance-attribute

epsilon = np.asarray(epsilon, dtype=float) instance-attribute

path = path instance-attribute

smoothed = bool(smoothed) instance-attribute

__init__(weight, scales=4, epsilon=0.01, path=None, smoothed=False)

value(model)

LogSum

Bases: _ImageRegulariser

Log-sum sparsity: weight * Σ log(1 + b / (ε b̄)).

b̄ is the mean pixel flux over the support. As ε → 0 the penalty counts the bright pixels: a smooth surrogate for the L0 "norm" of SQUEEZE (Candès, Wakin & Boyd 2008; Baron, Monnier & Kloppenburg 2010). It is not convex, so start the fit from a good image, and sweep its weight with l_curve(..., warm_start=False): pixels a strong weight has driven dark cannot recover in the weaker fits that follow.

Parameters:

Name Type Description Default
weight float

Strength of the penalty.

 required
epsilon float

The flux, as a fraction of the mean pixel flux, at which a pixel starts to count (default 1e-2).

 0.01
path str

Path of the Image in the model.

 None
smoothed bool

Penalize the Image's smoothed() pixels instead of the raw ones (default False).

 False

weight = np.asarray(weight, dtype=float) instance-attribute

epsilon = np.asarray(epsilon, dtype=float) instance-attribute

path = path instance-attribute

smoothed = bool(smoothed) instance-attribute

__init__(weight, epsilon=0.01, path=None, smoothed=False)

value(model)

CleanResult dataclass

The result of :func:clean.

Provisional

May change in any 0.4.x release without a deprecation cycle.

Attributes:

Name Type Description
model SourceModel

System(base=base, clean=Image(...)) (a System base's components and clean side by side). The Image is non-zero only on the components, and its flux is their total relative to the base. Without a base, the Image alone; with no components, the base alone.
components (array, shape(npix, npix))

The components' total flux in each pixel, Σ_s G_s * F_s over the scales, relative to the base's weight (without a base, unit sum).
pixel_scale_mas float

Pixel size in milliarcseconds.
chi2_red array

χ² per data point before each iteration; the last entry is the final model's.
stop str

"target" (target_chi2_red reached), "stalled" (χ² stopped falling, even after a major cycle) or "max_iterations".
components_by_scale (array, shape(n_scales, npix, npix))

The fluxes F_s at each scale, at the centre of their shape, in the units of components; each sums to that scale's share of the flux.
scales_mas tuple of float

The FWHM of each scale's Gaussian shape (0 for a point).

model instance-attribute

components instance-attribute

pixel_scale_mas instance-attribute

chi2_red instance-attribute

stop instance-attribute

components_by_scale = None class-attribute instance-attribute

scales_mas = (0.0,) class-attribute instance-attribute

__init__(model, components, pixel_scale_mas, chi2_red, stop, components_by_scale=None, scales_mas=(0.0,))

restored(beam)

The components convolved with beam: the "restored" image.

In the orientation of the pixel grid, in the units of components. Unlike radio astronomy's restored image, it does not include the residuals.

Centroid

Bases: _ImageRegulariser

A Gaussian prior centring an Image's flux on the origin.

The penalty is ½ |c / sigma_mas|², where c is the flux-weighted centroid of the Image on the sky, including its dra/ddec. This is a genuine prior density, so it may be used when sampling.

Parameters:

Name Type Description Default
sigma_mas float

Standard deviation of the centroid, in milliarcseconds.

 required
path str

Path of the Image in the model.

 None

probabilistic = True class-attribute instance-attribute

weight = 1.0 instance-attribute

sigma_mas = np.asarray(sigma_mas, dtype=float) instance-attribute

path = path instance-attribute

__init__(sigma_mas, path=None)

centroid(model)

The Image's flux-weighted centroid (dra, ddec) in mas.

residuals(model)

value(model)

Beam dataclass

A Gaussian approximation to the core of the dirty beam.

Attributes:

Name Type Description
major_mas, minor_mas float

Full widths at half maximum along the major and minor axes, in mas.
pa_deg float

Position angle of the major axis, North to East, in degrees.

major_mas instance-attribute

minor_mas instance-attribute

pa_deg instance-attribute

__init__(major_mas, minor_mas, pa_deg)

LCurve dataclass

The result of :func:l_curve.

Attributes:

Name Type Description
weights array

The regularization weights, in the order fitted (largest first).
chi2 array

Total χ² of each fit.
chi2_red (array, shape(n_weights, n_datasets))

χ² per data point of each dataset.
penalty array

The unweighted regularizer, value / weight, of each fit.
results list of FitResult

The fits.

weights instance-attribute

chi2 instance-attribute

chi2_red instance-attribute

penalty instance-attribute

results instance-attribute

__init__(weights, chi2, chi2_red, penalty, results)

corner()

The weight at the L-curve's corner, where it bends most sharply.

The curve is (log χ², log penalty) parametrized by log w; the corner is the interior point of largest curvature (Hansen & O'Leary 1993). Check it by eye: a sparse or noisy sweep has a noisy curvature. A fit with a non-finite χ² or penalty is left out.

discrepancy(target=1.0)

The weight at which χ² per data point reaches target.

Morozov's discrepancy principle: regularize as strongly as the data allow. With several datasets the binding one (the largest χ² per point) must reach the target. It interpolates linearly in log w and returns None if the sweep never crosses the target. The truth has χ² per point 1 ± √(2/N), so 1 is the natural target; it relies on correct error bars and over-regularizes if they are underestimated.

classic_maxent(data, path='env')

The maximum-entropy weight of Gull and Skilling's "classic MaxEnt".

For a sweep of MaxEntropy fits, the weight w at which -2 w S equals the number of well-measured directions in the image, N = Σ λ / (λ + w) (Gull 1989, Skilling 1989, both in Maximum Entropy and Bayesian Methods, Kluwer). S is the entropy (minus the sweep's penalty) and the λ are the eigenvalues of the data's Gauss–Newton curvature in the entropy metric diag(√b) JᵀJ diag(√b). It is the stationary point of the Laplace evidence in w, so it assumes correct error bars. Both sides are evaluated at each fit in the sweep and the crossing is interpolated linearly in log w; None if the sweep does not bracket it. Like log_evidence it uses the quoted errors, and rejects a sweep fitted with noise= terms or one model per dataset.

Parameters:

Name Type Description Default
data OIData or sequence of OIData

The data the sweep was fitted to.

 required
path str

Path of the regularized Image in the model (default "env").

 'env'

Diagnosis dataclass

The result of :func:diagnose.

Attributes:

Name Type Description
checks dict

Check name to value. Checks made per dataset or per Image are lists, in the order of the data or of the Images in the model.
warnings list of str

One sentence for each problem found, saying what to do about it. Empty if none.

checks instance-attribute

warnings instance-attribute

__init__(checks, warnings)

__str__()

starlet(image, scales=4)

The starlet (isotropic undecimated wavelet) transform of an image.

The à-trous algorithm with a B3-spline kernel (Starck, Murtagh & Fadili 2010): the image is smoothed scales times with the kernel's taps twice as far apart each time, and each detail plane is the difference between successive smoothings. Plane j holds structure about 2**j pixels across. Beyond the edges the image is zero.

Parameters:

Name Type Description Default
image (array - like, shape(ny, nx))

The image.

 required
scales int

Number of detail planes (default 4).

 4

Returns:

Name Type Description
details (Array, shape(scales, ny, nx))

The detail planes, finest first.
coarse (Array, shape(ny, nx))

What is left after the last smoothing; details.sum(0) + coarse is the image.

clean(data, npix, pixel_scale_mas, base=None, *, gain=0.1, max_iterations=1000, target_chi2_red=1.0, refit_every=50, stall_window=50, stall_tolerance=0.001, refresh_norms=False, base_priors=None, spectrum=None, support=None, init=None, rotation_deg=0.0, dtype='float64', scales_mas=(0.0,), scale_bias=0.0)

Build an image from components, one at a time: gradient CLEAN.

Provisional

May change in any 0.4.x release without a deprecation cycle.

Each iteration adds gain times the Gauss–Newton step at the pixel that lowers χ² most (Högbom's CLEAN for linear data, else matching pursuit); every refit_every iterations a major cycle refits all fluxes by non-negative least squares. See CLEAN. Source: Högbom (1974), references.

Parameters:

Name Type Description Default
data OIData or sequence of OIData

The data, fitted jointly.

 required
npix int and float

Grid size and pixel size (mas).

 required
pixel_scale_mas int and float

Grid size and pixel size (mas).

 required
base SourceModel and float

Fixed scene (e.g. a star) to add components to, not offset if a System; fraction of each step taken, in (0, 1].

 None
gain SourceModel and float

Fixed scene (e.g. a star) to add components to, not offset if a System; fraction of each step taken, in (0, 1].

 None
max_iterations (int, float and int)

Iteration limit; χ² per point to stop at; iterations per major cycle (0: none).

 1000
target_chi2_red (int, float and int)

Iteration limit; χ² per point to stop at; iterations per major cycle (0: none).

 1000
refit_every (int, float and int)

Iteration limit; χ² per point to stop at; iterations per major cycle (0: none).

 1000
stall_window int and float

Stop when χ² falls by less than this fraction in this many iterations.

 50
stall_tolerance int and float

Stop when χ² falls by less than this fraction in this many iterations.

 50
base_priors dict and bool

Base priors refitted, and pixel norms |J e_p| recomputed (costly), each major cycle.

 None
refresh_norms dict and bool

Base priors refitted, and pixel norms |J e_p| recomputed (costly), each major cycle.

 None
spectrum Spectrum and array - like

Shared spectrum (default grey; needs a base); boolean (npix, npix) mask, e.g. circular_support; non-negative starting fluxes (npix, npix) or (n_scales, npix, npix).

 None
support Spectrum and array - like

Shared spectrum (default grey; needs a base); boolean (npix, npix) mask, e.g. circular_support; non-negative starting fluxes (npix, npix) or (n_scales, npix, npix).

 None
init Spectrum and array - like

Shared spectrum (default grey; needs a base); boolean (npix, npix) mask, e.g. circular_support; non-negative starting fluxes (npix, npix) or (n_scales, npix, npix).

 None
rotation_deg float and str

Grid position angle (match data.uv_grid.rotation_deg, see Image); precision.

 0.0
dtype float and str

Grid position angle (match data.uv_grid.rotation_deg, see Image); precision.

 0.0
scales_mas sequence of float and float

FWHMs (mas) of the components, 0 for a point ((0.0, 3.0): both); in [0, 1), the down-weighting 1 - scale_bias s / max(scales_mas).

 (0.0,)
scale_bias sequence of float and float

FWHMs (mas) of the components, 0 for a point ((0.0, 3.0): both); in [0, 1), the down-weighting 1 - scale_bias s / max(scales_mas).

 (0.0,)

Returns:

Type Description
CleanResult

The model, components and χ² history.

starting_image(data, *, star=True, oversample=4.0, largest_mas=None, start='moments', hole_mas=None, field_fwhms=6.0, inner_elements=4.0, aperture_mas=None, pad_to_mas=None)

A starting model for an image fit, sized from the data.

Fits an analytic star (if star) plus a circular Gaussian envelope to the data, then builds an image of field_fwhms envelope FWHMs (capped by largest_mas and field_of_view), with pixels oversample times finer than Nyquist. Fit it with fit(start, image_priors(start), data, regularizers); see imaging with regularizers. With star=True the envelope's flux ratio is fixed unless you add a prior on "env.flux" to image_priors(start) (a lone Image is normalized, so has no flux).

Parameters:

Name Type Description Default
data OIData or sequence of OIData

The data.

 required
star bool

Whether an unresolved star at the origin is in the scene (it fixes the position; otherwise use a Centroid prior).

 True
oversample float

Pixels per Nyquist pixel.

 4.0
largest_mas float

Cap on the data-sized field, which never exceeds λ / B_min.

 None
start (moments, dirty)

Starting pixels: the fitted Gaussian, or the positive part of the dirty_image minus the star.

 "moments"
hole_mas float

Radius of a hole in the support under the star (0.5 * beam(data).minor_mas works well).

 None
field_fwhms float

Field in envelope FWHMs; least field in resolution elements λ / (2 B_max).

 6.0
inner_elements float

Field in envelope FWHMs; least field in resolution elements λ / (2 B_max).

 6.0
aperture_mas float

A telescope's field of view, as a ceiling on field and envelope width.

 None
pad_to_mas float

Least grid width (mas), padded with empty pixels; may exceed λ / B_min, but beyond it the data do not constrain the flux, so pair a large field with a support or priors.

 None

Returns:

Type Description
SourceModel

System(star=PointSource(), env=Image(...)), or the Image alone.

dirty_image(data, npix, pixel_scale_mas, flux_ratio=None)

The dirty image: direct synthesis of the visibilities, unregularized.

I(x) = Σ_k w_k Re[V_k exp(+2πi u_k · x)] / Σ_k w_k, with uniform weights on the informed samples: the image convolved with the dirty beam, plus noise; it peaks at one for a lone point source. For AMIGO DISCO data the visibilities are the least-squares estimate from the modes, with |V| = 1 on average on the shortest baselines (the modes are blind to the total flux).

Parameters:

Name Type Description Default
data OIData or sequence of OIData

The data; closure phases alone cannot make a dirty image (ValueError).

 required
npix int

Number of pixels on a side.

 required
pixel_scale_mas float

Pixel size in mas.

 required
flux_ratio float

Flux of extended emission relative to a star at the origin: the best-fitting point source at the origin is subtracted and the image scaled by (1 + flux_ratio) / flux_ratio, so the map shows the emission alone.

 None

Returns:

Type Description
(Array, shape(npix, npix))

In the orientation of render (East left, North up). It has negative sidelobes.

image_priors(scene)

Priors on the log-brightness of every Image in a scene.

Returns a priors dict for fit. An Image with a plain log-brightness array gets a flat prior, {"env.log_brightness": ImproperUniform(...)}, so that its pixels are constrained only by the data and the regularizers. An Image with a GaussianField gets standard-normal priors on the field's latents, {"env.log_brightness.latent": Normal(0, 1)}: the Gaussian-process prior, which needs no regularizer. A scene with several such Images (a compact and an extended component, say) gets one prior per Image. Add priors for any other free parameters (fluxes, offsets) to the dict.

nyquist_pixel_scale(data)

The largest pixel scale (mas) that samples the data's finest fringes.

That is λ / (2 B) for the longest baseline B in units of wavelength, over the unflagged samples of data (an OIData or a sequence of them; see OIData.unflagged). A flagged sample carries no signal, so its baseline does not count. Reconstructions normally use pixels 2–4 times smaller.

field_of_view(data, largest_mas=500.0, *, aperture_mas=None)

The largest field of view (mas) that the data measure.

The smallest of largest_mas (a cap you choose, 500 by default); the interferometric field, the largest λ / B over the unflagged samples with a non-zero baseline B (structure larger than that is not measured, and on a uv lattice a larger field would alias); the bandwidth-smearing field R λ / B_max, when the data carry channel widths (see smeared_model); and aperture_mas, if given, the field of a single telescope, about λ / D. A ceiling, not a recommended size: image_field sizes an image from the data and uses it as the limit. Flagged samples are left out (see OIData.unflagged).

image_field(data, fwhm_mas=None, *, field_fwhms=6.0, inner_elements=4.0, aperture_mas=None, largest_mas=500.0)

The width (mas) of an image field that suits the data.

The field is field_fwhms × fwhm_mas (without fwhm_mas, the second_moment_fwhm; with no visibility amplitudes to measure it, the ceiling), clipped to a floor of inner_elements × λ / (2 B_max) and a ceiling of field_of_view on the unflagged baselines (the ceiling wins if they cross).

Parameters:

Name Type Description Default
data OIData or sequence of OIData

The data.

 required
fwhm_mas float

The FWHM of the source in mas, to override the one from the visibilities.

 None
field_fwhms float

The field in FWHMs of the source (6 by default; 3 to 8 are reasonable). A larger field spends pixels on empty sky, where flux can scatter in islands.

 6.0
inner_elements float

The floor, in resolution elements λ / (2 B_max) (4 by default).

 4.0
aperture_mas float

A telescope's field of view, as a further ceiling.

 None
largest_mas float

A cap on the ceiling (500 by default, as for field_of_view).

 500.0

Returns:

Type Description
float

The field in mas.

second_moment_fwhm(data, v2_min=0.3)

The FWHM (mas) of the Gaussian that the short baselines see.

A circular Gaussian of FWHM θ has -ln V² = π² θ² ρ² / (2 ln 2) (ρ = B / λ, cycles per radian), and any normalized source has this form to second order in ρ. It is fitted by weighted least squares on the unflagged samples with V² ≥ v2_min, needing no image and no fit. It is the rms width of the whole source, a star included, so it is meaningless when a bright unresolved component (flux fraction above about 0.55) keeps V² above v2_min on every baseline: use envelope_fwhm there.

Parameters:

Name Type Description Default
data OIData or sequence of OIData

Data with visibility amplitudes (V², |V| or log |V|); projected (kernel, DISCO) observables and phase-only data have none.

 required
v2_min float

The smallest V² to use, 0.3 by default.

 0.3

Returns:

Type Description
float or None

The FWHM in mas: zero if V² does not fall with ρ (an unresolved source), and None if fewer than two unflagged samples have V² ≥ v2_min.

envelope_fwhm(data, star=True, aperture_mas=None)

The FWHM (mas) of the source's envelope, from a quick parametric fit.

This is the width that starting_image sizes its field from: a star (if star) plus a Gaussian envelope is fitted, and the FWHM is the envelope's or, if larger, the one the short baselines see. Unlike second_moment_fwhm alone it is not diluted by a bright star. Pass it to image_field.

beam(data)

The resolution of a dataset: the FWHM ellipse of its beam's core.

The dirty beam falls off near its peak as 1 - 2π² xᵀ M x, with M the weighted mean of u uᵀ over the samples; this matches a Gaussian of covariance M⁻¹ / 4π², whose FWHM ellipse is the usual "beam", roughly λ/B and elongated where the coverage is. Every informative sample has equal weight ("uniform weighting"): in AMIGO DISCO data, the uv points where the modes have weight (at least 1e-3 of the largest); otherwise the unflagged samples (see OIData.unflagged), since a flagged sample measures nothing and counting it would claim a resolution the data do not have.

Parameters:

Name Type Description Default
data OIData or sequence of OIData

The data.

 required

Returns:

Type Description
Beam

convolve_beam(image, pixel_scale_mas, beam)

An image convolved with a Gaussian beam: the image at the data's resolution.

The "restored" image of radio astronomy: it shows only what the data resolve, so convolve both reconstructions to compare them. The kernel is an elliptical Gaussian with the beam's FWHMs and position angle (North to East), normalized to unit sum on an odd grid centred on a pixel. Flux beyond the edge is taken to be zero, so structure within about a beam of the edge is dimmed.

Parameters:

Name Type Description Default
image (array - like, shape(ny, nx))

The image, in the orientation of render (East left, North up).

 required
pixel_scale_mas float

Pixel size in milliarcseconds.

 required
beam Beam

The beam (FWHMs in mas, position angle in degrees), usually beam(data).

 required

Returns:

Type Description
(Array, shape(ny, nx))

The convolved image.

l_curve(model, priors, data, regularizer, weights, others=(), *, warm_start=True, **fit_options)

Fit a model over a range of weights for one regularizer.

The weights are fitted from largest to smallest, each starting from the previous solution. Plot penalty against chi2 (log axes) to see the trade-off, and compare LCurve.corner and LCurve.discrepancy with the images either side: there is usually a wide range of good weights.

Parameters:

Name Type Description Default
model

As for fit.

 required
priors

As for fit.

 required
data

As for fit.

 required
regularizer (TSV, TV, MaxEntropy, Laplacian, StarletL1 or LogSum)

The regularizer whose weight is swept (its own is ignored).

 required
weights sequence of float

The weights to try (sorted internally, largest first).

 required
others sequence

Further regularizers kept fixed, e.g. a Centroid prior.

 ()
warm_start bool

Start each fit from the previous, stronger one (default), or every fit from model (and init). Use False for penalties that switch pixels off, such as LogSum: pixels a strong weight has driven dark cannot recover in the weaker fits that follow.

 True
**fit_options

Passed to fit.

 {}

Returns:

Type Description
LCurve

log_evidence(model, data, path='env', *, priors=None, regularizers=())

Laplace-approximated log evidence of a Gaussian-process fit.

At the MAP model, log Z ≈ -½ χ² - L - ½ |z|² - ½ log det(I + JᵀJ) up to a data-only constant, with z the whitened latents of the Image's GaussianField (or the Gaussian-process parts of a MultipoleDisk), J their Jacobian and L the likelihood normalizer; extra priors in priors enter through their log_prob at the MAP and their JᵀJ block of the Hessian. Compare it across hyperparameters: it is conditional on those held in model (or sample them: log-uniform priors from gp_priors for a MultipoleDisk, your own log-uniform priors on a field's sigma and length_mas for an Image). Exact only for a linear model, and it assumes correct error bars (see error_scale) and uses the quoted errors, so no noise= terms or per-dataset models. See Evidence. Source: MacKay (1992), references.

Parameters:

Name Type Description Default
model SourceModel or FitResult

The MAP model, or better the FitResult, so fitted noise= terms are caught (ValueError); a list of models raises TypeError.

 required
data OIData or sequence of OIData

The data it was fitted to.

 required
path str or tuple of str or None

Path of the Image (default "env"), a tuple of paths, or None for every such Image and disk (TypeError if there is none).

 'env'
priors dict

Priors of free parameters; non-latent parameters enter via log_prob and their JᵀJ block of the Hessian.

 None
regularizers sequence

Fit regularizers; probabilistic ones (Centroid) enter χ² and J.

 ()

Returns:

Type Description
float

laplace_samples(model, data, n, key, path='env', npix=None, fov_mas=None, *, priors=None, regularizers=())

Draws from the Laplace posterior of a Gaussian-field image.

The whitened latents z are drawn from N(z₀, (I + JᵀJ)⁻¹) about the MAP z₀, the approximation of log_evidence. With priors, the latents' covariance is the marginal one of the full Hessian diag(I, 0) + JᵀJ, with the other parameters' columns in J. Exact only for a linear model, assumes correct error bars, uses the quoted errors and is conditional on the hyperparameters held in model (or sample them, as for log_evidence).

Parameters:

Name Type Description Default
model SourceModel or FitResult

The MAP model, or the FitResult (fitted noise= terms raise ValueError, lists of models TypeError).

 required
data OIData or sequence of OIData

The data it was fitted to.

 required
n int

Number of draws.

 required
key Array

A JAX random key.

 required
path str or tuple of str or None

Path of the Image (default "env"), a tuple of paths, or None for every Image with a GaussianField (TypeError if there is none).

 'env'
npix int and float

Together, also render the whole model for each draw (see SourceModel.render); one alone raises ValueError.

 None
fov_mas int and float

Together, also render the whole model for each draw (see SourceModel.render); one alone raises ValueError.

 None
priors dict or sequence of str

Free parameters (priors or paths) besides the latents; they enter the Hessian through JᵀJ, with no unit-Gaussian prior.

 None
regularizers sequence

Fit regularizers; probabilistic ones (Centroid) enter J.

 ()

Returns:

Type Description
dict

"latents", an array (n, *latent.shape) (a dict by Image path if several are selected), and with npix and fov_mas "images", (n, npix, npix) unit-sum images of the whole model.

error_scale(model, data, path='env', *, by_observable=False, priors=None, regularizers=())

Re-estimate the scale of the error bars from a Gaussian-field fit.

Returns the factor s to multiply every error bar by for consistency with the fit (s < 1: bars too large), solving s² = χ² / (N − γ) with γ the effective number of parameters the data measure; see Error scale. Fit, call this, rescale with OIData.with_error_scale, refit, and do so before choosing hyperparameters (bars that are too large make the evidence over-regularize). s much greater than 1 means a failed or inadequate model, not a noise fix: unmodelled structure is absorbed into s. It uses the quoted errors, so no noise= terms or per-dataset models. Source: MacKay (1992), references.

Parameters:

Name Type Description Default
model SourceModel or FitResult

The MAP model, as for log_evidence.

 required
data OIData or sequence of OIData

The data it was fitted to.

 required
path str or tuple of str or None

Path of the Image (default "env"), a tuple of paths, or None for every such Image and disk (TypeError if there is none).

 'env'
by_observable bool

One scale per kind of observable instead of one for all. Data with gains, closure-phase offsets or other covariance that does not scale with the errors raise ValueError.

 False
priors dict or sequence of str

Free parameters (priors or paths) besides the latents; they enter through their JᵀJ block, with no unit-Gaussian prior.

 None
regularizers sequence

Fit regularizers; probabilistic ones (Centroid) enter J.

 ()

Returns:

Type Description
float or dict

The scale s, or with by_observable=True a dict from each observable kind ("vis", "phi", ...) to its scale, ready for OIData.with_error_scale. A warning is issued if the block scales do not converge.

chi2_by_observable(model, data, path=None)

Raw χ² of a model for each kind of observable, on the quoted errors.

The errors are the data's own (no noise= terms). Closure phases from four or more telescopes contribute their whitened combinations, and their periodic penalty rows (see whitened_residuals) add to the χ² but not to the count.

Parameters:

Name Type Description Default
model SourceModel

The model.

 required
data OIData or sequence of OIData

The data.

 required
path str or tuple of str

Path(s) of the Gaussian-process Images, as for log_evidence. If given, the result also holds n_gamma.

 None

Returns:

Type Description
dict

chi2, n and chi2_red (χ² per independent datum), each a dict from the kinds of observable ("vis", "phi", then any extras) and "total" to a number. With path, n_gamma: MacKay's number of well-determined latent parameters, Σ s²/(1 + s²) over the singular values s of the Jacobian of the whitened residuals by the latents.

short_baseline_residual(model, data, quantile=0.2)

Mean whitened visibility residual on the shortest baselines.

A reconstruction that loses flux on large scales leaves a coherent residual on the shortest frequencies though the total χ² looks fine (the test of MACIM, Ireland et al. 2006). The residuals are (model − data)/σ, so the mean of n of them has standard error 1/√n if the model is right.

Parameters:

Name Type Description Default
model SourceModel

The model.

 required
data OIData or sequence of OIData

The data. Visibilities in a linear combination (vis_mat, as in kernel-phase and DISCO data) and data with calibration gains are not supported: a ValueError.

 required
quantile float

The fraction of visibility samples, shortest first by B/λ, that are used (default 0.2).

 0.2

Returns:

Type Description
tuple of float

(mean, standard_error).

residual_image(model, data, npix, pixel_scale_mas, rotation_deg=0.0)

The residual dirty map of a model, for any observable.

The normalised χ² gradient -g_p / |J e_p| of adding a little flux at each pixel p (g the gradient of χ², J e_p the response of the whitened residuals), as clean uses. For linear data it is the dirty image of the residuals; for closure phases, the gradient the fit still feels. Positive pixels need more flux; pixels the model cannot respond to are zero.

Parameters:

Name Type Description Default
model SourceModel

The model, which needs a flux for the added point sources to be measured against.

 required
data OIData or sequence of OIData

The data.

 required
npix (int, float)

The centred grid of the map, as in clean.

 required
pixel_scale_mas (int, float)

The centred grid of the map, as in clean.

 required
rotation_deg float

The grid's rotation, as in clean.

 0.0

Returns:

Type Description
ndarray

Shape (npix, npix).

diagnose(model, data, regularizers=())

Check a model and its data for common imaging pitfalls.

Nothing is printed or warned: print the returned Diagnosis. Run it on the fitted model, with the fit's regularizers.

Parameters:

Name Type Description Default
model SourceModel

The model, normally containing at least one Image (the Image checks are skipped otherwise).

 required
data OIData or sequence of OIData

The data.

 required
regularizers sequence

The regularizers of the fit; a Centroid fixes the position.

 ()

Returns:

Type Description
Diagnosis

checks holds, per dataset, chi2_red (χ² per data point), phase_regime (fraction of model visibilities with |arg V| > 0.8π or |V| < 0.05, for projected phases only), chi2_red_by_observable (see chi2_by_observable) and short_baseline_mean (see short_baseline_residual; the shortest 20% of baselines, None for vis_mat operators such as kernel-phase and DISCO data, or calibration gains); per Image, pixel_scale_mas, edge_flux (fraction of the flux in the outer 2 pixels) and centroid_mas (East, North offset); and anchored (whether something fixes the position) and flip_dchi2 (Δχ² when every Image is rotated by 180° about the origin).

Sky-plane and uv-plane geometry shared by the source models.

Position angles follow the on-sky convention used throughout virgil: North is +y (up), East is +x (left, i.e. x decreases left-to-right across an image array's columns), and position angle is measured North to East, i.e. counter-clockwise from the top in a plot with East to the left. This matches the coordinate grid produced by :func:image_coordinates.

pixel_offsets(npix, pixel_scale_mas)

Sky offsets (mas) of the pixel centres along one image axis.

(centre - index) * pixel_scale_mas, centre at (npix - 1) / 2: dra along columns (East left), ddec along rows (North up).

set_ft_defaults(ft=None, cpu_min_points=None, cpu_points_per_side=None, cpu_min_side=None, gpu_min_points=None, gpu_min_side=None, eps_f32=None, eps_f64=None, toeplitz_min_points=None, toeplitz_points_per_pixel=None, toeplitz_max_bytes=None)

Set how ft="auto" chooses the Fourier transform of pixel images.

Only the arguments you pass change. On a CPU the NUFFT is used when there are at least max(cpu_min_points, cpu_points_per_side * N) points, N the image's larger side, and its smaller side is at least cpu_min_side; on a GPU when there are at least gpu_min_points points and the smaller side is at least gpu_min_side. eps_f32 and eps_f64 are its tolerances, relative to the zero-baseline visibility.

chi2="auto" radio terms use the dirty-beam form when they have at least max(toeplitz_min_points, toeplitz_points_per_pixel * n_ddid * N²) points, N the image side (ComplexVisData.term).

The settings are read when a function is traced: a jitted function that has already compiled keeps its old choice. The toeplitz_* ones are read at each fit.

Parameters:

Name Type Description Default
ft (auto, dft, nufft)

What models with ft="auto" use: "dft" always gives the exact sum, "nufft" always the NUFFT, "auto" chooses by size.

 "auto"
cpu_min_points float

The CPU thresholds above.

 None
cpu_points_per_side float

The CPU thresholds above.

 None
cpu_min_side float

The CPU thresholds above.

 None
gpu_min_points float

The GPU thresholds above.

 None
gpu_min_side float

The GPU thresholds above.

 None
eps_f32 float

NUFFT tolerances for float32 and float64 images.

 None
eps_f64 float

NUFFT tolerances for float32 and float64 images.

 None
toeplitz_min_points float

The size thresholds of chi2="auto" above; toeplitz_max_bytes caps the beams' device memory (about 24 N² n_ddid bytes).

 None
toeplitz_points_per_pixel float

The size thresholds of chi2="auto" above; toeplitz_max_bytes caps the beams' device memory (about 24 N² n_ddid bytes).

 None
toeplitz_max_bytes float

The size thresholds of chi2="auto" above; toeplitz_max_bytes caps the beams' device memory (about 24 N² n_ddid bytes).

 None

Returns:

Type Description
dict

The previous settings: set_ft_defaults(**old) restores them.