Skip to content

virgil.orbit_search

Scoring candidate orbits on multi-epoch data with the nuisances an orbit shares across epochs shared in the score: one companion flux per band (with an optional chromatic slope) integrated out on a log-uniform grid, each dataset's error scales marginalized, calibration gains and closure-phase offsets profiled, and extra terms such as radial velocities added. Positions are evaluated at every sample's own time.

Scoring candidate orbits on multi-epoch data, with shared nuisances.

The first stage of an automatic orbit search scores many candidate orbits on all the data at once. score_orbits does this exactly, on the visibilities and closure phases, with the nuisances that a single orbit shares across epochs shared in the score:

  • Flux. One companion flux per band (or instrument), the same at every epoch of that band, optionally with a chromatic power law f(λ) = f₀ (λ/λ₀)^β. It is integrated out under a log-uniform (Jeffreys) prior, adaptively about its profiled peak, and its profiled value and uncertainty are reported. Letting each epoch have its own flux is a strictly looser model, whose extra freedom moves each epoch's peaks and makes new ones that no orbit with one flux visits, so it is never used to rank.
  • Error scales. Each dataset's error scale is integrated out block by block, as in marginal_loglike.
  • Calibration gains set with OIData.with_gains and closure-phase offsets set with OIData.with_closure_offsets are truly per epoch. They enter the residuals linearly and are profiled analytically (each costs one degree of freedom), not marginalized, so that a block's χ² stays proportional to 1/s².

The companion's position is evaluated at every sample's own time, so motion within a night needs no threshold. Extra likelihood terms, such as radial velocities (RVData) or published positions (PositionData), are added to the score.

SharedFlux dataclass

The companion flux shared across epochs, per band.

Parameters:

Name Type Description Default
grid (tuple, array - like or dict)

The coarse flux grid at the reference wavelength λ₀: (lo, hi, n) for n points spaced uniformly in ln f from lo to hi, or the flux values themselves (one value fixes the flux), or a dict {band: either}. The default, (1e-3, 1.0, 16), is used for every band not in the dict. The prior is log-uniform (Jeffreys) between the first and last points. The grid only finds the peak: the marginal is integrated adaptively about the profiled peak (see order).

 (0.001, 1.0, 16)
bands dict

{epoch or dataset name: band}; a dataset's own entry overrides its epoch's. Every dataset of a band shares that band's flux. By default all datasets form one band, "all".

 None
reference optional

The band whose flux is at most 1 (default: the first band, in the order of the datasets). Bounding the flux in one band fixes which star is called the companion (a companion at r with flux f is the same in V² and closure phases as one at −r with flux 1/f), so that the labelling is the same in every band; the other bands may be given grids above 1.

 None
slope (tuple, array - like or dict)

The chromatic slope β, as (lo, hi, n) (uniform prior) or the values, or a dict per band. None (the default) fixes β = 0.

 None
wavel0 float or dict

The reference wavelength λ₀ (metres), one or per band; by default the geometric mean of each band's wavelengths.

 None
order int

Gauss–Legendre nodes per side of the adaptive marginal (default 16): each side of the profiled peak is integrated out to where the score has fallen by 12 nats (grown from the curvature's length scale, and at a prior bound also the slope's), cut at the prior's bounds, with the nodes crowded towards the peak. With a slope, order nodes per side over the whole β prior each integrate ln f from their own conditional peak.

 16
newton int

Newton steps from the best coarse grid point to the profiled peak (default 8), each at most one grid step and kept only if it raises the score.

 8

grid = (0.001, 1.0, 16) class-attribute instance-attribute

bands = None class-attribute instance-attribute

reference = None class-attribute instance-attribute

slope = None class-attribute instance-attribute

wavel0 = None class-attribute instance-attribute

order = 16 class-attribute instance-attribute

newton = 8 class-attribute instance-attribute

__init__(grid=(0.001, 1.0, 16), bands=None, reference=None, slope=None, wavel0=None, order=16, newton=8)

__post_init__()

OrbitScores dataclass

The scores of candidate orbits, in the order they were given.

Returned by score_orbits. Rank them with :meth:order (or rank_scores).

Attributes:

Name Type Description
score ndarray

