Skip to content

virgil.grid_fit

Grid searches. Contrast limits built on them are in virgil.limits.

Grid searches over model parameters.

Grids are built with indexing="ij": every output has one axis per grid key, in the order of grid. For {"dra", "ddec", ...} axis 0 is dra (East offset) and axis 1 is ddec (North offset), so 2D maps need a transpose to be shown as images with North up; plot_grid_map handles this.

The functions that optimize a flux at every grid position find it as the one key whose last part is flux (flux, comp.flux, ...), unless flux_param says otherwise. Fluxes are relative to the primary (at flux 1), so for a companion the flux is its companion/primary flux ratio; see virgil.limits for converting to contrast or Δmag. Contrast limits (Ruffio, Absil) are in virgil.limits.

LogUniform

Bases: NamedTuple

Log-uniform (scale-invariant) prior on the companion flux ratio.

p(f) = 1 / (f ln(f_max / f_min)) on f_min <= f <= f_max. The flux ratio is a scale parameter spanning decades, so this is the invariant measure of the scaling group (the Jeffreys prior under that group action), not the root-Fisher-information prior of the linearized likelihood, which has constant Fisher information in f and so would be flat. It is improper without bounds, so the evidence needs 0 < f_min < f_max (finite), in the units of the flux (companion/primary).

f_min instance-attribute

f_max instance-attribute

Gaussian

Bases: NamedTuple

Gaussian prior N(mean, sd**2) on the companion flux ratio.

Gaussian-prior evidence (a computational approximation; f may go negative): the prior has support on negative flux, so its evidence is only a convenient closed form. Prefer LogUniform.

mean instance-attribute

sd instance-attribute

likelihood_grid(model, data, grid, *, batch_size=None)

Evaluate the log likelihood at every point of a parameter grid.

Parameters:

Name Type Description Default
model SourceModel or class

Template model whose parameters at the paths in grid are varied (e.g. a System with paths such as "comp.dra"), or a model class called with grid's keys as keyword arguments (e.g. BinaryModelCartesian).

 required
data OIData

Data to fit.

 required
grid dict[str, array - like]

Grid axes, as a mapping from parameter name or path to 1D values (e.g. dra/ddec in milliarcseconds and flux as a companion/primary flux ratio). The output has one axis per key, in this order.

 required
batch_size int

Number of grid points evaluated at once. By default, enough for about 220 model visibilities on a CPU and 223 on other backends (GPU, TPU), and at least 256. Larger can be faster for small data; smaller bounds memory for large models.

 None

Returns:

Type Description
array - like

Log likelihood with shape tuple(len(v) for v in grid.values()). Axis k follows the k-th key (indexing="ij"), so for {"dra", "ddec", "flux"} axis 0 is dra: transpose a 2D slice before showing it as an image with North up.

optimized_likelihood_grid(model, data, grid, *, flux_param=None, batch_size=None)

optimized_flux_grid(model, data, grid, *, flux_param=None, batch_size=None)

linear_flux_grid(model, data, grid, *, flux_param=None, batch_size=None, n_iter=0, prior=None)

Linearized best-fit companion flux at every grid position, in closed form.

A fast first pass for companion searches, beside the iterative optimized_flux_grid, and the equivalent of fouriever's linear contrast map (lincmap). For a faint companion (flux f much smaller than 1) the whitened residuals r of virgil's likelihood are linear in f at fixed position: r(f) = r(0) + f g, where g = dr/df at f = 0 is computed exactly by forward-mode automatic differentiation of the model (no finite difference, no hand-derived closure-phase derivative). The weighted least-squares solution is then

f_hat = -(g . r(0)) / (g . g), sigma_f = (g . g) ** -0.5,

which is (gᵀ C⁻¹ (d - m₀)) / (gᵀ C⁻¹ g) with C the data covariance, because r is already whitened: correlated closure phases are handled exactly as in the likelihood (whitened_residuals), and every observable in the data (|V| or V², and phases) contributes. One model evaluation with its derivative per grid position, and no optimizer.

