Pipelines

A tutorial shows how to do an analysis; a pipeline does it the same way every time and leaves a folder that a person, a script or an agent can read without rerunning the code. virgil.pipeline wraps virgil's public functions in fixed, named stages. Each stage writes its own outputs, an interrupted run resumes, and every run ends with deterministic quality checks and an executed quicklook notebook.

A pipeline adds no science: every number it reports comes from a public virgil function that a tutorial also calls, so you can reproduce any result by hand. The binary-search, contrast-limits and limb-darkening tutorials end with an "As a pipeline" section that does exactly that and asserts that the numbers agree; the orbit tutorial leads with its pipeline and then shows the functions underneath.

Install the extra (it brings in the plotting extra, h5py, nbformat, nbclient and ipykernel):

pip install 'virgil-astro[pipeline]'

import virgil never imports the pipelines; virgil.pipeline loads on first use.

Running a pipeline

import virgil as vg
from virgil.pipeline import BinaryPipeline, load

data = vg.OIData("hd1234.oifits")
res = BinaryPipeline(data, output="runs/hd1234", sigma=3.0).run()

res.summary["companion"]["sep_mas"]
res.checks            # a list of Check(name, status, value, threshold, message)
print(res.describe())

res = load("runs/hd1234")   # later, in a fresh session; nothing is recomputed
model = res.model()         # the best-fit model, a real virgil object

A pipeline takes the data, a model template, a run folder and keyword settings. Settings resolve as the class defaults (BinaryPipeline.defaults()), then a JSON config file (BinaryPipeline.from_config(path)), then keyword arguments. run(through="limits") stops after a stage, which is how a non-detection skips the fit and the sampler.

A run resumes by default: a completed stage whose outputs exist is skipped, and once a stage reruns every later stage reruns. A resumed run raises ConfigMismatchError if the settings, the data hash or the model template differ from those stored in the folder. run(resume=False) starts again.

The stages of BinaryPipeline

The fit is on the quoted errors unless error_scale="fit", and chi-squared per independent datum is always reported on the quoted errors. Priors are group-invariant: uniform in position (dra and ddec on a box of half-width max_sep_mas) and log-uniform in flux and in error scales. Both templates sample the same dra, ddec and flux, so they give the same prior on the sky; a BinaryModelAngular template reports sep and pa derived from them, with no wrap at 0 or 360 degrees in the fit. The flux range is a setting (flux_range).

max_sep_mas is by default the coarsest λ/B. Where the uv samples of one wavelength lie on a lattice (aperture masks, DISCO products), the visibilities repeat when the companion moves by a vector of the reciprocal lattice, so the default box is limited to half an alias cell, the largest square with no two points that differ by such a vector. The summary records the limit as search.alias_cap_mas (null when it was not applied), and the quicklook notes it. An explicit max_sep_mas is never capped; OrbitPipeline limits its position grid in the same way.

BinaryPipeline settings

BinaryPipeline.defaults() lists the values:

The search and the limits always use a point companion in Cartesian offsets; the model template (BinaryModelAngular or BinaryModelCartesian) sets only the parameters of the fit and the posterior.

StarPipeline: stellar diameters

from virgil.pipeline import StarPipeline

res = StarPipeline(data, output="runs/star").run()   # uniform and limb-darkened, compared
res = StarPipeline(data, "uniform", output="runs/ud").run()
res.summary["star"]["diam_mas"], res.model(), res.model("uniform")