(n,): the log likelihood with every shared flux (and slope) integrated out, error scales marginalized (or on the quoted errors), gains profiled, plus the extra terms. -inf where not finite.
profiled ndarray

(n,): the same at the profiled flux (and slope) of each band, instead of integrated over them.
flux, flux_err ndarray

(n, n_band): the profiled flux f₀ of each band (Newton steps in ln f from the best coarse grid point), and its uncertainty from the curvature in ln f there, σ_f = f σ_ln f (NaN with one grid point, or where the curvature is not negative).
flux_at_edge ndarray

(n, n_band) bool: the profiled flux is at a prior bound.
slope, slope_err, slope_at_edge ndarray

The same for the chromatic slope β (0 and NaN without a slope grid). With both free, the uncertainties are marginal (from the inverse of the 2 × 2 curvature).
fallback ndarray

(n, n_band) bool: the curvature at the peak was not negative or not finite, so the marginal is the coarse grid's trapezoidal rule instead of the adaptive integral.
terms ndarray

(n, n_terms): each extra term's log likelihood (included in score).
scale ndarray

(n, n_dataset): each dataset's largest block error scale ŝ = √(χ²/ν) at the profiled flux (and slope), on the quoted errors.
scale_at_bound ndarray

(n, n_dataset) bool: a block's scale posterior presses against s_max (ŝ within one posterior width in ln s, 1/√(2 dof ν), of the bound), so its score is penalized by the bound: the candidate needs errors inflated past s_max there. Always False with s_max=None.
bands tuple

The band names, in column order; reference is the one whose flux is at most 1.
reference object
cost int

The work units spent: candidates × Σ over datasets of the score evaluations of the band's marginal (grid points, Newton steps with their derivatives, and nodes).

score instance-attribute

profiled instance-attribute

flux instance-attribute

flux_err instance-attribute

flux_at_edge instance-attribute

slope instance-attribute

slope_err instance-attribute

slope_at_edge instance-attribute

fallback instance-attribute

terms instance-attribute

scale instance-attribute

scale_at_bound instance-attribute

bands instance-attribute

reference instance-attribute

cost instance-attribute

__init__(score, profiled, flux, flux_err, flux_at_edge, slope, slope_err, slope_at_edge, fallback, terms, scale, scale_at_bound, bands, reference, cost)

__len__()

order(quantum=_QUANTUM)

Candidate indices, best first, ties broken by index (see rank_scores).

score_orbits(epochs, model, orbits, *, shared=None, terms=(), scales='marginal', dof=1.0, s_max=5.0, batch_size=None, max_evaluations=None, dtype='float64')

Score candidate orbits on all epochs, with the flux shared.

For each candidate, the companion's position is computed at every sample's own time, and the visibilities g of a unit companion there once. Every evaluation at another flux is then cheap arithmetic on g: for a scene linear in the companion's flux f (a component weight), the complex visibility is V(f) = (A + f B) / (T₀ + f ΔT), with A and the total fluxes computed once per dataset and B from g. With a slope, f differs per sample: f = f₀ (λ/λ₀)^β.

For each dataset and flux the score is the scale-marginalized log likelihood m = -Σ_b (ν_b/2) ln χ²_b of marginal_loglike (or, with s_max, its bounded form; with scales="quoted", -χ²/2). Gains (OIData.with_gains) and closure-phase offsets (OIData.with_closure_offsets) are profiled analytically inside χ²: their widths are ignored, and each independent mode removes one degree of freedom from its block (a single gain per dataset is with_gains(modes=onp.ones((1, n))), with n the number of samples).

Each band's flux is then integrated out under its log-uniform prior. The coarse grid only finds the peak: Newton steps in ln f (and β) profile it, and Gauss–Legendre nodes on each side of it, out to a 12-nat drop or the prior's bounds, integrate it, so that the marginal does not depend on where the peak falls between grid points (a well-measured flux, σ_ln f ≈ 0.02, is far narrower than any affordable grid step). Where the curvature at the peak is not usable, the coarse grid's trapezoidal rule is used and OrbitScores.fallback is set. The extra terms are then added.

Parameters:

Name Type Description Default
epochs Epochs

