API Documentation

@variables begin
    [prior_1] ~ [UnivariateDistribution]
    [prior_2] ~ [UnivariateDistribution]
    calculation_3 = obs.[prior_1] + obs.[prior_2]
end

System(; name, bodies, observations=(), variables)

A model of one system: its bodies (and any explicit Orbit rows), its observations, and its system-level variables.

sys = System(
    name="HD 12345",
    bodies=[A, b],
    observations=[astrom, rvs],
    variables=@variables begin
        plx ~ Uniform(1, 100)
    end
)

Frame

Which of plx, ra, dec, pmra, pmdec, rv, ref_epoch the system block defines chooses the frame: none gives physical-unit observables only, plx alone gives angular observables in mas, and the full set gives a rigorously propagated absolute frame. There is no basis= keyword.

A system may contain a single Body and no orbits: the lone body is the system barycentre, and its motion is purely the frame's. That is the natural model for fitting the absolute astrometry of an isolated star (e.g. the first HipparcosIADObs tutorial) — no zero-mass companion needed.

Deferred lines

A system line that mentions a body by name (mut_inc = b.i - c.i) is evaluated after every body block, so couplings, period ratios and global LL += constraints need no new syntax. Bodies look up and outward; the system block looks down and inward, after the fact; siblings never see each other.

Deferred lines also see system_interim: a PlanetOrbits.System built from this model's own bodies, which they may solve. It carries no absolute frame — which is the point, since what deferred lines most usefully compute is the frame:

plx_interim = plx_A                                          # not deferred
Δ    = anchor_offsets(system_interim, :A, 57388.0)           # deferred
pmra = pmra_A - Δ.pmra                                       # …and so is this

Body–barycentre kinematics need a distance, not a sky position, so the interim is complete at Parallax level and there is no circularity. Give it that distance with a non-deferred plx_interim (or plx); without one it is frame-free and only the physical-unit observables (posx, velx, …) work. AnchoredFrame is the packaged form of this pattern. system_interim is visible only in deferred system lines: a body block cannot see it (it is built from the bodies), and an observation block must not, because observations are evaluated against the final system rather than this one. Reading it anywhere else is an error naming the block.

Observations

All observations live here. Each one names its own references (target=, ref=), so there is nothing left for per-companion attachment to express.

Corrections

observing_geometry and barycentric_lighttime gate PlanetOrbits' two source-side precision corrections. Each takes :auto (the default), :on or :off — true/false are accepted as aliases for the latter two.

:auto measures rather than guesses: at build, it draws ~300 parameter sets from the priors, solves each with the correction on and off, and compares each observation's own model predictions against that observation's own tightest uncertainty. The correction is declined only if the resulting accumulated bias — per-point impact × √(number of data points), since these corrections are common-mode and so grow as √n rather than averaging down — stays under 0.1σ for every observation. Needs are OR'd, so one µas dataset keeps it on for the whole model. A model in which no observation type can report comparable predictions resolves :on without taking a single draw.

The decision is taken once, here — never per draw, which would make the log density discontinuous — and recorded on the built system as sys.corrections, from where it travels into chain metadata. Tune with correction_draws, correction_seed, correction_threshold, and silence the build log with verbosity=0.

einstein_rv (:on by default, or :off) is the third and is not tri-state: it chooses whether radial velocities are predicted from the spectroscopic radvel or the kinematic velz. That is a counterfactual, for measuring what the Einstein term does to a posterior — not something the data can rule out, and not a provenance declaration, since no pipeline can have removed the part of the term that varies with the orbit.

See the "How Octofitter Computes Orbits" manual page, and recheck_corrections for the after-sampling check.

Body(; name, about=nothing, variables)

One gravitating body in the model, and — when about is given — the orbit that places it. about takes the same grammar PlanetOrbits.Orbit uses:

A = Body(name="A", variables=@variables begin mass ~ … end)          # root
b = Body(name="b", about=A,      variables=@variables begin … end)   # astrocentric
c = Body(name="c", about=(A, b), variables=@variables begin … end)   # Jacobi

Exactly one body may omit about; it is the root of the hierarchy. The convention is never inferred from semi-major axis — about=A and about=(A, b) are materially different models under KeplerianApprox.

Variables

Names in the block are read by position in three groups:

• mass [M⊙] — the body's mass. Zero is allowed (a test particle).
• flux, or flux_<band> — per-band flux, used to form photocentres. Setting the host's flux to 1.0 makes every other body's flux a contrast ratio.
• orbital elements — any keyword PlanetOrbits.Orbit accepts (a or P; e/ω or secosω/sesinω or ecosω/esinω; i, Ω; tp or M0+epoch or θ+epoch; or a full Cartesian state). Supplying the wrong combination is a mechanical error from the constructor.

Every other name is an ordinary local, free to be used by later lines in the same block (e.g. mass_jup ~ LogUniform(…) then mass = mass_jup * mjup2msol). A body block sees system.*; it does not see its siblings.

Orbit(; name, exterior, about, variables)

One hierarchy row whose exterior is not a single body — the case Body's fused about= cannot express:

Orbit(name="AB", exterior=(Ba, Bb), about=(Aa, Ab),
      variables=@variables begin P ~ …; e ~ … end)

name gives the row its own chain columns (AB_P, AB_e), which is how binary-star people talk about that orbit. Its variables are read exactly as a Body's orbital elements are; it has no mass or flux of its own, since every mass in the row already belongs to a body.

Barycentre
Barycentre(A, b, …)

The mass-weighted barycentre of the whole system, or of the subsystem spanned by the given bodies. Use as an observation's target or ref:

RadialVelocityObs(rvdata; target=A, ref=Barycentre, …)

Why absolute astrometry also says ref=Barycentre

ref names the point an observation's modelled offsets are measured from. It is not a claim that the point is at rest. The system barycentre does move across the sky in general — and that motion is exactly what the absolute frame in the system block describes: ra, dec, plx, pmra, pmdec and rv at ref_epoch are the barycentre's own catalogue quantities, which PlanetOrbits propagates rigorously in 3D, so its track carries perspective acceleration and light-travel curvature rather than being a straight line in (α, δ).

An absolute-astrometry prediction is the sum of the two:

absolute position = the barycentre's propagated track   # from the frame
                  + (target − ref) at this epoch        # from `target`/`ref`

with the annual parallax supplied by the observation's own parallax factors. G23HObs composes them literally — frame_pmra(sol) + Δpmra, where Δpmra is the five-parameter refit of raoff(sol, target, Barycentre) over the transits. So ref=Barycentre says "this source's excursion is measured from the system's centre of mass"; the centre of mass moving is not an omission, it is the frame's job.

Naming a sub-barycentre is legal and is usually not what you want for a catalogue source in a multi-source system: ref=Barycentre(Aa, Ab) measures the A pair's photocentre against its own centre of mass, which removes that pair's motion about the system barycentre from the prediction — and for a wide pair that motion is the signal. Both sources take the whole-system ref=Barycentre, because both are predicted from the one shared frame; see the G23H tutorial.

Photocentre
Photocentre(:G)
Photocentre(:G, (Aa, Ab))
Photocentre((Aa, Ab))

The flux-weighted photocentre of the system — what a blended astrometric measurement actually tracks — or, given a member list, of that subset only. Bodies carry per-band fluxes (a flux or flux_<band> variable in their block); pass the band when more than one is defined.

A member list says this source blends exactly these bodies, in every draw. Use it when membership is structurally certain — two pairs several arcseconds apart, say, where only intra-pair blending is ever possible, and each catalog source is its own observation:

GaiaDR4AstromObs(dataA; target=Photocentre(:G, (Aa, Ab)), ref=Barycentre, name="srcA")
GaiaDR4AstromObs(dataB; target=Photocentre(:G, (Ba, Bb)), ref=Barycentre, name="srcB")

Member names are validated against the system's bodies at model-build time, and the whole spec lives in type parameters, so it constant-folds. A subset that is dark in the selected band is an error — a structural membership declaration over bodies with no flux is not a point on the sky.

Membership that varies per draw (a sampled resolved-flag) or per epoch (a scan-angle-dependent window) is not this: it belongs to the observation, which reads PlanetOrbits.fluxes(sys, band) and builds its own WeightedPoint. See the PlanetOrbits "Blended sources & photocentres" docs.

UniformCircular(domain=2π)

A variable on a continuous periodic domain. Introduces two normal variables <name>x, <name>y and derives the angle as atan(y, x), so the sampler can wrap freely instead of hitting a hard boundary.

@variables begin
    Ω ~ UniformCircular()
end

Sine()

A custom univariate distribution. The pdf is a sine function defined between 0 and π. This is a common prior distribution used when fitting orbits to astrometry.

The full Distributions.jl interface is not yet defined for this distribution, but the following methods work: pdf, logpdf, minimum, maximum, insupport, mean, var, cdf, quantile

kde = KDEDist(data)

A univariate distribution that obeys the Distributions.jl interface. Uses KernelDensity.jl to create a 1D kernel density estimator using the provided input data and optional bandwidth scale factor.

Appropriate to use as a prior in an Octofitter model.

AnchoredFrame(anchor; ref_epoch, variables, ra=true, dec=true, plx=true, pmra=true, pmdec=true, rv=true)

A variables= block for System in which the absolute frame is parameterized by an anchor source's observed catalogue solution rather than by the system barycentre's.

The frame still means the barycentre — nothing downstream changes — but what the sampler moves is the anchor's own (ra, dec, plx, pmra, pmdec, rv), with the model's motion of the anchor body about the barycentre subtracted to recover the barycentric values. That subtraction reads the interim system, so it sees the same bodies the final system is built from.

A = Body(name="A", variables=@variables begin mass ~ truncated(Normal(1.29, 0.02), lower=0) end)
b = Body(name="b", about=A, variables=@variables begin … end)
B = Body(name="B", about=A, variables=@variables begin … end)   # the wide companion

sys = System(
    name="ups And", bodies=[A, b, B], observations=[g23h_A, g23h_B],
    variables=AnchoredFrame(A; ref_epoch=mjd("2016-01-01"), variables=@variables begin
        ra_A     ~ Normal(24.19928, 1/3.6e6)
        dec_A    ~ Normal(41.40546, 1/3.6e6)
        plx_A    ~ Normal(74.19, 0.20)
        pmra_A   ~ Normal(-172.57, 0.05)
        pmdec_A  ~ Normal(-381.32, 0.04)
        rv_A     ~ Normal(-26900.0, 500.0)
    end))

The sampled names default to <quantity>_<anchor>; pass a Symbol to any of the six keywords to rename one, or false to leave that quantity unanchored — in which case your block must define the barycentric ra / dec / … itself, in the ordinary way. Mixing is legal and is the point of the per-quantity keywords: anchoring proper motion, where the degeneracy is, while sampling the barycentric parallax directly is a perfectly reasonable model.

ref_epoch [MJD] is both the frame's reference epoch and the epoch the anchor's motion is evaluated at, because that is what makes the two consistent: the catalogue solution being anchored to is itself referred to that epoch. Anchoring to several catalogue epochs is not a thing — one frame has one reference epoch, and a source's motion between catalogue epochs is what the observations are for.

What this expands to

Ordinary system lines. Nothing here is privileged, Base.show on the system prints them, and a model that wants something slightly different should write them out rather than reach for another keyword:

plx_interim = plx_A                                  # interim scale
_frame      = anchored_frame(system_interim, :A, ref_epoch,   # deferred
                 (true, true, true, true, true, true);
                 ra=ra_A, dec=dec_A, plx=plx_A,
                 pmra=pmra_A, pmdec=pmdec_A, rv=rv_A)
ra    = _frame.ra
dec   = _frame.dec
plx   = _frame.plx
pmra  = _frame.pmra
pmdec = _frame.pmdec
rv    = _frame.rv

_frame mentions system_interim, so it is deferred automatically, and so is everything that reads it. plx_interim is not: it has to be known before the bodies are, which is why the interim is built at the anchor's parallax rather than the barycentre's — a distinction worth ~1 part in 10⁵ (see barycentre_parallax).

A quantity passed false is left out of the tuple's correct mask and its barycentric value, which your block must define, is passed through the map untouched. All six are still needed, because the map's second pass builds a frame out of them; see anchored_frame.

Priors and the Jacobian

The map (anchored) → (barycentric) is triangular with unit diagonal in ra, dec, pmra, pmdec and rv: each barycentric quantity is its anchored counterpart minus a correction built from the body variables and plx_A only. So a prior declared on pmra_A transfers to pmra with volume factor 1, and needs no LL += correction.

That is exact for the one-pass map (refine=false). The default second pass — which is what makes pmra_A mean the anchor's proper motion to better than 10⁻⁶ mas/yr rather than 4 × 10⁻³ — couples each output weakly to all six inputs, perturbing the diagonal by 10⁻⁹–10⁻⁴ and the determinant by 10⁻⁹. That is a slowly-varying tilt five orders below anything a posterior resolves, and much smaller than the interpretation error it buys off; see anchored_frame.

Parallax is the one structural exception, and only because it composes as a reciprocal: ∂plx/∂plx_A = (plx/plx_A)², which differs from 1 by 2⋅dz⋅plx/2.06e8, with dz the anchor's line-of-sight offset from the barycentre — not the pair's separation, which is larger by the inverse mass fraction. That is 1e-8 for a planet-mass companion an AU out and 3e-5 for the 750 au stellar pair above. Add LL += 2*log(plx/plx_A) if you want the declared prior to be exactly a prior on the barycentric parallax; at these sizes it is a slowly-varying tilt across the prior's own support, not something a posterior can notice.

A physical prior on a derived barycentric quantity is one line, and needs no Jacobian either way — ~ on an expression is a likelihood term, not a change of sampling coordinates:

variables=AnchoredFrame(A; ref_epoch=…, variables=@variables begin
    …
    plx ~ Normal(74.05, 0.15)     # a prior on the *barycentre's* parallax
end)

Why not frame_shift=true

frame_shift reconditions the same degeneracy inside one G23HObs, by subtracting that observation's own Δpm from every channel. With two sources on one frame, each subtracts its own Δpm while redefining the same pmra, so both predict one number where the catalogue has two — for ups And, 13σ and 55σ apart. AnchoredFrame gets the conditioning without the redefinition, so it composes; G23HObs errors if frame_shift=true and more than one of them shares the frame.

See also anchor_offsets, System.

anchored_frame(system_interim, anchor, ref_epoch, correct; ra, dec, plx, pmra, pmdec, rv)

The barycentric frame quantities implied by an anchor source's catalogue solution, as (; ra, dec, plx, pmra, pmdec, rv) in System's units — the map AnchoredFrame emits, and the thing to call if you are writing the deferred lines yourself.

The six keywords are the inputs: the anchor source's observed values for whichever quantities are anchored, and the already-barycentric values for the rest. correct is a six-tuple of Bool in (ra, dec, plx, pmra, pmdec, rv) order saying which is which; an uncorrected quantity is returned unchanged. All six are required either way, because the second pass needs a complete frame.

Two passes, and why

The naive map subtracts anchor_offsets from the catalogue values. Those offsets come from the interim system, which is Parallax-level — it has a distance but no space motion — while the final system computes the same reflex through a full AbsoluteFrame. Two terms differ, and neither is negligible for a wide pair:

• proper motion. An angular offset ρ from the barycentre shrinks as the system recedes, at ρ⋅(ṙ/d). The interim's distance is constant, so it misses exactly that. Measured on a 750 au pair at 13.5 pc receding at 27 km/s: 4.0 × 10⁻³ mas/yr, which is 10% of a DR3 proper-motion σ.
• radial velocity. The interim's Einstein term is built from orbital velocities alone, missing ½|v_frame|²/c — 2.40 m/s for a 38 km/s space velocity, and a constant, since the frame's own speed does not vary with the orbit.

So the offsets are recomputed on a PlanetOrbits.reframed interim carrying the first pass's AbsoluteFrame — the same frame the final system will be built with, so the second pass computes the reflex on a system identical to the one the likelihood uses. What is left is the second-order feedback of the frame on its own offsets, ~10⁻⁴ of the first-pass residual: measured at 3 × 10⁻⁷ mas/yr and 2 × 10⁻⁴ m/s on the 750 au pair, i.e. 10⁻⁵ of a DR3 σ.

reframe rather than a second System(...) because nothing about the orbit may be recomputed here — see its docstring.

refine=false stops after the first pass. The reason to want that is the Jacobian: the one-pass map is exactly triangular with unit diagonal (the correction it subtracts is built from the body variables and the anchored parallax, and from nothing else), while the second pass makes each barycentric quantity depend weakly on all six inputs. Measured, that perturbs the diagonal by 10⁻⁹–10⁻⁴ and the determinant by 10⁻⁹ — far below anything a posterior can resolve, and much smaller than the 0.1σ interpretation error it removes, which is why refinement is the default. See AnchoredFrame's "Priors and the Jacobian".