model= is a short name, "uniform" (UniformDisk) or "limb_darkened" (QuadraticLimbDarkenedDisk with Kipping's q1, q2), or an instance of one of those classes (its q1 and q2 are the starting values). The default fits both and compares them. Each fitted model is saved under models/<name>/, with models/best the preferred one: the limb-darkened model only if its delta chi-squared exceeds the BIC penalty of its two extra parameters. res.model(), res.model_values(), res.fit_result(), res.samples() and res.sample_stats() take an optional name= to pick one.

Priors are group-invariant: log-uniform in the diameter over diam_range_mas (by default from a tenth of the resolution limit to four times the field of view) and uniform on [0, 1] in q1 and q2, which cover exactly the physical quadratic laws. The fit is on the quoted errors, and chi-squared per independent datum is reported on them for each model.

StarPipeline settings

ImagingPipeline: a model-independent image (provisional, first version)

from virgil.pipeline import ImagingPipeline

res = ImagingPipeline(data, output="runs/image").run()
img = res.image()          # a SkyImage: pixels, pixel_scale_mas, dra, ddec, beam

ImagingPipeline reconstructs a pixel image with a Gaussian-process prior whose amplitude σ and length ℓ are chosen by evidence (star="auto" also decides, by ΔlogZ against a threshold of 5 nats, whether a point source sits at the centre). The stages are load, overview, star, hyperparameters, members (default 8 fits at the shared hyperparameters, resumable one by one with run_member(k)), combine and quicklook. The result is images/mean.fits (a sky-oriented FITS with a WCS and the star's share of the flux) and images/sigma.fits (the spread of the members, which excludes hyperparameter uncertainty), models/best (res.model()) and a quicklook that plots the image.

This first version fits an ensemble of up to 8 members (members, default 8), which by default share the hyperparameters and differ only in start, pixel scale and field; ensemble.combine keeps those that fit V² and closure phases, and sigma.fits is their spread, which with shared hyperparameters leaves out hyperparameter uncertainty (member_grids=True gives each member its own). CLEAN starts, classical regularizers chosen by L-curve, the anneal option and per-channel imaging are to come.

virgil-pipeline image star.oifits -o runs/image --set star=auto

ImagingPipeline settings

DiskPipeline: rings and haloes (provisional)

from virgil.pipeline import DiskPipeline

res = DiskPipeline(data, output="runs/disk").run()                 # a Gaussian-process disk
res = DiskPipeline(data, "rings", output="runs/rings", companion=True).run()
res.summary["disk"], res.model(), res.model("geometry")

DiskPipeline fits a MultipoleDisk around one star. It is provisional: the settings, stages and summary keys may change in any 0.4.x release without a deprecation. model= picks one of two parallel options, never a chain:

• "gp" (the default, the option that assumes least about the disk): one RadialGP profile with low-degree harmonics (harmonic_degree) and, with more than spectral_index_min_channels channels, a spectral index that is itself a Gaussian process. At each length scale of ell_halffringes the profile is fitted with the geometry held at the single-ring fit's values (fast), then again with the inclination and position angle released, and the length of largest Laplace evidence is kept. A grid must not choose its own edge: if the best length is the smallest or the largest, the grid is extended once in that direction (never below ell_floor_halffringes, half a half-fringe), and the length_scale check warns if the best is still on an edge (on_grid_edge in the summary). With error_scale="rescale" (the default here) one MacKay error scale is estimated from a baseline fit at the middle length, with the companion if there is one and the geometry released, and every fit and evidence after it is on the rescaled data, so they share one footing; the scale is reported (chi2.error_scale_factor), and χ²/N is still quoted on the quoted errors.
• "rings": N GaussianRing profiles that share the inclination and position angle, added one at a time (each starting at the radius of a one-ring scan) until the BIC stops improving or max_rings is reached; after each fit every ring is re-placed in turn (dropped, rescanned on a finer radius grid, everything refitted) and the result kept if the BIC falls (not if it leaves two rings closer than twice the narrower one's width or a half-fringe). A count with such a pair is one ring fitted as two, or a ring still to be found: it is left out of the comparison and the search goes on; a ring whose removal lowers the BIC (one parked at the inner radius with no light, say) is dropped after the fit, which ends the search, the pruned fit standing as the candidate of its own count; ring_count warns when such a count had the lower BIC. The Laplace evidence and the BIC of every count are reported and the count of least BIC is kept.

The resolution scale is the half-fringe λ/2B of the longest baseline. The inner radius, the ring-width floor and the length-scale grid are stated in half-fringes, and the outer radius defaults to half the image_field: set r_out_mas if the disk is larger. A companion search is opt-in (companion=True): the residuals are searched for a point source on a grid, and a detection above companion_min_delta_chi2 is fitted jointly with the disk. A search of the best model's residuals always runs, to check that the disk has not absorbed a compact source.

Priors are group-invariant: uniform in cos i, an AngleVector for the position angle (no wrap), log-uniform in ring radii, widths and fluxes, and an isotropic Gaussian (scale 0.4) on the real and imaginary parts of each ring's harmonics, so that the phase of the lopsidedness is uniform. The rings are fitted on the quoted errors (error_scale="quoted", the default there), and error_scale="fit" is available for them only; "rescale" is for the Gaussian process (see above). The star is a grey point source unless star_diameter_mas (a uniform disk; the inner radius is then at least its radius), star_temperature_k or star_index says otherwise, and the disk's spectral index is relative to that star: it is conditional on the star's assumed spectrum, so read it with the model.star section of the summary beside it (on real discs a different star index moved the disk's index by 0.3 to 1.4). The rings' colour is a power-law index per ring by default; ring_colour="blackbody" (rings only) fits a temperature per ring and for the planet instead, which follows blackbody-like spectra across a broad band that a power law misses. A model saved under models/ is a real virgil object: res.model() is the preferred one, res.model("geometry") the single-ring fit and res.model("rings_2") or res.model("gp_ell1") a candidate.

DiskPipeline settings

OrbitPipeline: a binary seen on several nights (provisional)

OrbitPipeline fits the orbit of a binary observed on several nights, and it does not need to be told the period. The default run is a blind search over the period aliases: the window of periods comes from the data, and no published period, warm start or paper prior is used. Give it the OIFITS files (each is split into nights at gaps of gap_days), or an Epochs, or {night: data}:

from virgil.pipeline import OrbitPipeline

res = OrbitPipeline.from_oifits(files, output="runs/hd1234").run()

res.summary["modes"]["winner"]      # the winning alias band and its orbit
res.summary["modes"]["delta_log_z"] # its margin over the runner-up, in nats
res.aliases().table()               # every band: log Z, probability, flags

The orbit tutorial runs this on synthetic two-season data, reads the result and then calls the functions underneath (epoch_positions, fit_orbit_aliases) for the control the pipeline's settings do not give.

OrbitPipeline takes data that the instrument package has already prepared (the channel, band edges, outlier exposures and noise model): it neither selects a channel, trims the band nor screens exposures, and it has no instrument-specific settings, so it works the same for every instrument. Per-file rulings go in metadata (epoch, and mjd for a file whose time stamp is wrong or absent) and exclude ({file, night or dataset: reason}, recorded in run.json). The command line is virgil-pipeline orbit FILES... --output DIR [--metadata FILE.json].

A run that interrupts resumes: the bands already fitted are read back from their checkpoints, so only the missing bands are refitted. run(resume=False) starts over, bands included. The winner is the band of highest evidence; the summary's delta_log_z is its log Z over the runner-up's, and alias_ambiguity warns below 10 nats and fails below 1. no-orbit compares the best orbit's evidence (with its prior volume) with the best log likelihood of the star alone (a maximum, not an evidence, so the test leans towards no orbit), and a run that finds no orbit skips the posterior. The posterior is sampled with NUTS (dense mass matrix) on the engine's own tempered, scale-marginalized density, not with gauss_newton_mass: that builds its curvature from the Jacobian of a virgil model's whitened residuals, which this density (tempered, with the V² amplitude profiled) is not. On the two-season synthetic case a diagonal mass matrix gave a similar minimum effective sample size (225 against 156, 300 warm-up and 300 draws on 2 chains) at 60% of the cost, so the dense mass is a conservative choice for the narrow ridges that real, correlated orbits can have.

OrbitPipeline settings

The run folder

Everything is JSON, HDF5, NumPy .npz, PNG, OIFITS or one notebook, with no pickles. Files are written atomically, so a crash never leaves a half-written file under its final name, and nothing is written outside the run folder.

Reading a run: load and Result

load(path) returns a Result. Nothing is recomputed: each accessor reads the files the run wrote.

Checks

Checks are pure functions of plain numbers, so the same run always gives the same checks in a fixed order, each with a status of pass, warn or fail. The first is always the chi-squared per independent datum of the model on the quoted errors: a fitted error scale makes chi-squared close to one by construction, so it is never the headline number.

Warnings raised inside stages

A stage runs under warnings.catch_warnings(record=True) with the filter set to "always", so a warning raised inside a stage never reaches the terminal or the notebook. Each distinct message is recorded once per stage, with its stage, category and message, in stages/<name>/report.json and in the warnings list of summary.json; describe() counts them. The known messages are also turned into checks, matched by keywords so that rewording a message does not break the mapping. These checks appear only when the warning occurred, after the checks above, and always warn:

Deprecation notices are recorded in warnings but raise no check. The quicklook notebook's own kernel runs on the CPU (JAX_PLATFORMS=cpu, no memory preallocation), so it does not compete with the run for a GPU, and its stderr goes to quicklook.log.

Stability

The pipelines are provisional in 0.4: their constructor keywords, settings and defaults, stage names, summary.json keys, check names, thresholds and order, and file names may change in any 0.4.x release, without a deprecation. They are promoted to stable after a release of use. The tier is the class attribute STABILITY ("provisional"), written to run.json and summary.json; nothing warns at runtime. See Stability for the tiers.

What does not change under the provisional tier:

• the folder layout is versioned by its schema, virgil-pipeline-run-v1, and load refuses schemas it does not know;
• every number a pipeline reports comes from a public virgil function, so a result can always be reproduced by hand.

When a pipeline is promoted to stable, it keeps the following within a schema version:

• the constructor and run signatures, the stage names and their order, and the setting names and defaults;
• every key of summary.json and run.json, and its meaning;
• the check names, their order and their thresholds;
• the file names in the run folder.

After that, changes are additive: new settings whose defaults keep the old behaviour, new summary keys, new checks appended at the end, and new files. Anything else needs a new schema version. Snapshot tests already pin the signatures, defaults, summary keys and check order, so a change is visible in review; while a pipeline is provisional, updating the snapshot is allowed.

Open before promotion: global_nsigma is a Šidák estimate with a nominal number of resolution elements, not a simulated false-alarm probability (the grid maxima are correlated, and the resolution element is anisotropic). It must be calibrated by simulation, or dropped from the summary keys and the detection check, before the pipelines become stable. Until then it is labelled an approximation in summary.json (global_nsigma_method).

Further pipelines will be added to this table as they land, all provisional.

Command line

virgil-pipeline runs and inspects pipelines from a shell.

virgil-pipeline binary hd1234.oifits -o runs/hd1234 --set sigma=3 --set num_chains=2
virgil-pipeline binary hd1234.oifits -o runs/hd1234 --through limits --fresh
virgil-pipeline star star.oifits -o runs/star --model limb_darkened
virgil-pipeline disk disk.oifits -o runs/disk --model rings --set companion=true --set r_out_mas=25
virgil-pipeline orbit night*.oifits -o runs/orbit --metadata rulings.json
virgil-pipeline info runs/hd1234        # status and checks
virgil-pipeline validate runs/hd1234    # complete, and every file matches its hash
virgil-pipeline quicklook runs/hd1234   # rebuild the quicklook notebook

--set KEY=VALUE overrides a setting (the value is parsed as JSON where possible). info and validate read only JSON and hashes and do no JAX work, although import virgil still imports JAX.

API reference

The classes and functions are documented on the pipeline API page.