f_hat is not constrained to be positive (as with optimized_flux_grid and fouriever's lincmap), so noise gives negative values with SNR of either sign. Unlike fouriever, whose lincmap returns the variance 1 / (g . g) (and warns not to trust it), sigma_f here is the standard deviation.

Limitation, and n_iter. With n_iter=0 the linearization holds only for f much smaller than 1. The closure phase of a binary scales as f only to first order, with corrections of order f**2 (and f times the |V| change for amplitudes), so for a bright companion (for example f ~ 0.3) f_hat is biased, by tens of percent, and sigma_f is unreliable. n_iter Gauss–Newton steps soften this: each relinearizes at the current f_hat per pixel (g = dr/df at f_hat, then f_hat <- f_hat - (g . r(f_hat)) / (g . g)), and sigma_f comes from the final g, so a few steps (3 at f ~ 0.3) reach the optimizer's f_hat and the Laplace sigma_f, at n_iter + 1 model evaluations per pixel. The steps are not safeguarded, so for companions far brighter than the primary or strongly non-linear residuals they can fail to converge; use optimized_flux_grid to refine candidates. At Δ-phase residuals of order 1 rad the phase wrapping is not linear either. A companion at a position where g is nearly zero (for example at a null of the baselines) has a large sigma_f and so a small SNR.

Parameters:

Name Type Description Default
model SourceModel or class

Template model whose parameters at the paths in grid are varied (e.g. a System with paths such as "comp.dra"), or a model class called with grid's keys as keyword arguments (e.g. BinaryModelCartesian). At zero flux it must reduce to the primary alone (flux 1).

 required
data OIData

Data to fit.

 required
grid dict[str, array - like]

Grid axes, as for optimized_flux_grid: dra/ddec in milliarcseconds, plus a flux key, whose values are ignored (the flux is solved for, at f = 0), but which must be present to name the parameter. The output has one axis per coordinate key (every key except flux_param), in this order.

 required
flux_param str

The key of grid holding the flux, e.g. "comp.flux". By default, the one key whose last part is flux.

 None
batch_size int

Number of grid points evaluated at once, by default enough for a fixed number of model visibilities; see likelihood_grid.

 None
n_iter int

Number of Gauss–Newton refinement steps after the first linearization at f = 0 (default 0, the closed-form result).

 0
prior numpyro LogUniform or Normal

Prior on the flux ratio f. By default none, and the posterior and Bayes-factor fields of the result are None. With a prior they hold the posterior mean and sd and the marginal-likelihood detection map log_bayes_factor (see Returns); log B > 0 favours a companion at that pixel. All of these hold in the linear model about the final linearization point, i.e. exactly only where the residuals are linear in f over the posterior (f much smaller than 1, or after enough n_iter for the point to sit near the posterior); the position is not marginalized. Give numpyro's dist.LogUniform(low, high) or dist.Normal(mean, sd) with scalar parameters (the LogUniform and Gaussian named tuples of this module are still accepted); a bare (mean, sd) tuple is an error.

Recommended: dist.LogUniform(f_min, f_max), the scale-invariant (Jeffreys, under the scaling group) prior for a flux ratio, p(f) = 1 / (f ln(f_max / f_min)). The flux ratio is a scale parameter spanning decades, so the prior is the invariant measure of the scaling group, not the root-Fisher prior of the linearized likelihood (which would be flat). It is improper without bounds, and the evidence needs a proper prior, so both bounds are required (0 < f_min < f_max). The Bayes factor depends on them, as it must for a scale prior: for f_hat well inside the bounds, widening them by a factor changes log B by about -Δ ln(ln(f_max / f_min)) (the Occam factor). Choose them from the physics or the data, for example f_max the brightest companion you would entertain and f_min a little below the faintest contrast the data can reach (its dynamic range, e.g. the smallest flux_error over the grid). The likelihood in f is Gaussian with mean f_hat and sd sigma_f, so Z = ∫ N(f; f_hat, sigma_f**2) p(f) df / N(0; f_hat, sigma_f**2), computed by fixed 256-node Gauss–Legendre quadrature in ln f over the part of the bounds where the likelihood is not negligible; the posterior mean and sd come from the same quadrature. It is vmappable and jit-compatible.

dist.Normal(mean, sd) gives a Gaussian-prior evidence (a computational approximation; f may go negative). With P = g . g + 1 / sd**2 (g the final whitened derivative) the posterior is Gaussian with mean (g . (g f_hat) + mean / sd**2) / P and sd P ** -0.5, and the log Bayes factor against f = 0 is the closed-form Gaussian evidence ratio

log B = -0.5 log(sd**2 P) + (g.g f_hat + mean/sd**2)**2 / (2P) - mean**2 / (2 sd**2)

(Luger, Foreman-Mackey & Hogg 2017, arXiv:1710.11136).

 None

Returns:

Type Description
LinearFluxGrid

A named tuple, always of the same type, whose fields are arrays with one axis per coordinate key (axis 0 is the first, e.g. dra):

  • flux: best-fit flux ratio (companion/primary), unconstrained in sign.
  • flux_error: one-sigma uncertainty on flux, NaN where the model does not depend on the flux.
  • snr: flux / flux_error, the detection significance map.
  • posterior_mean, posterior_sd, log_bayes_factor: the posterior under the given prior and the log evidence ratio against f = 0; None if prior is not given, whatever the kind of prior.

Because there are six fields, unpack by attribute (res.flux) or take the first three with flux, error, snr = res[:3]; unpacking the result directly into three names fails.

Examples:

>>> grid = {
...     "dra": np.linspace(-300.0, 300.0, 61),
...     "ddec": np.linspace(-300.0, 300.0, 61),
...     "flux": np.array([1e-3]),  # ignored: only names the parameter
... }
>>> res = linear_flux_grid(
...     BinaryModelCartesian, data, grid
... )
>>> i, j = np.unravel_index(
...     np.nanargmax(res.snr), res.snr.shape
... )

laplace_flux_uncertainty_grid(model, data, grid, flux=None, *, flux_param=None, batch_size=None)

Laplace uncertainty of the flux at every grid position.

At each position the coordinates are held fixed and the uncertainty is the inverse square root of the curvature of the negative log likelihood along the flux.

Parameters:

Name Type Description Default
model SourceModel or class

Template model or model class, as for :func:likelihood_grid.

 required
data OIData

Data to fit.

 required
grid dict[str, array - like]

Grid axes, as for :func:optimized_flux_grid. The output has one axis per coordinate key (every key except the flux), in this order.

 required
flux array - like

Flux at which to evaluate the curvature, with one axis per coordinate key. By default this is the best fit from :func:optimized_flux_grid, which is also the mean that ruffio_upperlimit expects; pass it if you have already computed it.

 None
flux_param str

The key of grid holding the flux. By default, the one key whose last part is flux.

 None
batch_size int

Number of grid points evaluated at once, by default enough for a fixed number of model visibilities; see likelihood_grid.

 None

Returns:

Type Description
array - like

One-sigma flux uncertainty, with one axis per coordinate key. It is NaN where the curvature is not positive (the flux is not at a likelihood maximum).

best_grid_point(loglike_grid, grid)

Return the grid point with the highest log likelihood.

Parameters:

Name Type Description Default
loglike_grid array - like

Output of likelihood_grid for grid, with one axis per key. NaNs are ignored.

 required
grid dict[str, array - like]

The grid axes used to compute loglike_grid.

 required

Returns:

Type Description
dict[str, float]

{name: value} at the maximum, in the order of grid.