anchor_offsets(system_interim, anchor, epoch)

The anchor body's motion about the system barycentre, at epoch [MJD], in the units the catalogue quantities are in. Subtracting these from an observed catalogue solution for the anchor source gives the corresponding barycentric frame quantities — see AnchoredFrame, which emits exactly that.

anchor is a body name (a Symbol, or a Body model node). Returns

The first four are angular and so need the interim system to carry a parallax: give the model a plx_interim variable (or a non-deferred plx), or they fail with a MethodError on NoFrame. rv and dz are frame-free.

This is the one-pass primitive. The anchored map itself is anchored_frame, which calls this twice — once at Parallax level and once on the frame that produced — because two terms of the reflex are invisible to a system with no space motion. See its docstring for both, with their measured sizes.

method, observing_geometry and barycentric_lighttime are PlanetOrbits' solve settings and default to the same values orbitsolve uses. Note the default propagator is Keplerian regardless of what the model is configured with: pass method= to match an N-body model, or accept that the frame is derived from the Keplerian approximation of the anchor's reflex (which is a correction to a correction, and below any catalogue's precision for the hierarchies where the two propagators visibly differ).

Two notes on exactness, both immaterial at present accuracies but recorded so nobody has to re-derive them:

• rv is radvel, not velz — a catalogue radial velocity is the spectroscopic quantity, so the Einstein term belongs in the subtraction. On a Parallax-level system that term is built from orbital velocities alone; anchored_frame's second pass is what supplies the rest.
• the interim is solved with no observer position, so ra_cosdec/dec are offsets as seen from the solar-system barycentre. That is what a catalogue solution reports, having already removed the annual parallax, so this is the right quantity rather than an approximation to one.

See also AnchoredFrame, PlanetOrbits.reframe.

barycentre_parallax(plx_anchor, dz_au)

Barycentre parallax [mas], given the anchor source's parallax [mas] and the anchor's line-of-sight displacement from the barycentre [AU, +receding] — i.e. anchor_offsets(...).dz.

Exact rather than a series: d_bary = d_anchor − dz in AU, inverted back to mas. The correction is about dz⋅plx²/2.06e8 mas — 1 part in 10⁵ for a body an AU from the barycentre at 100 pc, and 1 part in 300 for a 750 au wide pair at 13 pc, which is where it starts to matter.

logjac_cartesian_to_campbell(a, e, Mtot)

Log-Jacobian converting a flat prior in Cartesian relative phase space (x, y, z, vx, vy, vz at a reference epoch) into a flat prior in Campbell elements (a, e, cos i, M, ω, Ω). Add it as a potential term next to a Cartesian-parameterized orbit:

LL += logjac_cartesian_to_campbell(a, e, Mtot)

a is in AU and Mtot in solar masses — the row's total gravitating mass, the same one Orbit uses. Only the parameter-dependent part is returned; additive constants are dropped, which is exactly right for MCMC.

What flat-Cartesian sampling actually implies

Sampling Cartesian phase space with no potential term is not an uninformative prior on the elements. It is the Jeans/Ambartsumian phase-space prior, and it is a specific, opinionated one.

Delaunay variables are canonical, so phase-space volume is preserved:

d³x d³v = dL dG dH dl dg dh

with L = √(μ a), G = L√(1−e²), H = G cos i, and (l, g, h) = (M, ω, Ω) already the angles wanted. Transforming (L, G, H) → (a, e, cos i) is triangular in that order — L depends only on a, G on a and e, H on all three — so the determinant is the product of the diagonal:

∂L/∂a      = ½√(μ/a)
∂G/∂e      = −√(μ a)·e/√(1−e²)
∂H/∂cos i  = √(μ a)·√(1−e²)

|det| = ½ · μ^{3/2} · √a · e

Hence

p(a, e, cos i, M, ω, Ω) ∝ μ^{3/2} · e · √a

i.e. a thermal eccentricity distribution p(e) ∝ e, isotropic orientation (flat in cos i), p(a) ∝ √a, and a time-uniform phase.

That is the joint density, and the distinction matters in practice: what you sample is a bounded region of phase space (a box on x…vz, say), and its image in element space is not a product. A velocity bound cuts off small-r states, and it does so harder at large a — so the observed marginal in a is √a times the region's own a-dependent volume, and it is not √a. e and cos i are far less sensitive to that truncation, so they show the predicted shapes directly. Reweighting by this term fixes the density pointwise; it cannot undo the region's boundary.

The surprise is the mass tilt. μ = G·Mtot, so flat-Cartesian sampling also multiplies whatever prior you declared on the total mass by Mtot^{3/2}. That factor is invisible in the model text — nothing in a @variables block mentions it — and it is the one people are caught by, because the other three consequences are familiar and defensible while this one is a silent reweighting of a quantity the user thought they had set. If you sample Cartesian and want your declared mass prior to be the mass prior, you need this term.

Uniform-in-Campbell

Adding this potential cancels the whole tilt, leaving a flat prior in (a, e, cos i, M, ω, Ω) over the sampled support. Declare the priors you want on top of it as ordinary constraints on the derived elements:

sys = System(name="…", bodies=[A, B], observations=[…],
    variables=@variables begin
        plx ~ Uniform(1, 100)
        M_A ~ truncated(Normal(1.29, 0.016), lower=0)   # the primary's mass
    end)

A = Body(name="A", variables=@variables begin
    mass = system.M_A
end)

B = Body(name="B", about=A, variables=@variables begin
    mass ~ Uniform(0.05, 1.0)
    # sampled in Cartesian phase space ...
    x ~ Uniform(-2000, 2000); y ~ Uniform(-2000, 2000); z ~ Uniform(-2000, 2000)
    vx ~ Uniform(-5, 5);      vy ~ Uniform(-5, 5);      vz ~ Uniform(-5, 5)
    epoch = 57388.0
    # ... with the tilt removed, leaving flat (a, e, cos i)
    LL += logjac_cartesian_to_campbell(x, y, z, vx, vy, vz, mass + system.M_A)
end)

Two scoping points that bite here, both consequences of the evaluation order at the top of model/codegen.jl:

• The seven-argument call is required. With x…vz sampled, a and e are not in scope in the block — they are derived inside PlanetOrbits.Orbit, downstream of the parameter map. The three-argument form is for models whose elements are variables.
• The row's total mass has to come through a system variable (M_A above), not system.bodies.A.mass. A body block looks up and outward: it sees system.*, and siblings never see each other.

Edge cases

e → 0 and a → 0 are genuine poles, not numerical accidents: the flat-Cartesian density vanishes there, so reweighting it to uniform needs unbounded weight. The term returns -Inf-ward values approaching them, which a sampler reads as a rejection. If your model genuinely needs e = 0, a uniform-in-Campbell prior is not what you want at that point — sample the elements directly instead.

See also Orbit's Cartesian initial conditions, which is the parameterization this exists to serve.

logjac_cartesian_to_campbell(x, y, z, vx, vy, vz, Mtot)

The same potential term, written in the coordinates actually being sampled. State in AU and AU/julian-year, Mtot in solar masses.

This is the method to use with a Cartesian-parameterized body, and the reason it exists is a scoping fact that is easy to trip over: when x…vz are the sampled variables, a and e are not in scope in the model block. They are derived inside PlanetOrbits.Orbit, downstream of the parameter map, so the three-argument form has nothing to read there. This method extracts the two it needs — a from vis-viva, e from the eccentricity vector — from the state itself:

B = Body(name="B", about=A, variables=@variables begin
    mass = 0.2256
    x ~ Uniform(-2000, 2000); y ~ Uniform(-2000, 2000); z ~ Uniform(-2000, 2000)
    vx ~ Uniform(-5, 5);      vy ~ Uniform(-5, 5);      vz ~ Uniform(-5, 5)
    epoch = 57388.0
    LL += logjac_cartesian_to_campbell(x, y, z, vx, vy, vz, mass + system.M_A)
end)

Mtot has to reach the block through a system variable (system.M_A), not system.bodies.A.mass: a body block sees system.* but never a sibling.

Unbound and near-parabolic states have no a to speak of; the term runs to ±Inf there, which a sampler reads as a rejection. Bound the state box, or add an explicit e ~ Uniform(0, 0.95)-style constraint, if that matters.

RelAstromObs(data; target, ref, name, variables=@variables begin end)

Relative astrometry: the sky-plane offset of target from ref [mas].

astrom = RelAstromObs(tab; target=b, ref=A, name="GPI")

data needs an :epoch column [MJD] plus either :ra, :dec, :σ_ra, :σ_dec (optionally :cor) or :pa, :sep, :σ_pa, :σ_sep; angles are mas, position angles radians.

Both references take the full grammar — a body, Barycentre(A, b), Photocentre — so an inner binary's photocentre versus an outer companion is spelled directly rather than assembled from per-planet terms.

Variables

• jitter [mas] — added in quadrature to both components.
• platescale — multiplicative scale on the measured separation.
• northangle [rad] — rotation of the measured position angle.

RadialVelocityObs(data; target, ref=Barycentre, name, variables=…,
                  trend_function=…, gaussian_process=…)

Radial velocities [m/s] of target measured against ref, positive receding.

rvs = RadialVelocityObs(tab; target=A, ref=Barycentre, name="HARPS",
    variables=@variables begin
        offset ~ Normal(0, 100)      # m/s, instrument zero point
        jitter ~ LogUniform(0.1, 100) # m/s, added in quadrature
    end)

data needs :epoch [MJD], :rv [m/s] and :σ_rv [m/s].

One RadialVelocityObs per instrument: the offset and jitter are per instrument, and that is exactly the granularity this object carries.

Variables

• offset [m/s] — instrument zero point, added to the model. Optional, for both absolute and relative RV. v8's StarAbsoluteRVObs injected offset ~ Uniform(-1000, 1000) when no variables block was given; v9 never invents a prior, so a model that relied on that now has no offset unless it declares one.
• jitter [m/s] — added in quadrature to σ_rv. Optional, same reasoning.

Trends

trend_function = (θ_obs, epoch) -> θ_obs.trend_slope * (epoch - 57000)

Evaluated per epoch and added to the model alongside offset. Its parameters are ordinary variables of this observation.

Corrections and data provenance

Two keywords declare what your reduction already did, and one is an experimentation hook. See the "How Octofitter Computes Orbits" manual page.

secular_acceleration = :model | :data_corrected     # default :model

The perspective-acceleration drift — the barycentre's line-of-sight velocity changing as its 3D motion swings the line of sight. :model adds it to the prediction (its constant part washes into offset; the drift and curvature are the payload: ~4.5 m/s/yr for a Barnard-class star, ~3 cm/s/yr for a generic 30 km/s star at 30 pc). :data_corrected says your pipeline already removed it, using frozen catalog values, so the model adds nothing.

Why this one is model-side by default while the barycentric correction is not: the Earth-side terms depend only on the Earth ephemeris and the catalog direction — zero fit parameters — so a pipeline removes them once, perfectly. Perspective acceleration depends only on μ, ϖ and rv, which are sampled parameters here, and its signature (a smooth trend with slight curvature) is degenerate with a long-period companion. :model is also what makes a joint RV + absolute-astrometry fit self-consistent: the same sampled (μ, ϖ) predict both data types.

It requires an AbsoluteFrame to be non-zero; with plx alone the term is definitionally zero and no error is raised. It is rejected at construction for a relative series, with an explanation — a blanket "it cancels" is false, and the half that does not cancel is large.

secular_acceleration is this observation's only provenance dimension. Whether the Einstein term is modelled is a property of the whole fit, not of one instrument — no pipeline can have removed the part of it that varies with the orbit — so it lives on the system as System(…; einstein_rv=:on|:off).

Correlated noise

gaussian_process = θ_obs -> GP(θ_obs.gp_η₁^2 * SqExponentialKernel() ∘
                               ScaleTransform(1/θ_obs.gp_η₂))

Fits the residuals with a GP rather than assuming they are independent. The callable is handed this observation's variables and must return a GP object; the residuals it sees are at this observation's own epochs, in table order, which the constructor has sorted.

Two backends are supported and neither is a dependency of Octofitter: an AbstractGPs GP works through the duck-typed default path (gp_condition), and OctofitterRadialVelocity adds methods for its vendored Celerite. Load OctofitterRadialVelocity for the Celerite kernels.

Cross-validation

likeobj_from_epoch_subset(obs, rows) returns an observation whose ln_like is the log-likelihood of rows alone. Without a GP that is just the subset; with one the remaining rows are kept as the conditioning set and rows become held_out_table, since a correlated model cannot score a point without the points it is correlated with. Prediction is implemented for Celerite only — the AbstractGPs case throws, exactly as in v1.

MarginalizedRVObs(data; target, ref=Barycentre, name, variables=…,
                  trend_function=…)

Radial velocities [m/s] of target measured against ref, with this instrument's zero point analytically marginalized out rather than sampled.

rvs = MarginalizedRVObs(tab; target=A, ref=Barycentre, name="HARPS",
    variables=@variables begin
        jitter ~ LogUniform(0.1, 100)   # m/s, added in quadrature
    end)

data needs :epoch [MJD], :rv [m/s] and :σ_rv [m/s]. Declare no offset — that is the parameter this type integrates over. Everything else (target/ref, jitter, trend_function) means what it does on RadialVelocityObs.

Marginalizing removes one dimension per instrument, which usually samples better. Reach for RadialVelocityObs instead when you want a specific prior on the zero point, correlations between instruments' zero points, a hierarchical model over them — or a Gaussian process, which this type does not support.

Renamed from v8

MarginalizedStarAbsoluteRVObs — "Star" and "Absolute" were both just ref choices, and are now spelled as ones.

PhotometryObs(data; target, band=:default, name, variables=@variables begin end)

Photometry of one body in one band: compares target's flux_<band> variable against the measurements in data.

phot = PhotometryObs(tab; target=b, band=:H, name="NIRC2")

data needs :phot and :σ_phot columns, in whatever flux units the model's flux_H variables are in (setting the host's flux to 1.0 makes every other body's a contrast ratio). Any Tables.jl source works. One band per observation; for a second band, add a second PhotometryObs.

band=:H reads the variable flux_H from target's @variables block; band=:default (the default) reads a plain flux. This is the same variable-name → band mapping Photocentre weights are built from, so photometry and photocentres cannot disagree about what a body's brightness is.

When to reach for this, and when not to

v9's @variables blocks take ~ over derived quantities, so the comparison itself is a one-liner in target's own block:

flux_H ~ Normal(15.0, 3.0)      # an ad-hoc constraint, not data

That is the right spelling for an ad-hoc constraint. Use PhotometryObs when the numbers are data, because of the general criterion this refactor sorts every borderline case by:

If it carries data you might want to hold out or simulate, it is an observation. If it only reshapes the prior, it is a ~ line or an _isprior term.

A ~ line becomes a UserLikelihood with _isprior = true: it is excluded from likelihood counts, has no table to subset, and has no generate_from_params. So it can never be held out by cross-validation and you can never simulate photometry from a fitted model — both of which people routinely do with real photometry.

Variables

None are required. The forward model is the body's flux variable, which lives in the body's block rather than here — that is the flux/band unification: v8 declared flux inside the PhotometryObs's own variables block, so two instruments observing the same body in the same band each carried their own independent flux parameter.

GaiaDR4AstromObs(data; target=Photocentre, ref=Barycentre, name="GaiaDR4",
                 detrend=false, variables=…)

Gaia DR4 individual (epoch) astrometry: one along-scan abscissa per transit.

data needs :epoch [MJD], :scan_pos_angle ψ [degrees], :parallax_factor_al [dimensionless], :centroid_pos_al [mas] and :centroid_pos_error_al [mas]; an :outlier_flag column is honoured if present.

Units

scan_pos_angle is the only angle-like column, and it is taken in degrees — the unit the Gaia archive publishes it in (the DR4 VOTABLE declares unit="deg"), so a table read straight out of the archive needs no conversion. The conversion to radians happens once here, at construction: obs.table keeps the data verbatim in degrees, while the sinψ/cosψ the likelihood projects with are radian-based derived quantities. Everything else is unchanged: epochs in MJD, along-scan positions and uncertainties in mas, parallax factors dimensionless.

The model for each transit is

η = Δα* sin ψ + Δδ cos ψ + ϖ · f_al

where Δα*, Δδ are the barycentre's own linear motion (from this observation's ra_offset_mas, dec_offset_mas, pmra, pmdec, ref_epoch variables) plus the offset of target from ref — normally the photocentre relative to the system barycentre.

Parallax factors

