`virgil.limits`

Contrast limits, significance, and conversions between flux ratios (companion/primary), contrasts (primary/companion) and magnitudes. A companion 100 times fainter than the primary has flux 0.01, contrast 100 and Δmag 5.

Contrast limits, significance, and flux/contrast/Δmag conversions.

virgil parameterizes a companion by its flux relative to the primary (companion/primary, so 0.01 for a companion 100 times fainter). Results are usually reported instead as a contrast, primary/companion (100 here), or as a magnitude difference Δmag = 2.5 log10(contrast) (5 mag here), which :func:flux_to_contrast and :func:flux_to_delta_mag compute.

:func:ruffio_upperlimit: Bayesian upper limits with a positive-flux prior (Ruffio et al. 2018).
absil_limits: frequentist limits from the chi-squared ratio to the no-companion model (Absil et al. 2011), using :func:nsigma.
injection_limits: limits by companion injection (Gallenne et al. 2015, as in CANDID): the flux at which an injected companion would be detected at sigma.
:func:radial_profile: azimuthal statistics of a limit map, for contrast curves.

`flux_to_contrast(flux)`

Contrast (primary/companion) of a companion/primary flux ratio.

float(flux_to_contrast(0.01)) 100.0

`contrast_to_flux(contrast)`

Companion/primary flux ratio of a contrast (primary/companion).

`flux_to_delta_mag(flux)`

Magnitude difference 2.5 log10(primary/companion) of a flux ratio.

float(flux_to_delta_mag(0.01)) 5.0

`delta_mag_to_flux(delta_mag)`

Companion/primary flux ratio of a magnitude difference.

`ruffio_upperlimit(mean, sigma, percentile)`

Percentile of a flux posterior truncated to non-negative values.

Following Ruffio et al. (2018, eqn 8), the flux posterior at each position is the Laplace Gaussian N(mean, sigma) of the unconstrained (possibly negative) best-fit flux, truncated to flux >= 0 by the positivity prior. This returns its percentile quantile, so percentile = norm.cdf(2) gives a 2-sigma-equivalent upper limit and 0.16, 0.5, 0.84 give a median and credible interval.

Parameters:

Name	Type	Description	Default
`mean`	`float or array - like`	Unconstrained best-fit flux, e.g. from `optimized_flux_grid`. It may be negative.	required
`sigma`	`float or array - like`	Laplace uncertainty of the flux, e.g. from `laplace_flux_uncertainty_grid`, broadcastable to `mean`.	required
`percentile`	`float or array - like`	Quantile(s) of the truncated posterior to return, between 0 and 1.	required

Returns:

Type	Description
`array - like`	Non-negative flux at each percentile, with shape `broadcast(mean, sigma).shape + percentile.shape`.

Notes

The flat flux >= 0 prior is the published convention (Ruffio et al. 2018) and is deliberately not the Jeffreys prior for a scale: a log-uniform prior on the flux makes the posterior improper as flux -> 0, so there would be no finite upper limit to quote.

The quantile is computed from the upper tail, mean + sigma * z with Q(z) = (1 - percentile) Q(-mean / sigma) and Q the standard normal survival function, evaluated in log space so that it stays accurate when the best fit is many sigma below zero. Where the tail underflows, the large-deviation limit z = sqrt(a**2 - 2 log(1 - percentile)) (a = -mean / sigma) is the starting point, and two Newton steps on log Q(z) polish the result.

`absil_limits(model, data, grid, sigma, *, flux_param=None, flux_bounds=(1e-06, 1.0), batch_size=None)`

Flux above which a companion is ruled out at sigma significance.

Following Absil et al. (2011), at each grid position this finds the flux at which the model fits the data worse than the no-companion model by a chi-squared ratio corresponding to sigma (see nsigma). Brighter companions at that position are excluded at sigma: with sigma=3 the result is a 3-sigma upper limit on the flux.

Parameters:

Name	Type	Description	Default
`model`	`SourceModel or class`	Template model or model class, as for `likelihood_grid`. The no-companion model sets every parameter in `grid` to zero.	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). The flux axis is only used with `flux_bounds=None`, where its smallest positive value starts the search (a single value is enough).	required
`sigma`	`float`	Exclusion significance. It must exceed the significance of a chi-squared ratio of 1 (about 0.67 for many degrees of freedom), and be below the largest significance the floating-point type can represent (about 12.9 in float32, 37 in float64).	required
`flux_param`	`str`	The key of `grid` holding the flux solved for at each grid position. By default, the one key whose last part is `flux`.	`None`
`flux_bounds`	`tuple[float, float] or None`	Search range of the flux (default `(1e-6, 1.0)`), searched upward from its lower end, as in `injection_limits`. Limits outside it are returned at the nearer bound, and a `RuntimeWarning` reports how many. Pass `None` to search without bounds, from the smallest positive value of the flux axis, e.g. for `System` weights that may exceed 1.	`(1e-06, 1.0)`
`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`	Flux limit (companion/primary), with one axis per coordinate key; see :func:`flux_to_contrast` and :func:`flux_to_delta_mag`.

Notes

The number of degrees of freedom is the number of data points; the fitted parameters are not subtracted.

The limit is the first flux, going up from the start of the search, at which the significance reaches sigma: the flux is stepped by decades until it does, and that decade is bisected in log flux. For a normalized scene the significance falls again once the companion outshines the primary (flux well above 1), so an unbounded search should start below that.

`injection_limits(model, data, grid, sigma, *, flux_param=None, flux_bounds=(1e-06, 1.0), batch_size=None)`

Flux at which an injected companion would be detected at sigma.