The data. Datasets without sample times (mjd) are evaluated at their snapshot time (Epochs.times).

 required
model callable

model(orbit, flux) returns the time-dependent scene, e.g. OrbitalBinary itself. It must be linear in flux (a component weight), and the scene at flux 0 must not depend on the orbit; both are checked on the first and last candidates.

 required
orbits sequence or orbit

A list of orbits (or starting_orbits pairs), or one orbit whose elements have a leading candidate axis. All share one class and t_ref; times are measured from t_ref, in float64 on the host, so float32 keeps them precise.

 required
shared SharedFlux

The shared flux grids, bands and slopes (default SharedFlux(): one band, f in [1e-3, 1] on 16 log-spaced points, no slope).

 None
terms sequence

Extra log likelihood terms: likelihood terms from RVData.term or PositionData.term, whose params (or orbit) callable receives {"orbit": orbit}, or callables orbit -> log likelihood. Each brings its own error model (e.g. a fixed RV jitter).

 ()
scales (marginal, quoted)

Integrate each dataset's error scales out (default), or use the quoted errors (the score is then -χ²/2, up to a constant).

 "marginal"
dof optional

The effective-dof fraction (one number, or a dict keyed by dataset or epoch name) and the bound on the scales, as for marginal_loglike (each block's likelihood integrated over ln s in [-ln s_max, ln s_max], log-uniform). The default s_max=5 leaves about twice the error inflation of typical interferometric data: a candidate that needs more than that on some dataset is penalized, and flagged in OrbitScores.scale_at_bound. The bound also keeps the score finite: with s_max=None a block's -(ν/2) ln χ² is unbounded as χ² → 0, so a block with few degrees of freedom (one triangle's closure phases, say) that a flux or slope fits almost exactly makes a narrow spike in the score, which can dominate the profile (profiled, flux, slope). A finite s_max bounds each block's gain at about ν ln s_max.

 1.0
s_max optional

The effective-dof fraction (one number, or a dict keyed by dataset or epoch name) and the bound on the scales, as for marginal_loglike (each block's likelihood integrated over ln s in [-ln s_max, ln s_max], log-uniform). The default s_max=5 leaves about twice the error inflation of typical interferometric data: a candidate that needs more than that on some dataset is penalized, and flagged in OrbitScores.scale_at_bound. The bound also keeps the score finite: with s_max=None a block's -(ν/2) ln χ² is unbounded as χ² → 0, so a block with few degrees of freedom (one triangle's closure phases, say) that a flux or slope fits almost exactly makes a narrow spike in the score, which can dominate the profile (profiled, flux, slope). A finite s_max bounds each block's gain at about ν ln s_max.

 1.0
batch_size int

Candidates evaluated together by jax.lax.map (default: all). The result does not depend on it.

 None
max_evaluations int

A budget in work units: candidates × Σ over datasets of the score evaluations of the dataset's band's marginal (coarse grid points, Newton steps counting each derivative as one, and nodes). The cost is predicted before anything is evaluated, and a cost over budget raises a ValueError.

 None
dtype (float64, float32)

Precision of the evaluation (default float64, in a local jax.enable_x64 context).

 "float64"

Returns:

Type Description
OrbitScores

Scores in the order of orbits; rank them with OrbitScores.order().

Examples:

>>> scores = score_orbits(epochs, OrbitalBinary, candidates,
...                       shared=SharedFlux((0.01, 1.0, 16)))
>>> best = scores.order()[0]
>>> scores.flux[best], scores.flux_err[best]

rank_scores(score, quantum=_QUANTUM)

Indices that sort score from best to worst, deterministically.

Scores are first rounded to multiples of quantum (nats), and equal rounded scores rank by index, so that two candidates differing only by rounding error (e.g. mirror-image grid cells) rank the same way on every platform and in either precision. Non-finite scores rank last.

Parameters:

Name Type Description Default
score array - like

One score per candidate, higher is better.

 required
quantum float

The rounding step (default 1e-3).

 _QUANTUM

Returns:

Type Description
ndarray

Integer indices, best first.

Examples:

>>> rank_scores([1.0, 3.0, 3.0 + 1e-9, -float("inf"), 2.0]).tolist()
[1, 2, 4, 0, 3]