f_al is Gaia's own parallax_factor_al, taken at face value: this observation consumes SSB observables plus explicit parallax factors, and never observer-aware observables. Every absolute-astrometry observation type does one or the other, never both, and which one is literal code in the type so that it can be reviewed per type.

What Gaia's factors omit is the annual–orbital (Kopeikin) coupling — the dependence of the parallax factor on the companion's line-of-sight depth, which is ≈ 4.85 · z[AU] / d[pc]² µas per AU of observer displacement. For any DR4 target that is sub-µas, well below the per-transit precision, so the face-value factors are exact enough. PlanetOrbits' observer-aware observables exist for the cases where it is not.

detrend

With detrend=true the linear (constant + slope) part of the target-versus-reference excursion is removed before it enters the model, so only the curvature does. That breaks the degeneracy between the fitted proper motion and a wide companion, and makes the observation's position and proper-motion variables describe the photocentre rather than the barycentre.

Variables

astrometric_jitter [mas] adds in quadrature to the per-transit formal error; ra_offset_mas, dec_offset_mas, pmra, pmdec, ref_epoch define the reference-point motion. The parallax comes from the system block.

Blended sources

target is what the catalog source is, and a catalog source is not generally a body: it is whatever flux the pipeline blended into one centroid. Photocentre (the default) is the whole system's flux-weighted point; Photocentre(:G, (Aa, Ab)) is the point over a named subset.

Two sources in a 2+2 quadruple — two tight pairs several arcseconds apart, so that only intra-pair blending is possible — are two instances of this observation, each with its own nuisance parameters, sharing the system's plx and frame:

System(name=:quad, bodies=(Aa, Ab, Ba, Bb, wide), observations=(
    GaiaDR4AstromObs(scans_A; target=Photocentre(:G, (Aa, Ab)),
                     ref=Barycentre, name="srcA", variables=@variables begin
                         ra_offset_mas ~ Normal(0, 100); dec_offset_mas ~ Normal(0, 100)
                         pmra ~ Normal(0, 100); pmdec ~ Normal(0, 100); ref_epoch = 57388.5
                     end),
    GaiaDR4AstromObs(scans_B; target=Photocentre(:G, (Ba, Bb)),
                     ref=Barycentre, name="srcB", variables=…),
), variables=…)

Each source's modelled signal then carries both its pair's wide-orbit motion and the intra-pair photocentric wobble, because a photocentre is one dot product over absolute body states — there is no per-level bookkeeping to get wrong. Bodies declare flux_G in their own blocks; setting the host's to 1.0 makes the others contrast ratios.

Membership that is not structurally fixed — a sampled resolved-flag, a scan-angle-dependent window — is not expressible as a static spec, and is not meant to be: an observation of that kind reads PlanetOrbits.fluxes(sys, band) and builds its own WeightedPoint per draw or per epoch.

G23HObs(; gaia_id, target, blends=(), …)

Joint Gaia DR2/DR3 + Hipparcos catalog astrometry from the G23H catalog: calibrated proper motions at five epochs, the Gaia DR3 astrometric excess-noise (UEVA/RUWE) channel, the Gaia RV variability channel, and optionally the Hipparcos per-transit abscissae.

Source membership

target=A, blends=(b, c, d)

target is the body the catalog source is centred on; blends are the other bodies whose light falls into the same source, in the order the flux-ratio variables are indexed. Both take Body model nodes or Symbols.

blends is photometry, not dynamics. A body left out of it still moves the target — every mass in the system does — it simply contributes no light to this source. So blends=() is the resolved source: body target alone, with its full orbital motion. It is not "no companions".

target/blends is not a single Photocentre spec, and cannot be: the Hipparcos branch is a grating response rather than a linear reduction (see _hippacentre!), and which bodies could ever blend is a structural, build-time statement while how much each one blends is a per-draw one.

Flux ratios

Each blend contributes light in proportion to its flux ratio against the target, in two bands: fluxratio (G, for the Gaia DR2/DR3 photocentre) and fluxratio_hip (Hp, for the Hipparcos abscissa).

Normally they come from the bodies themselves: a flux_G / flux_Hp variable in each body's block, read as f_k = flux_<band>(blend_k) / flux_<band>(target). Nothing is positional, and the same flux variable feeds every observation in that band. A model that declares no fluxes at all leaves every blend dark in every band — the right answer for a fit with no photometry — while a model that declares some band but not the one asked for is a name mismatch and errors.

Two ways to override that, both on the observation:

• fluxratio= / fluxratio_hip= on the constructor, for a constant. G23HObs(target=:B, blends=(:A,), fluxratio=0.0) says "A is resolved off B's source". The value is appended to the default variables block, so the defaults (σ_AL, transit_priorities, …) are kept.
• a fluxratio / fluxratio_hip entry in an explicit variables= block, for a ratio that is dynamic: derived from system variables, including deferred ones, so a resolved-flag latent at system level can gate a blend out of the photocentre for that draw. That is why blending state never has to round-trip through a body's own flux variable (which would be the body→deferred-system cycle codegen rejects), and it also covers a partial-blending weight that is not a physical flux ratio at all.

Either override must be a length-length(blends) container (a tuple, SVector, or vector) matching blends element for element; a bare scalar is accepted only when there is exactly one blend.

Data

catalog is the G23H catalog: a path to the Arrow file, or an already-loaded table (row selected by gaia_id), or a single row as a NamedTuple. Give either gaia_id or hip_id; hip_id is resolved to a Gaia source id through the catalog.

dr2_transits_catalog supplies the Gaia DR2 matched-transit count (astrometric_matched_observations_dr2), which sizes the DR2 epoch selection and is not carried by the published catalog. There is no fallback.

The Gaia scan geometry comes from one of:

• forecast_table — a prepared table with epoch [MJD], scanAngle_rad and parallaxFactorAlongScan. Nothing is fetched and no ephemeris is needed; this is also how the test suite runs offline.
• scanlaw_table — the scanninglaw Python package's output (times, angles); parallax factors are computed from the Earth ephemeris.
• neither — the GOST web service is queried for ra/dec.

hipparcos supplies the Hipparcos IAD as (; table, hip_sol) (see hipparcos_iad); by default it is loaded from the Hipparcos_IAD data dependency using the catalog's hip_id. A catalog row with hip_id = NaN drops every Hipparcos channel.

Channels

table.kind names the channels, one row each: :iad_hip, :ra_hip, :dec_hip, :ra_hg, :dec_hg, :ra_dr2, :dec_dr2, :ra_dr32, :dec_dr32, :ra_dr3, :dec_dr3, :ueva_dr3, :rv_dr3. likeobj_from_epoch_subset selects on those rows, so a caller can exclude a channel by dropping its row.

channels= restricts the set at construction, and does it by filtering the same table.kind rows — one code path, so the two spellings cannot diverge:

channels = (:ra_hip, :dec_hip, :ra_hg, :dec_hg, :ra_dr3, :dec_dr3)

is what HGCAObs is. Dropping the last Hipparcos channel drops the Hipparcos catalog distributions with it, and dropping :iad_hip drops the six abscissa nuisance parameters, exactly as a post-hoc subset would.

ueva_mode is :RUWE (default), :EAN, or :none. :none drops the :ueva_dr3 datum and the UEVA-driven DR3 covariance deflation — for stars whose sig_AL/sig_att_radec/sig_cal calibration is absent and which therefore cannot form the priors that channel needs — and leaves every other channel untouched.

frame_shift

Default true, and single-source behaviour is unchanged. Every channel — DR3's own included — is expressed relative to the model's DR3-epoch proper motion of the target body rather than of the barycentre, so the frame pmra stops meaning "the barycentre's" and starts meaning "the target's, as DR3 measured it". The Hipparcos–Gaia channel keeps its acceleration content because the shift is common-mode across epochs, and the reparameterization pays for itself: the frame proper motion is otherwise degenerate with the reflex, and this points the sampled direction at what the data pin hardest.

It is a reparameterization only under a wide prior on pmra/pmdec. Put an informative prior on either and frame_shift=true silently changes the posterior, because the prior is then a prior on a different quantity than it reads as. Nothing checks that for you.

With more than one G23HObs on the same system frame it is an error, not a warning. Each observation would subtract its own Δpm while redefining the one shared pmra, so every source's DR3 channel would predict the same number — for a wide pair like ups And, tens of σ from two catalogue values that genuinely differ. The relative proper motion is not merely uninformative under that shift; it is removed from the likelihood, and it is the entire wide-orbit signal. Use AnchoredFrame, which reconditions the same degeneracy in the model rather than inside one observation, and composes.

Variables

Defaults are derived from the catalog row unless variables is given: σ_AL, σ_att, σ_calib (mas, from the catalog's calibration), the epoch selections transit_priorities/transits/transits_dr2/transits_rv, the Hipparcos abscissa nuisance block (hip_iad_jitter, iad_Δra, iad_Δdec, iad_Δplx, iad_Δpmra, iad_Δpmdec), and σ_rv_per_transit. Optional σ_hip_pmra, σ_hg_pmra, σ_dr2_pmra, σ_dr32_pmra, σ_dr3_pmra (and the pmdec counterparts) add per-channel proper-motion jitter in quadrature.

The system's pmra/pmdec (or this observation's, if the system has no absolute frame) define the reference-point proper motion every modelled channel is expressed against.

Several sources in one system

Two Gaia sources belonging to one physical system are two G23HObs on one System, sharing one absolute frame. That shared frame is the whole point — it is what makes the wide pair's relative astrometry constrain the wide orbit — so give each observation its own target=/blends= and let the frame do the binding. Set frame_shift=false on all of them and parameterize the frame with AnchoredFrame.

Every nuisance parameter here stays per source, deliberately. transit_priorities, σ_AL, σ_att, σ_calib and u_dup_dr2 are declared in each observation's own namespace and are not shared, even for two sources a few arcseconds apart whose forecast pools are near-identical:

• the losses genuinely common to both (dead time, decontamination, the OBMT gap lists) are already removed by dr2_ok_mask/dr3_ok_mask before priorities are applied. What transit_priorities marginalizes is the residual, source-specific loss: gating, window class, saturation, per-source AGIS outlier rejection. For a pair 8.5 magnitudes apart these differ a lot, and not monotonically — the brighter star can have fewer usable transits, because bright-star handling discards them.
• sharing one priorities vector would force the smaller selection to be the top-k subset of the larger, which asserts a nesting the two selection mechanisms do not produce.
• σ_AL, σ_att and σ_calib are magnitude-dependent calibration terms and must not be shared under any scheme.

If you are tempted to share anyway, note the silent failure mode: transit_priorities indexes each observation's own forecast pool, built by a per-source GOST query. Two queries at different sky positions need not return the same number of rows, and if they differ by one, index i denotes a different transit in each — with no error, no warning, and a likelihood that is merely wrong.

HGCAObs(; gaia_id, target, blends=(), ref=Barycentre, kwargs...)

Proper-motion-anomaly astrometry in the style of the Hipparcos–Gaia Catalog of Accelerations: the Hipparcos and Gaia DR3 catalog proper motions and the long-baseline Hipparcos–Gaia scaled position difference.

This is a helper constructor, not a type. It builds a G23HObs restricted to the six channels the HGCA constrains:

HGCAObs(; gaia_id, target, blends=(), ref=Barycentre, kwargs...) =
    G23HObs(; gaia_id, target, blends, ref,
              channels    = (:ra_hip, :dec_hip, :ra_hg, :dec_hg, :ra_dr3, :dec_dr3),
              ueva_mode   = :none,
              include_iad = false,
              include_rv  = false,
              kwargs...)

so everything G23HObs documents — source membership through target/blends, flux ratios defaulting to the bodies' own flux_G / flux_Hp, variables=, the offline catalog=/forecast_table=/ hipparcos= inputs — applies unchanged, and likeobj_from_epoch_subset, generate_from_params and cross-validation come along for free.

Three parts of the mapping are load-bearing:

• ueva_mode = :none, not false. It is a three-valued symbol (:RUWE/:EAN/:none), and :none drops both the :ueva_dr3 datum and the UEVA-driven deflation of the DR3 covariance. The HGCA has no excess-noise channel, so both go.
• No DR2 channels. :ra_dr2/:dec_dr2 and the DR3−DR2 difference :ra_dr32/:dec_dr32 are dropped along with :rv_dr3: the HGCA is Hipparcos + Hipparcos–Gaia + DR3.
• Not bit-identical to v1. Those six channels are now modelled by G23H's treatment — its epoch-selection model, its five-parameter refits, its photocentre — rather than by the HGCA's own cross-calibration. The numbers move. This deserves one validation run against a published fit.

Example

A = Body(name=:A, variables=@variables begin
    mass ~ truncated(Normal(1.0, 0.1), lower=0.1)
    flux_G = 1.0                 # the target defines the contrast scale
end)
b = Body(name=:b, about=A, variables=@variables begin
    mass ~ Uniform(0, 100) * mjup2msol
    flux_G = 0.0                 # dark companion
    flux_Hp = 0.0
    a ~ LogUniform(1, 100); e ~ Uniform(0, 0.99)
    i ~ Sine(); ω ~ UniformCircular(); Ω ~ UniformCircular(); tp ~ Uniform(50000, 60000)
end)

sys = System(name=:HD1234, bodies=(A, b), observations=(
    HGCAObs(; gaia_id=756291174721509376, target=A, blends=(b,), ref=Barycentre),
), variables=@variables begin
    plx ~ truncated(Normal(24.0, 0.1), lower=0)
end)

See also G23HObs, which is what you want if you have the full G23H catalog row and can afford its extra channels.

HipparcosIADObs(; hip_id, target, blends=(), ref=Barycentre, …)

The Hipparcos per-transit abscissa residuals (van Leeuwen 2007 intermediate astrometric data) as a standalone likelihood.

G23HObs carries the same data as its :iad_hip channel, alongside the Gaia catalog channels; reach for this type when Hipparcos is all you have, or when you want the abscissae without the catalog proper motions.

Source membership

target=A, blends=(b, c)

target is the star the Hipparcos entry is centred on and blends are the other bodies whose light modulates its abscissa. ref is what the target's reflex is measured against — Barycentre by default. As in G23HObs, blends is photometry: a body left out of it still moves the target, so blends=() is a resolved source rather than a static one.

Hp-band flux ratios come from the same lookup G23HObs documents: the bodies' own flux_Hp variables as flux_Hp(blend) / flux_Hp(target), unless this observation declares its own fluxratio_hip, which wins. A model with no flux_Hp anywhere and no vector leaves every blend dark, which is the right answer when the other bodies contribute no light.

The forward model, and what it constrains

The measured abscissa carries the Hipparcos catalog's whole sky path, so the model has to reproduce it. That is the five-parameter frame offset of skypath.jl — FrameOffset — plus the source's own excursion:

b_i = (Δα + Δt·μα✱ + Δα_src)·cosϕ + (Δδ + Δt·μδ + Δδ_src)·sinϕ + ϖ·f_AL

with the residual |proj_meas_alongscan − b| compared against the renormalized per-transit σ, inflated by the BINARYS first-harmonic factor.

By default the parallax ϖ is the system's own plx with no offset (iad_Δplx = 0), while the position and proper motion are free nuisances with wide priors. So a Hipparcos-only fit constrains the parallax and the companion's reflex, and marginalizes the four linear terms — the standard setup, and the reason the frame block is offsets rather than fixed values. Override variables= to change that; iad_Δplx ~ Uniform(-10, 10) makes the parallax a nuisance too, exactly as G23HObs treats it.

Data

Give hip_id and the IAD is loaded through hipparcos_iad from the Hipparcos_IAD data dependency, or pass an already-loaded iad=(; table, hip_sol) (which is how the test suite runs offline). renormalize, attempt_correction and is_van_leeuwen are forwarded to the loader. recalibrate=true additionally applies the Brandt, Michalik & Brandt shift (+0.140 mas on the residuals, 2.25 mas of extra dispersion) that G23HObs applies unconditionally; it is off here so that this type reproduces the catalog data as published unless asked otherwise.

Variables

Defaults: hip_iad_jitter (excess per-transit dispersion, added in quadrature), iad_Δra, iad_Δdec, iad_Δpmra, iad_Δpmdec and iad_Δplx = 0, with iad_pmra/iad_pmdec derived by adding the offsets to the Hipparcos catalog proper motion.

ImageObs(data; targets, ref, band=nothing, name="images", variables=@variables begin end)

A set of reduced images, and the companions modelled in them.

ImageObs(image_dat; targets=(b, c), ref=A, band=:H, name="SPHERE")

data needs the columns (:image, :epoch, :platescale) — one row per image — where image is an AstroImage recentred so that index [0,0] is ref, epoch is MJD, and platescale is that image's plate scale in mas/pixel.

Sources