This is the injection method of Gallenne et al. (2015, section 3.2), as in CANDID's detectionLimit(methods=["injection"]). At each grid position a companion of flux f is added to the data (to every observable, including the extras: the observables of the model with the companion minus those of the no-companion model), and the no-companion model is fitted to the result. The limit is the flux at which the null model fits the injected data worse than the companion model does, by a chi-squared ratio corresponding to sigma (see nsigma). The companion model fits the injected data exactly as the null model fits the original data, so the ratio is chi2(data + signal(f)) / chi2(data) for the null model. (With gains or OI_FLUX data, whose whitening depends on the model, the companion model's chi-squared on the injected data is computed in full.)

Compare absil_limits, which uses chi2(data - signal(f)) / chi2(data). The two differ by the sign of the cross term between the data's residuals and the signal.

Parameters:

Name	Type	Description	Default
`model`	`SourceModel or class`	Template model or model class, as for `likelihood_grid`. The no-companion model sets every parameter in `grid` to zero.	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). The flux axis is only used with `flux_bounds=None`, where its smallest positive value starts the search (a single value is enough).	required
`sigma`	`float`	Detection significance; see `absil_limits` for its range.	required
`flux_param`	`str`	The key of `grid` holding the flux that is solved for. By default, the one key whose last part is `flux`.	`None`
`flux_bounds`	`tuple[float, float] or None`	Search range of the flux (default `(1e-6, 1.0)`), searched upward from its lower end. Limits outside it are returned at the nearer bound, and a `RuntimeWarning` reports how many. Pass `None` to search without bounds, from the smallest positive value of the flux axis, e.g. for `System` weights that may exceed 1.	`(1e-06, 1.0)`
`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`	Flux limit (companion/primary), with one axis per coordinate key; see :func:`flux_to_contrast` and :func:`flux_to_delta_mag`.

Notes

The limit is the first flux, going up from the start of the search, at which the significance reaches sigma, found as in absil_limits. The significance rises with flux once the signal exceeds the noise, but the cross term can make it dip at very faint fluxes, and for a normalized scene it falls again once the companion outshines the primary (flux well above 1), so an unbounded search should start below that.

Differences from CANDID: CANDID refits the primary's diameter to the injected data before computing the null chi-squared (for V² and T3), which virgil does not; the null model is exactly the one with every grid parameter set to zero. CANDID also brackets the limit in steps of 1.4 in flux and interpolates linearly in significance, where virgil solves the criterion by bisection. The chi-squared and the number of degrees of freedom are virgil's, as for absil_limits.

Examples:

>>> import jax.numpy as np
>>> from virgil import PointSource, System, UniformDisk, injection_limits
>>> template = System(star=UniformDisk(0.8), comp=PointSource(0.01))
>>> grid = {
...     "comp.dra": np.linspace(-10, 10, 21),
...     "comp.ddec": np.linspace(-10, 10, 21),
...     "comp.flux": np.array([0.01]),
... }
>>> limits = injection_limits(template, data, grid, 3.0)
>>> limits.shape
(21, 21)

`chi2ppf(p, df)`

Percentile function for chi-square.

It is 2·gammaincinv(df/2, p), inverting jax.scipy.special.gammainc by Halley's method from a Wilson–Hilferty start, with the residual taken in the upper tail (gammaincc) for p > 1/2. It is accurate to about 1e-14 relative in float64 for every df and p in [1e-10, 1 - 1e-10], needs no optional dependency, and is differentiable (by the implicit function theorem) and jit-able in both arguments.

Parameters:

Name	Type	Description	Default
`p`	`array - like`	Percentile value.	required
`df`	`array - like`	Degrees of freedom.	required

Returns:

Type	Description
`array - like`	Corresponding chi2 value to the percentile.

Notes

p is clipped to [eps, 1 - eps] of its own floating-point type, so the result stays finite. Near p = 1 this loses precision; to convert small tail probabilities, use nsigma, which works with the upper tail directly.

`nsigma(chi2r_test, chi2r_true, ndof)`

Convert a reduced-chi-squared ratio to a Gaussian-equivalent significance.

The statistic x = ndof * chi2r_test / chi2r_true is compared with a chi-squared distribution of ndof degrees of freedom, and its upper-tail probability is expressed as the equivalent two-sided Gaussian significance (as in Absil et al. 2011).

Parameters:

Name	Description	Default
`chi2r_test`	Reduced chi-squared of test model.	required
`chi2r_true`	Reduced chi-squared of true model.	required
`ndof`	Number of degrees of freedom.	required

Returns:

Name	Type	Description
`nsigma`	`float`	Detection significance in Gaussian sigma.

Notes

The upper tail is computed directly (with the regularized incomplete gamma function) rather than as 1 - cdf, so significances stay finite and accurate far beyond 8σ, including in float32 (up to about 13σ).

`radial_profile(values, dra, ddec, center=(0.0, 0.0), r_max=None, bins=20)`

Azimuthal statistics of a map in bins of separation.

Parameters:

Name	Type	Description	Default
`values`	`array - like`	Map with shape `(len(dra), len(ddec))` (axis 0 is `dra`), as returned by the grid and limit functions, e.g. flux limits.	required
`dra`	`array - like`	Grid axes in milliarcseconds.	required
`ddec`	`array - like`	Grid axes in milliarcseconds.	required
`center`	`tuple[float, float]`	Centre `(dra, ddec)` of the annuli in milliarcseconds.	`(0.0, 0.0)`
`r_max`	`float`	Outer radius in milliarcseconds (default: the largest separation on the grid).	`None`
`bins`	`int`	Number of annuli.	`20`

Returns:

Type	Description
`dict`	`r` (annulus centres, mas), `mean`, `std`, `median`, `q16`, `q84` and `count` per annulus. Non-finite values are ignored; empty annuli are NaN.