targets lists the bodies whose signal this likelihood models; ref is the point the images are centred on (usually the host). targets may be a single body or a tuple, and both take the usual grammar — a Body node, a Symbol, or a Photocentre for a pair the imager does not resolve.

targets is a structural statement about which sources the forward model includes, so it is deliberately not inferred from which bodies happen to have a flux in this band: "fit the image term with c alone" and "include b but let its flux go to zero" are different models, and only the first can be said by leaving b out.

Each target's brightness is its own flux_<band> variable, in the image's units (contrast, Jy, counts — whatever the pixels are, consistently). band selects among them and may be omitted when the bodies declare exactly one.

Contrast

Pass a contrast column (a callable mapping separation in pixels to the 1σ flux uncertainty there, e.g. contrast_interp) or a contrastmap column (an image of the same geometry). With neither, a contrast curve is measured from each image by contrast.

Variables

• platescale — multiplier on the plate scale of every image [default 1].
• northangle [rad] — rotation of the images relative to true north [default 0].

Neither is a flux: flux on an ImageObs is a v8 spelling and is rejected.

Example

image_dat = Table(;
    epoch      = [mjd("2016-01-01"), mjd("2017-01-01")],
    image      = [AstroImages.recenter(img1), AstroImages.recenter(img2)],
    platescale = [19.4, 19.4],
)

ImageObs(image_dat; targets=(b, c), ref=A, band=:H, name="SPHERE",
    variables=@variables begin
        platescale = 1.0    # or: ~ truncated(Normal(1, 0.01), lower=0)
        northangle = 0.0    # or: ~ Normal(0, deg2rad(1))
    end)

with flux_H declared on b and on c.

LogLikelihoodMapObs(data; target, ref, name="likemap", variables=@variables begin end)

One or more maps of log-likelihood versus Δright-ascension and Δdeclination, computed with some other tool, as a function of where one companion is.

LogLikelihoodMapObs(likemap_dat; target=b, ref=A, name="GRAVITY")

data needs the columns (:map, :epoch, :platescale) — one row per map — where map is an AstroImage recentred so that index [0,0] is ref, epoch is MJD, and platescale maps its pixels to mas.

target is the companion the map describes and ref the point it is measured from; both take the usual grammar. The map is interpolated at the modelled position of target relative to ref and the value added straight to the log-density, so whatever normalization, jitter and correlated noise the upstream tool applied is inherited verbatim.

An optional fillvalue column gives the value used outside the map and wherever it is not finite; the minimum finite value of each map is used by default.

Variables

• platescale — multiplier on the plate scale of every map [default 1].
• northangle [rad] — rotation of the maps relative to true north [default 0].

Example

likemap_dat = Table(;
    epoch      = [mjd("2016-01-01"), mjd("2017-01-01")],
    map        = [AstroImages.recenter(map1), AstroImages.recenter(map2)],
    platescale = [19.4, 19.4],
)

LogLikelihoodMapObs(likemap_dat; target=b, ref=A, name="GRAVITY")

InterferometryObs(data...; targets, ref=Barycentre, band, name, …)

Interferometric observables — closure phases, optionally squared visibilities, optionally projected onto a kernel-phase basis — of a collection of point sources.

vis = InterferometryObs(tab;
    targets = (A, b, c),
    ref     = A,
    band    = :K,
    name    = "GRAVITY",
    variables = @variables begin
        σ_cp_jitter ~ LogUniform(0.1, 100)
        platescale = 1.0
        northangle = 0.0
    end)

The forward model

V(u,v) = Σ_j f_j exp(−2πi (u Δα✱_j + v Δδ_j)) / Σ_j f_j

targets names the bodies in the sum, ref is the phase centre the offsets are measured from, and f_j is each body's flux_<band> variable — the host included, so a model reproducing v8's contrast-ratio convention gives the host flux_K = 1.0 and the companions their ratios.

Shifting the phase centre multiplies V by a global phase, which every supported observable (|V|², closure phases, kernel phases) is invariant to, so ref is a free choice; Barycentre is the default and A matches v8's spelling exactly.

targets is a structural declaration of which sources this observation's forward model contains, and is deliberately not defaulted to "every body with a flux in this band": "fit with c only" and "include b, but let its flux go to zero" are different models, and you must be able to say the first.

Data

Each row is one exposure. Give either a filename column pointing at an OI-FITS file (read by _prepare_input_row, with optional wavelength_min_meters / wavelength_max_meters cuts) or the prepared columns directly: epoch [MJD], u, v [inverse wavelengths, baseline × channel], cps_data, dcps [degrees, triangle × channel], index_cps1, index_cps2, index_cps3, and eff_wave [m]. Add vis2_data, dvis2 and a true use_vis2 to include squared visibilities. Rows are sorted by epoch.

Options

• kernel_phases=true selects KernelPhases over the default ClosurePhases; kp_correlation=false drops the spectral correlation within a kernel phase. In kernel-phase mode a jitter column and an optional kp_Cy column name, per exposure, the variables holding the kernel-phase jitter [degrees] and the spectral correlation coefficient; they may live in this observation's variables or the system's.
• fiber_coupling=true models single-mode fibre injection losses. fiber_pointing is where the fibre is pointed, as a reference — Photocentre(band) by default, or a body, or Photocentre(band, (A, b)) for an explicit subset — and each source's throughput is a function of its own offset from that point.

Variables

• σ_cp_jitter [degrees] — added in quadrature to the closure-phase uncertainties (closure-phase mode only; kernel-phase mode uses the per-exposure jitter variable named by the table).
• platescale, northangle [rad] — the instrument calibration, applied as in Octofitter.sky_offset.

GRAVITYWideKPObs(data...; targets, ref, band, name="GRAVITY-WIDE", …)

GRAVITY-WIDE preset: InterferometryObs with kernel_phases=true and fiber_coupling=true.

This is a preset, not the v8 type — v8's GRAVITYWideKPObs took its companion fluxes from a positionally-indexed flux vector on the observation and had no notion of which bodies it was modelling. targets and band are required, and the fluxes come from the bodies' flux_<band> variables.

ClosurePhases()

Model closure phases directly, with the per-observable uncertainties from the data and an optional σ_cp_jitter added in quadrature. Optionally also models squared visibilities, for rows whose use_vis2 is true.

Works for any array geometry.

KernelPhases(; correlated=true)

Project the closure-phase residuals onto the kernel-phase basis (the orthonormal basis of the closure design's row space built by the internal OctofitterInterferometry.kernel_phase_basis) and evaluate them under the block-structured correlation model of CKP. correlated=false keeps the projection but assumes independent spectral channels.

Requires GRAVITY's array geometry (6 baselines, 4 closure triangles).

ObsPriorONeil2019(wrapped; orbit=nothing, name=…)

The observable-based priors of K. O'Neil et al. 2019, "Improving Orbit Estimates for Incomplete Orbits with a New Approach to Priors: with Applications from Black Holes to Planets", applied on top of wrapped.

astrom = RelAstromObs(tab; target=b, ref=A, name="GPI")
System(…, observations=[ObsPriorONeil2019(astrom)], …)

Pass only the wrapper, not the wrapped likelihood as well: this object evaluates wrapped's own log-likelihood and adds the Jacobian term to it, so listing both double-counts the data.

The correction is only correct if you put Uniform priors on all Campbell orbital elements and a Uniform prior on period (not on semi-major axis). The period prior's range has a large effect on the fit and no recommendation was published with the paper.

The orbit

The Jacobian is a property of one orbit evaluated at the wrapped likelihood's epochs. In v8 that orbit was implicit — the likelihood was attached to a planet — and in v9 nothing attaches, so it is named:

• orbit=nothing (the default) takes the wrapped observation's own target, which is right for relative astrometry (target=b, ref=A → b's orbit) and for relative RVs.
• orbit=b names it explicitly. Stellar-reflex radial velocities need this: a RadialVelocityObs(…; target=A, ref=Barycentre) measures the host, whose motion is caused by the companion whose orbit the prior is about, and A has no orbit of its own.
• orbit=(b, c) sums the term over several orbits, which is what v8's system-attached method did over every planet.

The orbit meant by a body is the hierarchy row that places it, so any hierarchy convention may be named.

What changed inside

v8 recovered each epoch's solution by scanning the planet's solution vector for sol.t == epoch, falling back to re-solving the orbit when no match was found — O(N) per epoch, and a silent fallback for AbsoluteVisual orbits, whose stored solution times are light-travel corrected and so never compare equal. v9 reads the row's own elements and the trajectory's emission epoch, which is exactly what the propagator's Kepler solve used, at O(1).

Note the anomalies come from the row's Keplerian elements rather than from an integrated trajectory, deliberately: this is a prior over the sampled elements, so it means the same thing under KeplerianApprox and under a full N-body propagator, where the row only sets initial conditions.

OrbitOrderPrior(b, c, d, …; name="OrbitOrderPrior")

Assert that the named bodies' orbits are ordered in semi-major axis: a_b < a_c < a_d < …. Anything else gets -Inf.

System(…, observations=[astrom_b, astrom_c, OrbitOrderPrior(b, c)], …)

This breaks the labelling degeneracy of a multi-companion fit. Without it a sampler is free to swap which companion is called b and which c, and the posterior is a superposition of every permutation — technically correct, and useless to summarize.

Each argument is a Body model node or a Symbol naming one; the orbit meant is the row that places that body (see _placing_row), so a body may be named whatever hierarchy convention the model uses. v8's PlanetOrderPrior took Planets and looked their orbits up positionally with findfirst over keys(θ_system.planets).

Carries no data, so it cannot be held out or simulated: it is a prior term (_isprior = true), not an observation.

NonCrossingPrior(; bodies=(), name="NonCrossingPrior")

Forbid successive orbits from crossing: the zero-threshold case of LimitClosestApproachAUPrior, whose docstring documents both.

LimitClosestApproachAUPrior(soft_closest_approach_au; bodies=(), name=…)
LimitClosestApproachAUPrior(hard_closest_approach_au, soft_closest_approach_au; …)
NonCrossingPrior(; bodies=(), name=…)

Require successive orbits to stay apart: sorting the selected orbits by semi-major axis, the gap between one's apoapsis and the next one's periapsis must exceed hard_closest_approach_au [AU], and is softly penalised (a 1/x² repulsive term) below soft_closest_approach_au.

NonCrossingPrior() is the both-zero case — orbits may not cross — and is the usual thing to reach for.

bodies= restricts the prior to the rows placing those bodies; with no list, every hierarchy row in the system is included, which is what v8 did. Name the bodies in any system where the rows are not all planetary orbits of one star — comparing an inner pair's apsides with a wide binary's is not a meaningful test.

Carries no data: a prior term (_isprior = true), not an observation.

HillStabilityPrior(; bodies=(), name="HillStabilityPrior")

Require successive orbits to be separated by more than 2√3 mutual Hill radii — the Gladman (1993) criterion for a two-planet system on nearly circular, nearly coplanar orbits. Anything closer gets -Inf.

Sorting the selected orbits by semi-major axis, each adjacent pair (inner, outer) must satisfy

a_out − a_in  >  2√3 · R_H,     R_H = a_out · ((m_in + m_out) / 3M★)^(1/3)

where m_in/m_out are the total masses of the two rows' exterior bodies and M★ is the mass interior to the outer row, excluding the inner row's bodies.

bodies= restricts which bodies the constraint applies to.

octofit(
    [rng::Random.AbstractRNG],
    model::Octofitter.LogDensityModel,
    target_accept::Number = 0.8,
    ensemble::AbstractMCMC.AbstractMCMCEnsemble = MCMCSerial();
    adaptation = 1000,
    iterations = 1000,
    drop_warmup = true,
    max_depth = 12,
    initial_samples = pathfinder ? 500 : 250_000,  # deprecated
    initial_parameters = nothing,                  # deprecated
    step_size = nothing,
    verbosity = 2,
)

Sample from the posterior defined by model using Hamiltonian Monte Carlo with the No U-Turn Sampler from AdvancedHMC.jl.

octofit(model, 0.6, iterations=2000, adaptation=2000)   # correct
octofit(model, iterations=2000, target_accept=0.6)      # MethodError

For posteriors with widely separated modes, or with a discrete variable, reach for octofit_pigeons instead — NUTS cannot cross a low-density gap and cannot move a discrete parameter at all.

See also initialize!, octofit_rejection, octofit_pigeons.

octofit_rejection(
    [rng::Random.AbstractRNG],
    model::Octofitter.LogDensityModel;
    draws=100_000,
    verbosity=2,
)

Sample from the posterior defined by model using rejection sampling with the prior as the proposal distribution.

This sampler draws draws samples from the prior, evaluates the likelihood at each point, and accepts each sample with probability proportional to its likelihood. The accepted samples are independent (no autocorrelation), but the method can be very inefficient for high-dimensional problems or when the posterior is much narrower than the prior.

Returns an MCMCChains.Chains object, consistent with octofit and octofit_pigeons.

initialize!(model::LogDensityModel, fixed_params=nothing; kwargs...)

Find starting points for sampling, optionally guided by values you supply as a named tuple.

The values you give are a starting guess, not a constraint. They pin the parameters they name during the global optimization stage that looks for the posterior mode, and a guess covering every free variable skips that stage altogether and goes straight to the pathfinder variational approximation, using your values as its starting point. Nothing is held fixed after that: pathfinder varies every parameter, and sampling then proceeds normally over all of them. To fix a variable for real, give it a value in the model itself (a = 3.2 rather than a ~ … in a @variables block).

See startingpoints! to set the starting points outright, with no optimization or pathfinder at all.

The fixed_params guess can include:

• System-level variables: (; plx=24.4, pmra=10.2, ...)
• Body variables: (; bodies=(; b=(; a=1.5, e=0.1, ...), ...))
• Observation variables: (; observations=(; ObsName=(; var1=val1, var2=val2, ...), ...))

Available keyword arguments include:

• verbosity=1: control extra logging, can be 0 for silent, up to 4 for debugging info
• pathfinder_autodiff=AutoForwardDiff(): what autodiff backend to use for initialization (not necessarily the same one used for the model in general)
• nruns=8: how many runs of multi-pathfinder to use
• ntries=2: how many times can pathfinder fail and restart
• ndraws=1000: how many draws to return from the pathfinder approximation

Example:

init_chain = initialize!(model, (;
    plx=24.4,
    pmra=10.2,
    bodies=(;
        b=(;
            a=1.5,
            e=0.1,
        )
    ),
    observations=(;
        GaiaRV=(;
            offset_gaiarv=-50.0,
            jitter_gaiarv=0.1,
        ),
        GaiaDR4=(;
            astrometric_jitter=0.05,
        )
    )
))

startingpoints!(model::LogDensityModel, point::NamedTuple; ndraws=1000)
startingpoints!(model::LogDensityModel, point1::NamedTuple, point2::NamedTuple, ...)

Set the model's starting points to exactly the values you provide, instead of having initialize! choose them by global optimization and pathfinder.

Each point is a named tuple in the natural (constrained) parameter space, with the same shape as the optional second argument of initialize! — system variables at the top level, body variables under bodies, and observation variables under observations. Unlike initialize!, which takes a partial named tuple as a guess and searches for the rest, startingpoints! needs a value for every free variable in the model, since it is setting a complete starting point rather than seeding a search. Each point is transformed into the unconstrained space the sampler works in (via model.link) before being stored.

Given a single point — the common case — every starting point is set to that same value: model.starting_points = fill(mapped_point, ndraws). Given several points, one starting point is stored per named tuple, in the order given.

Returns an MCMCChains.Chains of the starting points, as initialize! does, so you can inspect what was set.

Example:

startingpoints!(model, (;
    M=1.05,
    plx=50.0,
    bodies=(;
        b=(;
            a=5.0,
            e=0.3,
            i=0.6,
            ω=1.2,
            Ω=2.4,
            θ=0.8,
        )
    ),
))

chain = octofit(model)

See also initialize!.

initialize_from_nss!(model; gaia_id, body=:b, catalog=:dr3, kwargs...)

Convenience function that queries the NSS catalog for gaia_id, converts the orbital solution to model parameters, and calls initialize! with those as a starting guess.

The NSS parameters anchor the global optimization search, but are not used as priors. All model parameters remain free during sampling.

Example

model = Octofitter.LogDensityModel(sys)
init_chain = initialize_from_nss!(model; gaia_id=4295745059252873600, body=b)

Any additional keyword arguments are forwarded to initialize!.

The method signature of Octofitter.hmc is as follows:

chain = advancedhmc(
    [rng::Random.AbstractRNG],
    model::Octofitter.LogDensityModel
    target_accept::Number=0.8,
    adaptation=1000,
    iterations=1000,
    drop_warmup=true,
    max_depth=12,
)

Sample from the posterior defined by model using Hamiltonian Monte Carlo with the No U-Turn Sampler from AdvancedHMC.jl.

using Pigeons
chain, pt = octofit_pigeons(model; n_rounds, n_chains=16, n_chains_variational=16)

Sample with parallel tempering (Pigeons.jl), for posteriors that are multimodal or that HMC explores badly — widely separated orbit families, detection-limit and completeness models, and anything with a discrete variable. Pigeons must be loaded by the user; the methods live in a package extension.

n_rounds is required: Pigeons runs 2^n_rounds scans, doubling each round.

using Octofitter, Pigeons
model = Octofitter.LogDensityModel(sys)
chain, pt = octofit_pigeons(model, n_rounds=10)

For expensive models on a multi-core machine, pass cores=N to run the sampler in N separate worker processes instead of threads:

chain, pt = octofit_pigeons(model, n_rounds=10, cores=8)

This is often about twice as fast for models that are slow to evaluate (many RV epochs, Gaia-Hipparcos absolute astrometry, images), at the cost of a minute or two of startup per run while the workers load packages and compile the model — so it is not worth it for small models, and a hint is printed when one path or the other looks clearly better. Results checkpoint each round under results/ in the current directory, and the same (chain, pt) is returned. Packages beyond Octofitter and its companion packages that are needed to reconstruct the model in a worker (rare) can be listed with dependencies=[SomePackage].

cores is a total budget. By default each worker process uses one core; for models with thousands of epochs, threads_per_process=2 (or 4) splits the budget into fewer workers that each also thread the trajectory solve — cores=16, threads_per_process=2 runs 8 workers × 2 threads. Chains within a worker sample one after another, so only reach for this once there is already one worker per chain: trading workers for threads on a chain-starved budget loses more parallelism than the threaded solve wins back.

pt is the Pigeons PT object, so Pigeons.stepping_stone(pt) gives the log evidence ratio against the reference (the prior-only model built by prior_only_model). For a log evidence, add the reference's own normalization — see prior_only_model's exclude_all=true note.

Also accepts a Pigeons.Inputs or an existing Pigeons.PT (to continue a run after Pigeons.increment_n_rounds!).

See also octofit.

Octofitter.CorrectionReport

The CorrectionDecisions for a built model, plus informational advisories that carry no decision — the size of the radvel Einstein term and of the secular-acceleration drift against each radial-velocity series' own uncertainty. Stored on the System and propagated into chain metadata.

Octofitter.CorrectionDecision

What :auto decided for one correction flag, and the evidence for it.

Fields: flag, resolved (the Bool actually used), source (:user when given explicitly, :auto when measured), impacts (a worst-first vector of ObsImpact), ndraws/nfailed, threshold (the accumulated-bias limit, in σ of the posterior mean), seed, and a human-readable note.

An impact of NaN means that observation type does not report predictions the test can compare, which is read as "keep the correction on".

Octofitter.ObsImpact

What one observation had to say about one correction flag.

impact is the per-point figure — the 99th percentile over draws of max_epochs |Δ prediction| / σ. bias is impact · √n, which is the quantity the decision is actually made on: a common-mode per-point shift of b·σ across n points moves a posterior mean by about b·√n·σ, so the same per-point impact matters far more to a 10⁵-transit dataset than to a 20-epoch one. NaN means the observation could not say.

recheck_corrections(model, chain; ndraws=300, seed=…, verbosity=1)

Re-run the build-time impact test on draws from the sampled posterior, and report whether the corrections the model was built with are still the right ones. Returns a CorrectionReport.

Two directions, and they mean different things:

• Escalation (warning). A flag resolved off under the priors would resolve on under the posterior — the posterior concentrated somewhere the correction matters, so the results may be biased and the fit should be re-run with the flag on. Expected to be rare.
• De-escalation (hint). A flag resolved on — usually because broad priors covered parameter space the posterior never visits — is comfortably unneeded under the posterior. Nothing is wrong; the next fit of the same data can set it :off and sample faster. This is the common case.

Draw from the target chain only. Tempered chains sit closer to the prior and would re-trigger exactly the broad-prior inclusions the posterior has ruled out.

construct_system(model, chain, i)
construct_system(model, θ)

The PlanetOrbits.System for posterior sample i — every body, every hierarchy row, and the frame, exactly as the likelihood saw them.

This replaces v8's construct_elements(chain, :b, i). There is no per-planet orbit object to hand back any more: a companion's motion is a property of the system, which is the whole point of the change. Query it the same way a likelihood does:

sys  = construct_system(model, chain, 1)
traj = orbitsolve(sys, epochs)
raoff.(traj, :b, :A)

drawfrompriors(system)
drawfrompriors(system; overrides, rng=Random.default_rng())

One draw from the model's priors, already expanded into the nested NamedTuple structure — system variables at the top level, then bodies, then observations.

overrides pins chosen free (~) variables to values you supply, in the same nested shape, and is the right way to build a parameter set for generate_from_params:

θ = drawfrompriors(sys; overrides=(;
    plx = 24.5,
    bodies = (; b = (; mass = 85mjup, a = 45.0, e = 0.15)),
))

The values are written into the flat parameter vector before it is expanded, so every derived (=) variable is recomputed from them. That is what mergeing into a template cannot do: merge replaces an entry in the already-expanded structure, leaving any derived variable that was computed from it stale — and generate_from_params reads the derived elements, not the sampled ones they came from. Overriding a derived variable is an error naming the model's free variables, rather than a silent no-op.

generate_from_params(obs::HipparcosIADObs, ctx; add_noise)

A new HipparcosIADObs whose abscissa residuals are the ones this sample predicts. res is shifted by the model residual and proj_meas_alongscan is rebuilt from it, so the regenerated observation is evaluated by exactly the same code path — the reconstructed reference sky path and the scan geometry are data and do not move.

generate_from_params(obs::G23HObs, ctx; add_noise)

A new G23HObs whose catalog values are those this sample predicts. Every covariance the likelihood will assemble at these parameters is mirrored here, including the DR2/DR3 cross-correlation and the UEVA-driven deflation — independent draws would leave the two conditional DR2/DR3 directions over-dispersed by 1/(1−ρ²) in whitened space, and a fit reads that excess as astrometric acceleration, i.e. spurious decades-period companions.

generate_from_params(obs::AbstractObs, ctx::ObsContext; add_noise)

Fallback: a prior-shaped term has no data to regenerate, so it passes through unchanged. An observation that does carry data and has not implemented forward simulation is an error rather than a silent pass-through — otherwise a "simulated" system would quietly retain the real measurements, and every downstream check (SBC rank statistics, completeness) would be measuring the wrong thing.

generate_from_params(system, θ=drawfrompriors(system); add_noise=false)

A copy of system whose observations hold data simulated at the parameters θ — the model's own generative direction, used for prior-predictive checks, posterior-predictive checks, simulation-based calibration and injection-recovery.

θ is a nested parameter NamedTuple, as produced by drawfrompriors, model.arr2nt(θ_flat), or mcmcchain2result.

With add_noise=true each observation adds a draw from its own noise model (including any jitter θ specifies); with add_noise=false the data are the noiseless model prediction, which is what you want to check that a fit recovers a known truth exactly.

The system is solved once over the union of every observation's epochs and each observation is handed the same ObsContext its likelihood would see, so simulated data and the likelihood that will be fit to them cannot drift apart.

generate_from_params(obs::ImageObs, ctx; add_noise)

A replicate image set drawn from the model — for posterior-predictive checks and simulation-based calibration.

The replicate is generated in the space this likelihood actually reads: each image is treated as the matched-filter output map the estimator above assumes, in which the expected value at a position is the flux of a source there and the uncertainty is that position's contrast. So each replicate is a fresh map holding the modelled sources, and nothing else — the original pixels are not reused.

Each source is injected as a point: the four pixels bracketing its modelled position are raised so that interpolating there returns exactly its flux. With add_noise=false the replicate therefore reproduces its own generating parameters exactly. With add_noise=true every pixel gets an independent draw of scale σₓ, the same contrast the likelihood reads.

generate_from_params(obs::LogLikelihoodMapObs, ctx; add_noise)

A replicate set of maps, each the measured surface translated so that its peak sits at the modelled position of target.

A log-likelihood surface is a reduction output, not a measurement with a noise model: there is no sampling distribution to draw a new one from. What can be done — and what this does — is preserve the surface's measured shape and put it where the model says the companion is, which is what a posterior-predictive overlay of a localization map wants. add_noise=true therefore has nothing to add and is an error rather than a silent no-op.

generate_from_params(obs::InterferometryObs, ctx; add_noise)

Replace this observation's closure phases (and squared visibilities, where present) with the model's, optionally perturbed by their own uncertainties.

v8 accepted add_noise and ignored it, which made simulation-based calibration on this type meaningless; the noise is drawn here.

prior_only_model(system; exclude_all=false)

A copy of system whose data likelihoods have been stripped out, so that sampling it samples the prior. Used for prior-predictive checks, for tempering, and — with exclude_all=true — for the reference log_Z0 that turns a Pigeons log evidence ratio into a log evidence.

Each real observation is replaced by a BlankLikelihood carrying the same name and the same @variables block, so the parameter vector keeps its shape and a chain from the prior-only model lines up column for column with one from the full model.

exclude_all=false (the default) keeps prior-shaped terms — lhs ~ dist and LL += lines, and the UnitLengthPrior behind each UniformCircular. They reshape the prior rather than adding data, so a "prior only" model is the one that still has them.

exclude_all=true drops those as well, leaving nothing but the declared prior distributions. That is what you want when the point is to normalize:

prior_model = Octofitter.LogDensityModel(Octofitter.prior_only_model(model.system, exclude_all=true))
_, pt_prior = octofit_pigeons(prior_model, n_rounds=10)
log_Z0 = stepping_stone(pt_prior)

pointwise_like(model, chain; verbosity=1) -> (likelihood_mat, epochs)

The log-likelihood of each individual data point under each posterior sample: an N_sample × N_data matrix, plus the epoch labelling each column.

Columns are in model order — observation by observation, then table row by table row. Prior-shaped terms (lhs ~ dist lines, LL += lines, the UnitLengthPrior behind a UniformCircular) are not columns: they carry no data, so including them would both invent data points and add the same constant to every column. Observations that carry data but have no epoch table contribute one column each, labelled NaN.

Feed the result to ParetoSmooth.jl for leave-one-out cross-validation:

likelihood_mat, epochs = Octofitter.pointwise_like(model, chain)

using ParetoSmooth
result = psis_loo(collect(likelihood_mat'), chain_index=ones(Int, size(chain, 1)))

Implementation

v8 built one whole System per data point and compiled a fresh RuntimeGeneratedFunction likelihood for each — N_data model builds before the first number came out. v9 needs one: the per-row observations are evaluated directly against a single trajectory solved over the union of their epochs, exactly as make_ln_like does it, so the cost is one solve per sample rather than one model compile per data point.

The consequence worth knowing about is that the columns sum to the model's total log-likelihood minus its prior-shaped terms. That is the quantity PSIS-LOO wants; v8's columns did not have that property.

calibrationhmc(system; θ=sample_priors(rng, system), target_accept=0.85, add_noise=true, kwargs...)

One simulation-based-calibration trial (Talts et al. 2018): draw parameters from the prior, simulate a data set at them, refit, and report where the true value falls in the resulting posterior.

Returns (priorsampledict, rdict, chains) — the drawn parameter values, the rank of each within the posterior as a percentage, and the chain itself. Repeated over many trials, a well-specified model and an unbiased sampler give a flat histogram of ranks for every parameter.

Keywords

• θ — the flat prior draw to calibrate against; defaults to a fresh one.
• add_noise — whether the simulated data get a draw from the noise model.
• target_accept, verbosity and anything else are forwarded to octofit.

sbctrial(system, chainparams, saveas)

Run one calibrationhmc trial and write it to disk as four files:

• <saveas>_sampler_parameters.toml — the sampler settings used
• <saveas>_parameters.toml — the parameter values drawn from the prior
• <saveas>_rank_stats.toml — the rank of each true value in the posterior
• <saveas>_chains.fits — the full chain

chainparams is a NamedTuple or Dict of keyword arguments; θ (a flat prior draw) and any sampler settings go in it. Run this script many times with different seeds, then histogram the rank statistics — see the "Simulation Based Calibration" tutorial.

completeness_jobs(; masses, separations, n_trials=5)

Generate a list of CompletenessJob descriptions for a completeness grid.

Each job specifies a (mass, separation) grid point and a trial index. Jobs are independent and can be dispatched to separate processes or cluster nodes.

Arguments

• masses — iterable of companion masses in M⊙ (v9 has one mass unit; write 5mjup for a Jupiter-mass grid point)
• separations — iterable of semi-major axes [AU]
• n_trials::Int=5 — number of independent trials per grid cell

Returns

Vector{CompletenessJob} — one job per (mass, separation, trial) combination.

Example: cluster dispatch

jobs = completeness_jobs(masses=10 .^ range(-1, 2, 15), separations=10 .^ range(-0.3, 1.7, 15), n_trials=10)
# Write job index from SLURM_ARRAY_TASK_ID
job = jobs[parse(Int, ENV["SLURM_ARRAY_TASK_ID"])]
result = run_completeness_trial(job, system, sampler; inject=my_inject)
# Save result...

run_completeness_trial(job, system, sampler; inject, add_noise=true, verbosity=0)

Execute a single injection-recovery trial: inject a companion, simulate observations, and sample the posterior. No detection decision is made here — that happens in assemble_completeness.

The returned CompletenessResult stores the full posterior chain and the true injected parameters, allowing detection criteria to be applied and iterated on after the fact.

Arguments

• job::CompletenessJob — job description (grid point + seed)
• system::System — template system with priors, observations, and bodies
• sampler — callable (model) -> chain; e.g. m -> octofit(m, iterations=5000)

Keyword Arguments

• inject — callable (mass, separation) -> NamedTuple; maps grid values to parameter overrides applied to the drawn prior sample. Must return overrides for free (prior) parameters only, not derived parameters. Example: (m, s) -> (; bodies=(; b=(; mass=m, a=s)))
• add_noise::Bool=true — whether to add measurement noise to simulated data
• verbosity::Int=0 — logging verbosity (0=silent, 1=info, 2=debug)

Returns

CompletenessResult containing the job, posterior chain, and true parameters.

Details

1. Draws parameters from system's priors using a seeded RNG
2. Overrides parameters using inject(job.mass, job.separation)
3. Simulates observations via generate_from_params
4. Builds a LogDensityModel from the simulated system
5. Initializes the sampler at the true parameters (see the note at the top of this file)
6. Calls sampler(model) to obtain a posterior chain
7. Returns the chain and true parameters for later analysis

Example

result = run_completeness_trial(job, system,
    model -> octofit(model, iterations=5000, verbosity=0);
    inject = (mass, sep) -> (; bodies=(; b=(; mass=mass, a=sep))),
)

# Inspect the posterior chain
result.chain
result.θ_true

assemble_completeness(results, detection_criterion; masses, separations)

Apply a detection criterion to a collection of CompletenessResults and assemble a CompletenessMap.

Detection is applied here — not during the trial — so you can call this function multiple times with different criteria to iterate on thresholds without re-running the sampler.

Arguments

• results — iterable of CompletenessResult
• detection_criterion — callable (chain, θ_true) -> Bool; returns whether the injected companion was recovered in a given trial
• masses — the mass grid used to generate the jobs
• separations — the separation grid used to generate the jobs

Returns

CompletenessMap with completeness fractions on the mass × separation grid.

Example

# Assemble with a Bayes factor threshold:
cmap_bf3 = assemble_completeness(results,
    (chain, θ) -> let p = mean(chain["b_planet_present"]); p/(1-p) > 3 end;
    masses=masses, separations=seps,
)

# Try a stricter threshold on the same results:
cmap_bf10 = assemble_completeness(results,
    (chain, θ) -> let p = mean(chain["b_planet_present"]); p/(1-p) > 10 end;
    masses=masses, separations=seps,
)

# Or a simple mass recovery criterion:
cmap_mass = assemble_completeness(results,
    (chain, θ) -> quantile(vec(chain["b_mass"]), 0.05) > 0.1;
    masses=masses, separations=seps,
)

completeness_map(system, sampler, detection_criterion; inject, masses, separations, n_trials=5, add_noise=true, verbosity=1) -> (CompletenessMap, Vector{CompletenessResult})

Compute a completeness map by running injection-recovery trials locally.

This is a convenience wrapper that calls completeness_jobs, run_completeness_trial, and assemble_completeness in sequence. Returns both the map and the full results vector, so you can re-apply different detection criteria without re-sampling.

For cluster-scale work, use the three-phase API directly.

Arguments

• system::System — template system
• sampler — callable (model) -> chain
• detection_criterion — callable (chain, θ_true) -> Bool

Keyword Arguments

• inject — callable (mass, separation) -> NamedTuple of parameter overrides
• masses — grid of companion masses [M⊙]
• separations — grid of semi-major axes [AU]
• n_trials::Int=5 — trials per grid cell
• add_noise::Bool=true — add measurement noise to simulated data
• verbosity::Int=1 — logging level

Returns

(cmap::CompletenessMap, results::Vector{CompletenessResult}) — the assembled map and the raw results for re-thresholding.

Example

using Octofitter, Distributions

# Run completeness map
cmap, results = completeness_map(
    sys,
    model -> octofit(model, iterations=5000, verbosity=0),
    (chain, θ) -> quantile(vec(chain["b_mass"]), 0.05) > 0.1;
    inject = (mass, sep) -> (; bodies=(; b=(; mass=mass, a=sep))),
    masses = 10 .^ range(-1, 2, length=12),
    separations = 10 .^ range(-0.3, 1.7, length=12),
    n_trials = 5,
)

# Plot
using CairoMakie
completenessplot(cmap)

# Try a different threshold without re-running:
cmap_strict = assemble_completeness(results,
    (chain, θ) -> quantile(vec(chain["b_mass"]), 0.05) > 1.0;
    masses = 10 .^ range(-1, 2, length=12),
    separations = 10 .^ range(-0.3, 1.7, length=12),
)

CompletenessJob

A lightweight, serializable description of a single injection-recovery trial. Contains the grid indices, physical values, and RNG seed — everything needed to reproduce the trial deterministically.

Fields

• i_mass::Int — index into the mass grid
• i_sep::Int — index into the separation grid
• i_trial::Int — trial number within this grid cell
• mass::Float64 — companion mass [M⊙]
• separation::Float64 — semi-major axis in AU (or period, depending on usage)
• seed::UInt64 — RNG seed for reproducibility

CompletenessResult

The output of a single injection-recovery trial. Stores the full posterior chain and the true injected parameters so that detection criteria can be applied (and re-applied) after the fact.

Fields

• job::CompletenessJob — the job description that produced this result
• chain::Chains — full posterior chain from the sampler
• θ_true::NamedTuple — the true injected parameter values

CompletenessMap

Assembled completeness results on a 2D grid of mass × separation.

Fields

• masses::Vector{Float64} — mass grid values in M⊙ (v9 has one mass unit; write 5mjup for a Jupiter-mass grid point)
• separations::Vector{Float64} — separation grid values [AU]
• completeness::Matrix{Float64} — fraction of trials detected (mass × sep)
• n_detected::Matrix{Int} — number of detections per cell
• n_total::Matrix{Int} — number of trials per cell

octoplot(model, chain; kwargs...) -> OctoPlotResult

One figure summarizing a fit: a sky panel of every orbit (when angular observables exist) and one time-series panel per data channel, with residuals and marginal histograms, all sharing a calendar-date epoch axis.

Requires a Makie backend to be loaded (e.g. using CairoMakie).

Returns a result whose fields are the figure, the named axes (for direct annotation: text!(res.axes.sky, ...)), and the underlying PosteriorSeries. Keywords are documented in the extension method; fname="..." saves the figure (nothing is written by default), and channels= restricts it to some of the data — an observable function, a channel or observable name, or a collection of either. A channels= the model has no data for is drawn as a prediction: the model curves alone, over the queries default_queries picks.

tmin=/tmax=/ts= set the epoch grid the curves are drawn over (see PosteriorSeries) — that, not xlims!, is how a figure is extended past the data — and curvecolor=/datastyle= override the model-curve colours and the data marks.

This is the many-draw figure, and three of its defaults follow from that: residuals are whitened and each point carries the boxplot of its z-score over the draws; phase-folded panels are off (a fold needs one ephemeris, and a posterior need not have one); and a fitted correlated-noise model rides in the curves rather than in a band (gpcurve), each draw carrying the Gaussian process conditioned on its own residuals. The first two turn back on for ndraws=1 or via show_phase=, and with one draw the noise model becomes the band gpband= draws instead. rvplot is the single-draw figure where all three are the point.

octocorner(model, chains...; small=false, includecols=[], excludecols=[],
           labels=Dict(), truth=(), fname=nothing, kwargs...)

Corner (pair) plot of the fit parameters. Labels, units, and radian→degree conversions come from PlanetOrbits.paraminfo — the same resolver table the axis labels use — keyed by the flat <owner>_<var> chain naming, so custom parameters simply show their column name. small=true keeps only each body's a, e, i, mass. UniformCircular helper pairs, fixed values, and tp duplicated by a sampled θ/M0 are dropped; includecols forces columns in, excludecols out. Extra keywords pass through to PairPlots.pairplot. Nothing is written unless fname is set.

Requires both a Makie backend and PairPlots to be loaded.

rvplot(model, chain, [sample_idx]; kwargs...) -> OctoPlotResult

The radial-velocity summary figure for a single posterior draw (sample_idx, by default the highest-posterior-density sample): one time-series panel carrying every instrument at once, with a residual strip and marginal histogram, plus one phase-folded panel per hierarchy row that moves the star.

This is the one figure allowed to put several RV instruments on one axis (see sharepanel), and showing a single draw is what buys that. A calibrated RV series is only defined per draw — the zero points, the jitters, the trend and the other rows' subtracted signals all move between samples — so everything here belongs to one sample and nothing is misrepresented. octoplot is the many-draws view and gives each instrument its own panel instead, with the data left uncalibrated and each draw's own offset, trend and activity model carried by its model curve.

The red binned points on a phase panel carry the uncertainty of the binned mean, max(1/√Σw, s_w/√n): the analytic weighted-mean error as a floor, and the bin's own bias-corrected weighted scatter over √n where the points disagree by more than the noise model allows. It is not the scatter of the points, which is what v8 drew. See phasebinmeans.

Requires a Makie backend. rvplot_animated records the same figure over successive single-draw slices of the chain.

Called rvpostplot before v9. The old name still works and forwards here: what the figure shows is one draw, not the posterior.

rvplot_animated(model, chain; N=50, seed=0, framerate=4,
                fname="rv-posterior.mp4", kwargs...)

rvplot recorded over N single-draw slices of the chain — v8's "sweep through the posterior" animation. Each frame is the whole figure rebuilt for one draw, so every panel (including the phase folds, which move with the drawn period) is that draw's own. Returns fname.

The extension is anything Makie can write from a VideoStream: .mp4, .gif, .mkv.

dotplot(model, chain; mode=:separation, epoch=nothing, kwargs...)

Mass against separation (or period) for every body in the fit, coloured by eccentricity, with marginal histograms. A posterior summary — no data, no observations — so it works for any model whose bodies have a sampled mass and e.

mode=:separation uses each draw's semi-major axis, or the instantaneous 3-D separation at epoch= when one is given; mode=:period uses the orbital period. Masses are M⊙ (v8 plotted Mⱼᵤₚ; v9 has one mass unit throughout). Requires a Makie backend.

gaiastarplot(model, chain, sample_idx=MAP; kwargs...)

The host's reflex orbit in the Gaia frame for one posterior draw, with each transit's along-scan residual re-projected into the sky plane along its own scan angle and drawn as a segment through the modelled track. This is the "is there a wobble, and does the orbit fit it?" picture for a GaiaDR4AstromObs; the along-scan-versus-time half of it is the generic panel octoplot already draws. Requires a Makie backend.

Draw several draws side by side with gaiastarplot!, which takes a grid cell instead of making its own figure.

gaiastarplot!(gridposition_or_axis, model, chain, sample_idx=MAP; kwargs...)

gaiastarplot into a cell of a figure you already have — gaiastarplot!(fig[i, j], model, chain, idx) — or into an axis you made yourself. Returns the Axis, so a grid of draws can be linked and have its interior decorations hidden in the usual Makie way. Requires a Makie backend.

gaiatimeplot(model, chain; kwargs...)

Along-scan abscissa against time for a GaiaDR4AstromObs: the posterior cloud of modelled abscissae over the measurements, with a per-epoch boxplot of the residuals against the quoted formal errors below.

This is the same data as octoplot's generic :along_scan panel, drawn in v8's per-epoch-boxplot idiom, which answers a different question: not "are the residuals normal" but "at which epochs is the posterior spread larger than the measurement error". Requires a Makie backend.

skytrackplot(model, chain, sample_idx=MAP; ra=nothing, dec=nothing,
             gaia_id=nothing, ts=nothing, keplerian_mult=1, kwargs...)

The star's whole path on the sky for one draw: parallactic loops, proper motion, and the orbital wobble superimposed — the picture of why the wobble is hard to extract. keplerian_mult exaggerates the orbital term.

The parallax ellipse needs a sky direction to project onto. It is taken from the system's own ra/dec frame variables when the model declares an absolute frame; otherwise give ra=/dec= in degrees, or gaia_id=, which reads them from the published solution via gaia_dr3_solution. Requires a Makie backend, and the DE440 ephemeris data dependency for the Earth's position.

hipparcosplot(model, chain, sample_idx=MAP; kwargs...)

Hipparcos intermediate astrometry for one draw, in its own geometry: the catalog's five-parameter sky path, the modelled path with the companion's perturbation, each transit's abscissa line, and the perpendicular residual and formal error drawn against it — plus a residual-versus-time strip.

Works with a HipparcosIADObs or with a G23HObs that keeps its :iad_hip channel. Requires a Makie backend.

completenessplot(cmap::CompletenessMap; kwargs...)

Heatmap of a completeness map over mass and separation. Requires a Makie backend to be loaded (e.g. using CairoMakie).

completenessplot!(gridposition, cmap::CompletenessMap; kwargs...)

Draw a completeness heatmap into an existing figure layout. Requires Makie. See completenessplot.

PlotChannel

One 1-D data channel of an observation, as declared by plotchannels:

• name — key into the NamedTuple residuals returns.
• label, unit — axis labelling.
• scale — display scale factor (e.g. rad2deg(1) for position angles); residuals already applies it, and a smooth model curve from query must be multiplied by it too.
• wrap — display period the quantity wraps at (in display units, e.g. 360 for position angle in degrees), or nothing.
• query — an ObservableQuery whose curve is this channel's smooth model prediction, or nothing when the channel has no meaning off the data epochs (e.g. Gaia along-scan abscissae, which depend on per-transit scan angles).
• derived — this channel is a deterministic re-expression of the observation's native measurement rather than a measurement in its own right: a sep/pa table's (ra, dec), say. The likelihood scores the native channels; a derived one exists so that a mixed dataset can put every point on every panel. octoplot only draws one when some other observation declares the same channel natively, or when channels= asks for it by name, and a consumer that wants the measurements themselves (a goodness-of-fit table) should filter these out.

plotchannels(obs) -> Tuple{Vararg{PlotChannel}}

The 1-D data channels this observation exposes for plotting; () if it is not plottable as data-vs-model series (prior terms return the default).

Implement this alongside ln_like and residuals when adding an observation type; a generic time-series panel then works with no plot-side code.

plotobs(obs) -> AbstractObs

The observation whose measurements a plot draws for obs — itself for every observation that carries its own data, and the wrapped observation for a wrapper like ObsPriorONeil2019.

A wrapper delegates the whole plotting protocol (plotchannels, residuals, sharepanel, datacalibration, noisemodel) to what it wraps, so the generic panels need no unwrapping. This exists for the few places that still ask what kind of observation is this — the sky panel's relative-astrometry overlay, the bespoke hipparcosplot/gaiastarplot selectors — where a isa test on the wrapper would silently drop the dataset. Test plotobs(obs) isa …, never obs isa ….

The wrapper itself, not the result of this, stays the object handed to obscontext and looked up in a series' data_maps: the fitted jitter/platescale/northangle are registered under the wrapper's likelihoodname, so unwrapping before building a context would silently fall back to the defaults.

defaultpanels(obs) -> Tuple

The escape hatch for observations that are not epoch series (an HGCA row, a catalog solution): bespoke panels instead of the generic time-series ones. Return () (the default) to use the generic panels derived from plotchannels; return (name => build, …) pairs to opt out — octoplot then calls each build(gridposition, series) in the panel stack and merges its returned NamedTuple of axes under name. Include a timeaxes key (a tuple of epoch axes) in that NamedTuple to have them linked with the figure's shared time axis.

A bespoke panel special-cases the drawing, not the plumbing: it still receives the shared PosteriorSeries (same draws as every other panel) and should source its residuals/whitening from residuals/ln_like machinery, never re-derive it.

sharepanel(obs) -> Bool

May this observation's data share a panel with another observation measuring the same quantity?

Only one draw's parameters can calibrate a panel, so several instruments on one axis means every instrument but that draw's is drawn slightly wrong. The question is whether "slightly" is small enough to be worth the far more readable figure, and the answer is a property of the observation type:

• false (the default). Radial velocity: an instrument zero point is an unconstrained free parameter of order the data range, so the calibrated series moves visibly from draw to draw. Each instrument gets its own panel, its data are drawn uncalibrated, and each draw's model curve carries that draw's own offset and trend (datacalibration) — which is the only way many draws and one dataset can appear together without misrepresenting either.
• true. Relative astrometry: platescale and northangle are calibration constants pinned to within a fraction of a percent, so every instrument's points land in the same place under any draw. Merging them onto one sky, separation and position-angle panel — calibrated by the maximum-posterior draw — is what makes the figure legible, and it is what Octofitter has always done.

rvplot is the deliberate exception on the other side: it puts every RV instrument on one panel because it shows a single draw, so there is no inconsistency to hide.

datacalibration(obs, ch::PlotChannel, ctx::ObsContext, epochs) -> Vector | nothing

The additive term that carries channel ch's pure-observable model curve into this observation's raw measured frame under the parameters in ctx — an RV instrument's zero point plus its trend, evaluated at epochs. nothing (the default) means the channel needs no additive calibration.

This is the inverse of what residuals applies to the data, and the two must agree: residuals subtracts it so the points lie on a pure model curve, and a panel drawing uncalibrated data adds it to the curve instead. Both spellings come from here, so they cannot drift apart.

noisemodel(obs, ctx::ObsContext, epochs) -> (; mean, var) | nothing

The observation's correlated-noise model — a Gaussian process fitted to the residuals — conditioned on this draw's residuals and evaluated at epochs; nothing (the default) when the observation has none.

Plots use it twice: what it predicts joins the model curve — as a band around one draw's, or, over many draws, added into each draw's own curve (noisecurves) — and the residual strip subtracts its mean and adds its variance to σ_eff. That second use is what makes a whitened residual meaningful for a GP fit at all: without it the strip shows exactly the correlated structure the GP was fitted to explain, and the z-scores are not standard normal even for a perfect fit.

A figure is not a fit, so this returns nothing — with a warning — rather than throwing where the noise model cannot be evaluated: a backend with no gp_predict method, or a draw whose covariance will not factorize (the same draw ln_like scores -Inf). The rest of the figure is still worth drawing, and it is then the Keplerian model alone. gp_predict itself still errors, so cross-validation — which has no partial answer to fall back on — fails loudly.

default_queries(sys, f) -> Vector{(ObservableQuery, Symbol)}

The natural queries of observable f for a system that has no observation declaring them — what octoplot draws when channels= names a quantity the model was not fitted to, so that a fit can predict a not-yet-observed signal.

Two conventions, by the kind of observable:

• separations (raoff, projectedseparation, posx, …): one query per hierarchy row, the exterior side about the interior side — exactly the relationship that row parametrizes, and the same set the sky panel draws.
• reflex signals (radvel, pmra, the velocities): one query per root body — the bodies no row places — against the whole-system barycentre.

The second element of each pair names the row (or body) the query belongs to, which is what the panel labels and colours itself by.

predictedchannels(sys, f) -> Vector{Tuple{Nothing,PlotChannel}}

Model-only channels for observable f — one per query default_queries picks, labelled from PlanetOrbits.plotinfo, with no observation behind them.

This is how a fit draws a quantity it has no data for: channels=radvel on a relative-astrometry fit predicts the reflex RV curve the orbit implies, which is the figure you want when deciding whether a target is worth spectroscopic time. Every observable is plottable this way whether or not it was observed; what makes a panel is the model, not the dataset.

ObservableQuery(f, target, ref)

The quantity f(sol, target, ref) as a value: an observable function (or its name as a Symbol) plus a target and a reference in the usual grammar (a Body node, a Symbol, Barycentre/Barycentre(A, b), Photocentre).

ObservableQuery(radvel, :A, Barycentre)     # stellar reflex RV
ObservableQuery(:radvel, :c, :b)            # planet–planet relative RV
ObservableQuery(raoff, Barycentre(:Aa, :Ab), Barycentre)

Anywhere a query is accepted, a plain tuple (f, target, ref) works too.

PosteriorSeries(model, chain; N=250, seed=0, ii=nothing,
                points_per_period=30, max_points=1000,
                tmin=nothing, tmax=nothing, ts=nothing)

Everything the plot layer needs from a fit, computed once and shared by every panel:

• a draw selection ii (default: N samples without replacement, seeded, so panels agree and reruns reproduce), plus the maximum-a-posteriori draw;
• the rebuilt PlanetOrbits.System per draw (construct_system);
• a dense epoch grid ts sized per hierarchy row (PlanetOrbits.plot_epochs), and each draw's trajectory solved over it;
• each draw's trajectory at the union of the data epochs, and ObsContexts for every observation — with method, observing_geometry and barycentric_lighttime forwarded from the model, so curves cannot silently disagree with the likelihood.

Accepts a chain from octofit or any MCMCChains.Chains with matching parameter names.

The epoch grid

By default the grid spans the data, padded 1.5 %, widened to at least the 35th-percentile orbital period over the draws, and clamped to 1900–2100. Since every panel draws its curves over exactly this grid — and clips its axis to it — it is also what sets how far past your data the figure extends.

• tmin= / tmax= replace one end of that span, keeping the automatic point density (points_per_period, max_points). Each is an MJD number, a date string (tmax="2035-01-01"), or a Date/DateTime; see mjd. tmin >= tmax is an error.
• ts= replaces the grid outright with the epochs you give (sorted, and duplicates dropped). The data epochs are not merged in as they are for the automatic grid, so a coarse ts= can step over fine structure — pass the epochs you want the curve evaluated at, and nothing else.

Give either ts= or tmin=/tmax=, not both. Data outside an explicitly requested window is still drawn as points (its residuals do not depend on this grid); only the model curves stop at the window's edge.

OctoPlotResult

What octoplot returns: the figure, the named axes (a nested NamedTuple — res.axes.sky.sky, res.axes.rv.main, res.axes.rv.resid, …) for direct annotation with ordinary Makie calls, and the underlying PosteriorSeries for further panels. Displays as its figure.

obscontext(series, obs; draw=nothing) -> ObsContext

The evaluation context for obs under posterior draw index draw (into series.ii), or under the MAP sample when draw === nothing. This is the same context the likelihood saw, so simulate, residuals and ln_like all evaluate consistently.

modelcurves(series, query) -> Vector{Vector{Float64}}

The query evaluated over the dense epoch grid series.ts, one vector per posterior draw. Display scaling (e.g. rad → deg for position angles) is the caller's job, via the channel's scale.

mapcurve(series, query) -> Vector{Float64}

The query evaluated over series.ts for the MAP draw.

noisecurves(series, obs, epochs; ndraws=length(series)) -> Vector | nothing

obs' correlated-noise model (noisemodel) at epochs, once per posterior draw: a vector of (; mean, var), the d-th conditioned on draw d's own residuals — its offsets, its trend, its jitter, its kernel hyperparameters — never on the MAP draw's.

This is what lets a many-draw figure draw the model each draw actually implies. A Gaussian process conditioned on one draw's residuals is part of that draw's model, so the honest ensemble is one orbit + activity curve per draw; the alternative is 250 Keplerian curves over residuals that visibly are not white.

Returns nothing — for the whole family, not just the offending draw — when the observation declares no noise model, or when any draw's cannot be predicted (noisemodel says why, once). A panel then falls back to the orbit alone, which is what it drew before there was an activity model to add.

timeseriespanel!(gridposition, series, channelgroup; kwargs...)

Generic data-vs-model time-series panel: posterior model curves, calibrated data with errorbars, a residual strip, and a marginal residual histogram. Requires Makie. See octoplot for the assembled default.

skypanel!(gridposition, series; kwargs...)

Sky-plane panel: phase-coloured orbit tracks for the default row queries (or explicit ones), overlaid relative-astrometry data, star marker. Requires Makie. See octoplot.

phasefoldpanel!(gridposition, series, entries; row, kwargs...)

Data-vs-model panel folded on hierarchy row row's orbital phase: the row's isolated signal (rowsignal) per posterior draw, calibrated data with the other rows' signals removed, noise-weighted binned means, and a phase-folded residual strip. Requires Makie. See octoplot.

photometrypanel!(gridposition, series, obs; kwargs...)

Bespoke panel for a PhotometryObs: the posterior of the modelled flux as a band, and the measurements with their errorbars. Photometry has no epoch axis, so the x axis is the measurement index. Requires Makie.

likemappanel!(gridposition, series, obs; kwargs...)

Bespoke panel for a LogLikelihoodMapObs: per epoch, how far below that epoch's map maximum the modelled position falls, over the posterior draws. Requires Makie (and OctofitterImages). See defaultpanels.

phasebinmeans(x, y, w, nbins) -> (; centre, mean, sigma)

Noise-weighted means of y in nbins equal bins of phase across [-0.5, 0.5), with w the weights a phase-folded panel pools over its instruments (1/σ_eff², so jitter and any fitted correlated-noise term are in them).

sigma is the error bar drawn on the binned mean: the larger of two estimates of the same quantity, the uncertainty of the mean itself.

σ_bin = max( 1/√(Σᵢ wᵢ),  s_w/√n )

The first term is the analytic error on a weighted mean of independent measurements with wᵢ = 1/σᵢ² — the RadVel/juliet convention, and exactly right in the white-noise limit. The second is the bin's own weighted sample standard deviation over √n: what the points in the bin actually do, rather than what the noise model says they should. The two agree when the noise model is honest and the residuals are white; the empirical term takes over when the residuals inside a bin are correlated (a fitted Gaussian process doing real work will do exactly this) or the quoted uncertainties are too small. Taking the maximum keeps the analytic value as a floor — a mean of n measurements cannot be more precise than that — while still letting real unmodelled structure inflate the bar.

s_w carries the frequency-weights bias correction,

s_w² = (Σᵢ wᵢ (yᵢ − μ)² / Σᵢ wᵢ) · n/(n−1)

which is what StatsBase's ProbabilityWeights correction (corrected=true) computes, and is identical to StatsBase's FrequencyWeights 1/(Σw − 1) correction applied to weights renormalised to sum to n. Renormalising is the whole point: w here is 1/σ², whose absolute scale is an inverse variance and not a sample size, so FrequencyWeights(w) reads Σw (≈ 0.19 for three 4 m/s points) as the number of measurements and returns a negative variance. Only the relative weights carry information; only n counts points. Before v9 this term used StatsBase's uncorrected default — the plain weighted rms, biased low, and drawn as though it were the error on the mean rather than the spread about it.

Bins holding fewer than two points are not returned at all. A mean of one point is that point: drawing it as a binned mean moves it from its own phase to the bin centre, and gives it the scatter of a single value — zero — in place of the error bar it already has on the axis. That is invisible when the bins are full and is the entire plot when they are not, so a fold with fewer points than bins yields no binned series rather than a red copy of its own data, one mark per measurement, each shouldered sideways and each claiming perfect precision.

Returns three same-length vectors, one entry per drawn bin.

residuals(obs, ctx::ObsContext) -> NamedTuple

Calibrated data, model predictions, residuals and uncertainties for each channel of plotchannels, at the observation's epochs, under the parameters in ctx — using exactly the likelihood's math (jitters in quadrature, platescale/northangle applied, instrument offsets removed, outlier masks respected).

Returns (; <channel name> = (; epoch, data, model, resid, σ, σ_eff, use), …) where data is calibrated into the model frame (so it overlays a pure-orbit model curve), resid is data - model in the same (display-scaled) units, σ is the measurement uncertainty and σ_eff includes fitted jitter. use is a Bool vector; false marks points the likelihood excluded.

Two optional keys carry per-channel extras when the observation has them: epoch_lo/epoch_hi bound an averaging window (a catalog proper motion is not measured at an epoch), and gp_mean/gp_var are this draw's correlated-noise prediction at the data epochs — see noisemodel. resid is always the plain data − model; a consumer that wants the residual the fit is left with subtracts gp_mean itself.

Not exported to avoid clashing with StatsBase.residuals; call it as Octofitter.residuals.

rvpostplot(model, chain, [sample_idx]; kwargs...)

Deprecated alias for rvplot, which is the same figure under a name that does not claim to show the posterior. See docs/src/v9-migration.md ("Migrating to Octofitter v9").

rvpostplot_animated(model, chain; kwargs...)

Deprecated alias for rvplot_animated. See docs/src/v9-migration.md ("Migrating to Octofitter v9").

savechain("saved-chain.fits", chain)

Save an MCMCChains.Chains to a FITS binary table, together with the run metadata in chain.info and the chain's section map.

Unicode parameter names are stored as their LaTeX spellings and restored by loadchain.

loadchain("saved-chain.fits")
loadchain("saved-chain.fits"; model)

Load an MCMCChains.Chains from a FITS binary table written by savechain.

Pass model to have the loaded chain checked against that model's parameter list (checkchain) before it is returned. That is the check worth doing: a chain and a model that disagree about a parameter's name do not error, they silently produce missing.

Files written by Octofitter v8 load, with a warning — their column names use v8's <planet>_<observation>_<variable> spelling for observations that hung off a companion, which no v9 model will match.

checkchain(model, chain; strict=true)

Assert that chain carries every free parameter model expects, and explain the mismatch if it does not.

Worth doing explicitly because the failure it catches is quiet: mcmcchain2result looks each parameter up by name and yields missing for anything absent, so a chain from a different model — or from Octofitter v8, whose planet-owned observations were named <planet>_<obs>_<var> — flows through plotting and post-prediction producing nonsense rather than an error.

With strict=false the mismatch is a warning instead. Returns chain.

savehdf5("filename.hdf5", model, chain)
savehdf5("filename.hdf5", model, chain, :b)

Save an Octofitter chain in the HDF5 layout used by orbitize! and by whereistheplanet.com.

Only the eight columns of orbitize!'s standard basis are written — sma, ecc, inc, aop, pan, tau, plx, mtot — so this exports one companion's visual orbit and nothing else. No data are exported.

The companion defaults to the model's first non-root body. mtot is the total mass of the bodies its orbit binds — <body>_mass plus the mass of everything the model places it about. For the usual about=A companion that is A_mass + b_mass; for a Jacobi chain (about=(A, b)) it is the sum over the whole interior, which is the same convention orbitize! uses for its own multi-planet fits.

loadhdf5("fname.h5")
loadhdf5("51erib")
loadhdf5(fname, numchains; host=:A, bodynames=("b","c","d","e"))

Load an orbitize! posterior from an HDF5 file and convert it into an Octofitter chain. Both tools use the same orbital-element conventions, so this is a rename plus one change of phase variable (orbitize!'s tau, a fraction of a period past a reference epoch, becomes Octofitter's tp).

Passing a name rather than a filename looks the target up on whereistheplanet.com.

numchains interprets the stored array as that many chains concatenated together.

host and bodynames give the v9 body names the columns are written under — m0 becomes <host>_mass and sma1, ecc1, … become <bodynames[1]>_a, <bodynames[1]>_e, …. Masses are solar masses, as orbitize! stores them and as v9 uses throughout.

Whereistheplanet_astrom("targetname"; target=:b, ref=:A)

Load relative astrometry hosted on whereistheplanet.com by a given target name, as a vector of RelAstromObs. If the name is not found, a list of similar target names is reported.

orbitize! stores separation/position-angle and RA/Dec rows in one table; Octofitter needs one likelihood object per format, so up to two are returned:

seppa, radec = Octofitter.Whereistheplanet_astrom("51erib"; target=b, ref=A)

target and ref are the model references the astrometry measures — the companion and the host in the usual case. They take the full v9 grammar (a Body, a Symbol, Barycentre(…), Photocentre(…)), because whereistheplanet's "object 1" is only a companion by convention.

The returned objects are named "<name>_seppa" and "<name>_radec"; both can go into the same System, and their names must differ for that to be legal.

sonora_photometry_interpolator(:Keck_L′, [metalicity="+0.0"])

Given a supported photometric band and [M/H] metalicity (default=solar), return a function of temperature (K) and mass that gives the absolute magnitude of the object in that bandpass.

absmag_L = sonora_photometry_interpolator(:Keck_L′)
absmag_L(1200.0, 12mjup)      # mass in M⊙, like every mass in v9

Out-of-grid inputs give NaN rather than an extrapolation.

Feeding a flux_<band> variable

The result is a magnitude. A body's flux_<band> variable is a linear flux — photocentres weight bodies by it — so convert:

flux_L = 10^(-0.4 * $absmag_L(tempK, mass))

or, if the host's flux is pinned to 1.0 so that companion fluxes are contrast ratios, take the magnitude difference to the host first.

Mass units

mass_unit selects how the second argument is interpreted: :Msol (default, matching v9 body mass variables), :Mjup (what v8 passed), or :Mearth.

Supported bands: :MKOY, :MKOZ, :MKOJ, :MKOH, :MKOK, :MKOL′, :MKOM′, :TwoMASSJ, :TwoMASSH, :TwoMASSKs, :KeckKs, :KeckL′, :KeckMs, :SDSSg′, :SDSSr′, :SDSSi′, :SDSSz′, :IRAC36, :IRAC45, :IRAC57, :IRAC79, :WISEW1, :WISEW2, :WISEW3, :WISE_W4

Supported metalicities: "+0.0", "-0.5", "+0.5"

itp = sonora_cooling_interpolator()

Create a function mapping (ageMyr, mass) -> tempK using Sonora Bobcat cooling model grids.

cooling = sonora_cooling_interpolator()
cooling(15.0, 12mjup)         # mass in M⊙, like every mass in v9

Out-of-grid inputs give NaN rather than an extrapolation. mass_unit selects the input unit: :Msol (default), :Mjup (the v8 convention), or :Mearth.

itp = bhac15_mass_age_interpolator(; key=:G)
itp = bhac15_mass_age_interpolator("BHAC15_iso.GAIA"; key=:G)

Create a function mapping (age_Myr, mass) -> absolute magnitude in the column named by key, using the BHAC15 model grids. With no filename the BHAC15_GAIA DataDep is used (and downloaded on first use).

itp = Octofitter.bhac15_mass_age_interpolator(key=:G)
itp(15.0, 0.08)               # mass in M⊙, like every mass in v9

key names a column of the isochrone file (:Teff, :G, :G_BP, … for the GAIA tables). Out-of-grid inputs give NaN rather than an extrapolation. mass_unit selects the input unit: :Msol (default), :Mjup (the v8 convention), or :Mearth.

Like the Sonora interpolators, this returns an absolute magnitude. A body's flux_<band> variable is a linear flux, so convert before assigning it:

flux_G = 10^(-0.4 * $itp(system.age, mass))

gaia_plx(; gaia_id)

A truncated Normal prior on parallax [mas] read from that source's Gaia DR3 astrometric solution, for use as plx ~ gaia_plx(gaia_id=…).

The truncation at ±10σ is deliberate and inherited: it keeps a sampler from walking into the negative-parallax tail, where the AU→mas conversion is undefined.

gaia_dr3_solution(; gaia_id) -> NamedTuple

That source's row of the Gaia DR3 gaia_source catalog, as a NamedTuple keyed by the catalog's own column names (parallax, pmra, pmdec, ra, dec, ref_epoch, phot_g_mean_mag, ruwe, …).

Queried from the ESA TAP service and cached at _gaia_dr3_final/source-<gaia_id>.csv in the working directory, so a repeated call — or a re-run on a cluster node with the directory copied across — is offline.

Useful for seeding an absolute-frame system block from the published solution:

cat = gaia_dr3_solution(gaia_id=756291174721509376)
sys = System(name="s", bodies=(A, b), observations=(pma,), variables=@variables begin
    plx   ~ truncated(Normal(cat.parallax, cat.parallax_error), lower=0)
    ra    = $(cat.ra)
    dec   = $(cat.dec)
    pmra  ~ Normal(cat.pmra, 10)
    pmdec ~ Normal(cat.pmdec, 10)
    rv    = 0.0
    ref_epoch = Octofitter.jd2mjd(2457388.5)
end)

(Note the asymmetry: $ interpolation is needed on = lines, whose right-hand sides are quoted and evaluated later inside the model, and is rejected on ~ lines, which already see the enclosing scope.)

query_nss(; gaia_id, catalog=:dr3)

Query the Gaia Non-Single Star (NSS) two-body orbit table for a given source ID. Returns a named tuple of the NSS solution columns, or nothing if no solution is found.

Responses are cached per source, the same way gaia_dr3_solution's are (see Octofitter._gaia_cache_path): a _gaia_nss_dr3/ directory in the working directory wins if it already exists, otherwise a per-user depot scratchspace.

Arguments

• gaia_id: Gaia source ID (integer)
• catalog: :dr3 or :dr4 (default :dr3)

nss_to_starting_point(nss_sol, model; body=:b)

Convert an NSS orbital solution (as returned by query_nss) into a named tuple suitable for passing to initialize! as a starting guess.

This function inspects the named body's free variables and maps NSS values to whichever parameterization you chose (Thiele-Innes or Campbell). Derived variables are skipped: they are computed from the free ones, so setting them would be silently discarded.

Mapped parameters

For Thiele-Innes bodies (free variables A, B, F, G):

• A, B, F, G from NSS a_thiele_innes, b_thiele_innes, etc.

For Campbell bodies (free a or P, and i, Ω, ω):

• Converts NSS Thiele-Innes constants to Campbell elements
• Sets a (semi-major axis in AU) or P (period in days)
• Sets i, Ω, ω — including the Ωx/Ωy pair a UniformCircular() uses

Common parameters mapped in both cases:

• e (eccentricity)
• tp (periastron time in MJD, if the body samples tp)

Arguments

• nss_sol: Named tuple from query_nss
• model: A LogDensityModel
• body: the Body to set, or its name (default :b)

Returns

A named tuple like (; bodies=(; b=(; e=0.3, A=5.2, ...))) ready for initialize!.

nss_model, nss_chain = nss_to_model_chain(nss_sol; plx=nothing, gaia_id=nothing, N=10_000)

Build a minimal model and chain from an NSS orbital solution, suitable for comparison plotting against your own posterior (e.g. via PairPlots or octoplot).

The returned nss_chain contains N draws from Normal distributions centred on the NSS best-fit values with the NSS-reported uncertainties. The returned nss_model is a two-body model — a host A and a companion b — whose companion is parameterized by the Thiele-Innes constants, exactly as the NSS table is. It therefore reproduces the NSS photocentre track, which is what the solution describes.

The total system mass is derived automatically from the NSS period, Thiele-Innes constants, and parallax via Kepler's third law, so you do not need to provide it. Parallax is taken from the NSS table if available, otherwise looked up from Gaia DR3.

Arguments

• nss_sol: Named tuple from query_nss
• plx: Parallax in mas. If nothing, uses the NSS table value or queries Gaia DR3.
• gaia_id: Gaia source ID (used to look up parallax if not in nss_sol)
• N: Number of draws (default 10_000)

Returns

(nss_model, nss_chain) — a LogDensityModel and MCMCChains.Chains object.

Example

nss_sol = query_nss(gaia_id=4295745059252873600)
nss_model, nss_chain = nss_to_model_chain(nss_sol)

# Compare with your posterior in a pair plot
using PairPlots
pairplot(
    "Posterior" => chain,
    "NSS"       => nss_chain,
)

# Or overlay on octoplot
octoplot(nss_model, nss_chain)

g23h_scan_uncertainty(; gaia_id=nothing, hip_id=nothing, catalog=nothing)
    -> (; gaia_id, phot_g_mean_mag, σ_AL, σ_att, σ_calib,
          σ_formal, n_ccd, σ_transit_formal, σ_transit_true)

The measured Gaia along-scan noise budget for one source, read from the G23H catalog (Thompson et al. 2026) — the same three numbers G23HObs builds its σ_AL/σ_att/σ_calib priors from, so a simulation based on these is consistent with what G23H assumes for that star.

Unlike a hand-picked "0.04 mas per scan", these are calibrated per source against real Gaia performance, and they vary by a factor of a few from star to star with magnitude, colour and scan geometry.

catalog is anything G23HObs accepts: nothing for the ~14 GB G23H_Catalog DataDep, a path to an Arrow file, a Tables.jl table, or a NamedTuple standing in for a single row.

The three terms, and where each applies

Gaia measures a source ~9 times per field-of-view transit (SM, AF1–AF9), and n_ccd is that number for this source: astrometric_n_good_obs_al_dr3 / astrometric_matched_transits_dr3. So

• σ_formal = √(σ_att² + σ_AL²) is the per-CCD formal uncertainty;
• σ_transit_formal = σ_formal / √n_ccd is the formal uncertainty of one transit-level abscissa;
• σ_transit_true = √(σ_transit_formal² + σ_calib²) is its actual scatter, including the calibration term Gaia's formal errors omit — for most sources σ_calib dominates it.

Consistency with G23HObs

G23HObs weights each transit of its GOST pool by σ_formal (fit_5param_prepared(…, σ_formal)) and then multiplies the resulting χ² by N_AL = astrometric_n_good_obs_al_dr3 / astrometric_matched_transits_dr3 before comparing it to the catalog's CCD-level statistics. The effective per-transit variance in that likelihood is therefore σ_formal² / n_ccd exactly, which is what σ_transit_formal is.

The same identification holds against the real data: for the three Gaia DR4 pre-release sources shipped with the documentation, σ_formal reproduces the median published per-CCD centroid_pos_error_al to 4–13% (Gaia-4 0.085 vs 0.082 mas, Gaia BH3 0.090 vs 0.083 mas, HD 114762 0.108 vs 0.124 mas).

See also gaia_dr4_transit_template, which turns these into a simulated DR4 table.

gaia_dr4_transit_template(; σ_al, gaia_id=nothing, ra=nothing, dec=nothing,
                          forecast_table=nothing, baseline=:dr4)

A transit-level table in Gaia DR4 epoch-astrometry format for any star, to hand to GaiaDR4AstromObs and fill in with generate_from_params.

One row per forecast field-of-view transit, with the columns that observation needs: epoch [MJD], scan_pos_angle ψ [degrees, the archive's own unit], parallax_factor_al, centroid_pos_al [mas] and centroid_pos_error_al [mas], plus outlier_flag. The measurements themselves are zero — this is a template, and generate_from_params(system, θ; add_noise=true) replaces them with the modelled abscissae plus a draw of centroid_pos_error_al.

Where the scan geometry comes from

Either GOST (queried and cached by GOST_forecast, the same service G23HObs uses) or a forecast_table you supply. A supplied table may be in either of the two forms already in circulation:

• GOST's own columns, i.e. what GOST_forecast returns; or
• G23HObs's forecast_table contract — epoch [MJD], scanAngle_rad, parallaxFactorAlongScan — so the same table can drive a G23HObs and this template, which is the cleanest way to keep a DR2/DR3+Hipparcos fit and a simulated DR4 fit on identical scan geometry.

ra/dec (degrees) are needed only to query GOST; pass gaia_id instead and the DR3 solution is looked up (and cached) for them.

The parallax factors are GOST's own, not recomputed from an Earth ephemeris. Against the Gaia-4 DR4 pre-release, GOST's forecast reproduces the published per-transit parallax_factor_al to 5e-5 rms and scan_pos_angle to 0.004° rms; computing the parallax factor from an Earth ephemeris instead is a 0.005 rms error, a hundred times larger, because Gaia observes from L2, ~0.01 AU from the geocentre. That is 0.07 mas at ϖ = 13 mas — around twice the per-transit precision.

σ_al

The per-transit along-scan uncertainty [mas], scalar or one value per transit. For a star with a G23H calibration, g23h_scan_uncertainty gives the measured value:

σ = g23h_scan_uncertainty(; gaia_id)
transits = gaia_dr4_transit_template(; gaia_id, σ_al=σ.σ_transit_true)

σ_transit_true (the actual scatter, calibration term included) is the right choice when you are both simulating and fitting these data, since it is then also the uncertainty the likelihood is told about. σ_transit_formal is what a real DR4 table would quote — smaller, with the difference left for the observation's astrometric_jitter to absorb.

What this does not model

GOST forecasts every scheduled transit; real DR4 loses some to dead time, and more to AGIS's own outlier rejection, which is harshest for bright stars. For Gaia-4 the chain is 122 forecast → 109 in the pre-release → 93 used by AGIS; HD 114762 (G = 7.2) keeps 63 of 89. Nor is any noise correlated between transits here. Drop rows and inflate σ_al if you want either.

forecast_table = GOST_forecast(ra_deg,dec_deg;baseline=:dr3)

Given an Ra and Dec position, retreive a forecast of Gaia observations from the GOST tool automatically. See tool URL here: https://gaia.esac.esa.int/gost/

Please be aware that others might be able to discover the target coordinates you searched for (though not who performed the search) via information leaked to the external service.

Baseline can be :dr3, :dr4, or :dr5.

mjd("2020-01-01")

The modified julian day of a calendar timestamp, for writing epochs down readably: mjd("2020-01-01") == 58849.0.

The timestamp is read on the TT (Terrestrial Time) scale, which is the scale epochs are wanted on — TT and TDB differ by periodic terms below 2 ms, which no orbit fit can see. The Conventions page of the manual has the full statement.

mjd(Date("2020-01-01"))

The modified julian day of a Date or DateTime, read on the TT scale. Same caveats as the string method above: no timescale conversion is performed.

mjd()

The modified julian day of right now, on the TT scale.

years2mjd()

Convert from decimal years (e.g. 1995.25) into modified julian date, rounded to closest second

mjd2date(modified_julian)

Get a Date value from a modfied julian day, rounded to closest day

Examples

julia> mjd2date(59160.8)
2020-11-08

ObsContext

What a likelihood is handed: the sample's parameters, its own variables, the PlanetOrbits.System built from those parameters, and the Trajectory solved once for the whole model at the union of every observation's epochs.

Use solutionat rather than indexing traj directly — the trajectory is over the deduplicated, sorted epoch union, not this observation's table order.

refspecs(obs)

Tuple of reference specs (see Barycentre) this observation resolves against the solved system. Used to validate names at model-build time and to document the observation in show.

epochs(obs)

The epochs [MJD] this observation needs the system solved at, in table order. Observations with no table return an empty vector.

solutionat(ctx, i)

Per-epoch solution for row i of this observation's table.

resolverefs(ctx, specs::Tuple)

Resolve a tuple of reference specs in one go, returning a tuple of BodyRef/WeightedPoint values in the same order. For observations that name several targets (ImageObs, InterferometryObs, G23HObs), call it once outside the epoch loop and index the result.

It accepts only a Tuple, deliberately. The obvious spelling —

refs = [ref(ctx, s) for s in obs.targets]     # don't

builds a Vector whose element type is the join of the resolved types, so a mixed (BodyRef, WeightedPoint) list widens to Vector{Any}; it allocates once per likelihood evaluation, and the widening propagates into every raoff call downstream. map over a tuple keeps each element's concrete type, constant-folds (the specs carry their content in type parameters), and allocates nothing. Restricting the signature to Tuple makes the fast path the easy path rather than a thing to remember.

likelihoodname(obs)

Name of an observation. It labels the observation's variables in the chain, so it must be unique within a system.

likeobj_from_epoch_subset(obs::RadialVelocityObs, rows)

An observation whose ln_like is the log-likelihood of rows of obs.

Without a GP the rows are independent, so that is literally the subset — identical to what RelAstromObs does. With a GP they are not: the score of a held-out point depends on the points it is correlated with, so the other rows are kept as the conditioning set and rows move to held_out_table. Both spellings answer the same question, which is what cross-validation asks.

rows = : means "all of them", i.e. the ordinary full-data likelihood, with nothing left to condition on.

sky_offset(sol, target, reference; platescale=1, northangle=0) -> (Δα✱, Δδ)

Single-epoch sky-plane offset [mas] of target from reference, carried onto the detector by the instrument's calibration. sol is a per-epoch solution (solutionat(ctx, i)) and the two references are already resolved — by ref(ctx, spec) or resolverefs, outside the epoch loop. Anything raoff accepts works, including a body name.

Convention

platescale and northangle mean exactly what they mean in RelAstromObs: they describe the instrument's calibration error, so that a reported measurement is corrected to the true sky by

sep_true = sep_reported × platescale
pa_true  = pa_reported + northangle        (position angle: N through E)

This function applies the inverse — it takes the model, which lives on the true sky, and produces what the instrument would have reported:

Δα✱ = ( Δα✱_model·cos(northangle) − Δδ_model·sin(northangle) ) / platescale
Δδ  = ( Δα✱_model·sin(northangle) + Δδ_model·cos(northangle) ) / platescale

which is a rotation of the position angle by −northangle and a shrink by platescale. So the values it returns are directly comparable with the numbers in the data table — the pixel grid of an image, the phase reference of a visibility — with no further correction.

That is the direction the four consumers want, and it is the opposite of what RelAstromObs does internally: RelAstromObs moves the data onto the true sky instead, because its residual and its σ live in the data's own (sep, pa) or (ra, dec) frame. The two are not interchangeable — a residual formed in the model frame is the data-frame residual rotated and scaled, which reweights it against an untransformed σ — so do not "unify" them without deciding which frame the σ belongs to.

sky_offset!(Δα✱, Δδ, ctx, target, reference; platescale=1, northangle=0)

Fill Δα✱ and Δδ with the sky-plane offset [mas] of target from reference at every epoch of this observation's table, in table order, with the instrument calibration applied. Returns (Δα✱, Δδ).

The buffers are overwritten, not accumulated into — unlike accumulate_offsets!, whose callers lay down a reference-point linear motion first. One entry is written per row of the observation's table; pass @alloc'd storage from ctx.buf rather than fresh vectors.

target and reference may be specs or already-resolved references; either way they are resolved once, outside the loop, and sincos(northangle) is computed once. Nothing is allocated and nothing is typed to Float64, so this is safe in the hot loop and under ForwardDiff.

See sky_offset for the sign and scale convention — it is stated there in full, because a sign error here propagates into four likelihoods.

sky_calibration(ctx) -> (platescale, northangle)

The two instrument-calibration nuisance parameters, read from this observation's own variables with their identity defaults (platescale = 1, northangle = 0) and promoted to the sample's number type so they stay ForwardDiff-clean when only one of the two is being sampled.

Every likelihood with a plate scale spells the lookup the same way; this is that spelling, so a new one cannot quietly pick a different default.

FrameOffset(Δra, Δdec, plx, pmra, pmdec)

The five-parameter astrometric solution of an instrument's own frame, in the tangent plane about that instrument's catalog position: a position offset [mas], a parallax [mas], and a proper motion [mas/yr].

Built by frame_offset from an observation's variables and consumed by frame_offset_alongscan. isbits, so building one per evaluation is free.

frame_offset(θ_obs, plx_anchor, ::Type{T}) -> FrameOffset{T}

Read an observation's frame-offset block from its own variables. The names are fixed, and are the ones G23HObs has used since the port:

Each is optional and defaults to zero, so an instrument that fits only some of the five declares only those. The proper motions are absolute because their natural anchor is the instrument catalog's own solution, which a @variables block expresses directly (iad_pmra = 4.53 + iad_Δpmra).

plx_anchor is supplied by the caller rather than looked up, because whose parallax it is, is a modelling decision: G23HObs anchors on the Hipparcos catalog value (its abscissa channel exists for the companion curvature, and the frame is pure nuisance), while a Hipparcos-only fit anchors on the system's own plx so that the data actually constrain it.

frame_offset_alongscan(off, Δt_yr, cosϕ, sinϕ, parallax_factor_al,
                       Δα=0, Δδ=0)

Project a FrameOffset onto one scan: the along-scan coordinate the instrument would have measured for a source at Δα/Δδ [mas] off the frame's own path, Δt_yr Julian years after the frame's reference epoch.

b = (Δra + Δt·pmra + Δα)·cosϕ + (Δdec + Δt·pmdec + Δδ)·sinϕ + ϖ·f_AL

Δα/Δδ carry whatever the source itself is doing — an orbital reflex, a photocentre wobble, the Hipparcos grating response. Compare b against the measured abscissa (proj_meas_alongscan in a Hipparcos IAD table).

gp_condition(gp, epochs, σ²) -> fx

Attach the observation's epochs and per-point white-noise variances to the user's GP, giving whatever object gp_ln_like and gp_predict consume.

The default is AbstractGPs' finite-GP call, gp(epochs, σ²). A backend that computes a factorization in place (Celerite) returns the GP itself.

gp_ln_like(fx, residuals) -> Real

Log-likelihood of residuals under the conditioned GP from gp_condition. Defaults to logpdf.

gp_predict(fx, residuals, epochs) -> (mean, var)

Posterior predictive mean and variance at epochs, given residuals at the conditioning epochs. Two callers: cross-validation scores held-out points against it, and noisemodel turns it into the correlated-noise band a plot draws.

OctofitterRadialVelocity implements both shipped backends — Celerite, and AbstractGPs via posterior/mean_and_var. (v8 had the Celerite case only, so cross-validating an AbstractGPs-correlated fit threw; that hole is closed.) There is no default, because there is no way to guess it: a backend that implements gp_condition and gp_ln_like must implement this too if its fits are to be cross-validated or plotted.