API Documentation

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

Planet([derived,] priors, [astrometry,], name=:symbol)

A planet (or substellar companion) part of a model. Must be constructed with a block of priors, and optionally additional derived parameters and/or sastrometry. name must be a symbol, e.g. :b.

System([derived,] priors, [images,] [propermotionanom,] planets..., name=:symbol)

Construct a model of a system. Must be constructed with a block of priors, and optionally additional derived parameters. You may provide ProperMotionAnomLikelihood() and/or Images() of the system. Finally, planet models are listed last. name must be a symbol e.g. :HD56441.

data = Table(
    (epoch = 5000, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 5050, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
    (epoch = 5100, ra = -505.7637580573554, dec = -66.92982418533026, σ_ra = 10, σ_dec = 10, cor=0),
)
PlanetRelAstromObs(data)

Represents relative astrometry observations between a host star and a secondary body. :epoch is a required column, in addition to either :ra, :dec, :σ_ra, :σ_dec or :pa, :sep, :σ_pa, :σ_sep. All units are in milliarcseconds or radians as appropriate.

In addition to the example above, any Tables.jl compatible source can be provided.

StarAbsoluteRVObs(
    (;epoch=5000.0,  rv=−6.54, σ_rv=1.30),
    (;epoch=5050.1,  rv=−3.33, σ_rv=1.09),
    (;epoch=5100.2,  rv=7.90,  σ_rv=.11);

    name="inst name",
    variables=@variables begin
        offset ~ Normal(0, 100)           # RV zero-point (m/s)
        jitter ~ LogUniform(0.1, 100.0)  # RV jitter (m/s)
    end
)

# Example with trend function and Gaussian Process:
StarAbsoluteRVObs(
    (;epoch=5000.0,  rv=−6.54, σ_rv=1.30),
    (;epoch=5050.1,  rv=−3.33, σ_rv=1.09),
    (;epoch=5100.2,  rv=7.90,  σ_rv=.11);

    name="inst name",
    trend_function = (θ_obs, epoch) -> θ_obs.trend_slope * (epoch - 57000),  # Linear trend
    gaussian_process = θ_obs -> GP(θ_obs.gp_η₁^2 * SqExponentialKernel() ∘ ScaleTransform(1/θ_obs.gp_η₂)),
    variables=@variables begin
        offset ~ Normal(0, 100)             # RV zero-point (m/s)
        jitter ~ LogUniform(0.1, 100.0)    # RV jitter (m/s)
        trend_slope ~ Normal(0, 1)          # Linear trend slope (m/s/day)
        gp_η₁ ~ LogUniform(1.0, 100.0)      # GP amplitude
        gp_η₂ ~ LogUniform(1.0, 100.0)      # GP length scale
    end
)

Represents a likelihood function of absolute radial velocity of a host star. :epoch (mjd), :rv (m/s), and :σ_rv (m/s) are all required.

In addition to the example above, any Tables.jl compatible source can be provided.

The offset and jitter variables should be defined in the variables block and represent the RV zero-point and additional uncertainty to be added in quadrature to the formal measurement errors.

When using a trend function, it should be a function that takes θ_obs (observation parameters) and epoch and returns an RV offset. Trend parameters should be defined in the variables block.

When using a Gaussian process, the gaussian_process parameter should be a function that takes θ_obs (observation parameters) and returns a GP kernel. GP hyperparameters should be defined in the variables block and accessed via θ_obs.parameter_name.

PlanetRelativeRVObs(
    (;epoch=5000.0,  rv=−6.54, σ_rv=1.30),
    (;epoch=5050.1,  rv=−3.33, σ_rv=1.09),
    (;epoch=5100.2,  rv=7.90,  σ_rv=.11);

    name="inst name",
    variables=@variables begin
        offset ~ Normal(0, 100)           # RV zero-point (m/s)
        jitter ~ LogUniform(0.1, 100.0)  # RV jitter (m/s)
    end
)

# Example with trend function and Gaussian Process:
PlanetRelativeRVObs(
    (;epoch=5000.0,  rv=−6.54, σ_rv=1.30),
    (;epoch=5050.1,  rv=−3.33, σ_rv=1.09),
    (;epoch=5100.2,  rv=7.90,  σ_rv=.11);

    name="inst name",
    trend_function = (θ_obs, epoch) -> θ_obs.trend_slope * (epoch - 57000),  # Linear trend
    gaussian_process = θ_obs -> GP(θ_obs.gp_η₁^2 * SqExponentialKernel() ∘ ScaleTransform(1/θ_obs.gp_η₂)),
    variables=@variables begin
        offset ~ Normal(0, 100)             # RV zero-point (m/s)
        jitter ~ LogUniform(0.1, 100.0)    # RV jitter (m/s)
        trend_slope ~ Normal(0, 1)          # Linear trend slope (m/s/day)
        gp_η₁ ~ LogUniform(1.0, 100.0)      # GP amplitude
        gp_η₂ ~ LogUniform(1.0, 100.0)      # GP length scale
    end
)

Represents a likelihood function of relative radial velocity between a host star and a secondary body. :epoch (mjd), :rv (m/s), and :σ_rv (m/s) are all required.

In addition to the example above, any Tables.jl compatible source can be provided.

The jitter variable should be defined in the variables block and represents additional uncertainty to be added in quadrature to the formal measurement errors.

An offset variable can optionally be included in the variables block to fit an RV zero-point. Unlike StarAbsoluteRVObs, the offset is not included by default.

When using a trend function, it should be a function that takes θ_obs (observation parameters) and epoch and returns an RV offset. Trend parameters should be defined in the variables block.

When using a Gaussian process, the gaussian_process parameter should be a function that takes θ_obs (observation parameters) and returns a GP kernel. GP hyperparameters should be defined in the variables block and accessed via θ_obs.parameter_name.

data = Table(
    (phot=15.0, σ_phot=3.0),
    (phot=14.8, σ_phot=0.5),
)
PhotometryObs(
    data,
    name="INSTRUMENT",
    variables=@variables begin
        flux ~ Uniform(0, 10)
    end
)

An observation type for comparing measured photometry points in a single filter band to data (provided here). Requires the :phot and :σ_phot columns. Can be provided with any Tables.jl compatible data source.

For multiple bands, create separate PhotometryObs objects.

The flux variable should be defined in the variables block rather than in the planet definition. This can be derived from physical models that take planet mass and other system parameters as input.

The name is used for variable naming in the chain output.

HGCAObs(;
    gaia_id=1234,
    variables=@variables begin
        fluxratio ~ [Uniform(0, 1), Uniform(0, 1)]  # array for each companion
    end
)

Model Hipparcos-Gaia Catalog of Accelerations (Brandt et al) data using a full model of the Gaia and Hipparcos measurement process and linear models.

The fluxratio variable should be an array containing the flux ratio of each companion in the same order as the planets in the system.

Upon first load, you will be prompted to accept the download of the eDR3 version of the HGCA catalog.

ObsPriorAstromONeil2019(astrometry_likelihood, period_prior)

Given a an astrometry likelihood (PlanetRelAstromLikelihood), apply the "observable based priors" of K. O'Neil 2019 "Improving Orbit Estimates for Incomplete Orbits with a New Approach to Priors: with Applications from Black Holes to Planets".

This prior correction is only correct if you supply Uniform priors on all Campbell orbital parameters and a Uniform prior on Period (not semi-major axis). This period prior has a significant impact in the fit and recommendations for its range were not published in the original paper.

Examples

astrom_like = PlanetRelAstromLikelihood(astrom_table)

# Apply observable based priors ontop of our uniform Campbell priors:
obs_prior = ObsPriorAstromONeil2019(astrom_like)

# The astrometry lieklihood object is passed as a first parameter
# since the obserable-based priors depend on the observation 
# epochs.

@planet b Visual{KepOrbit} begin
    # Instead of a prior on sma
    # a ~ Uniform(0.001, 10000)

    # Put a prior on period:
	P ~ Uniform(0.001, 2000) # yrs
    a = cbrt(system.M * b.P^2)

    # Keep sine prior on inclination
    i ~ Sine()

    # Rest are uniform
    e ~ Uniform(0.0, 1.0)
    ω ~ UniformCircular()
    Ω ~ UniformCircular()
    τ ~ UniformCircular(1.0)
end astrom_like obs_prior

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

mjd("2020-01-01")

Get the modfied julian day of a date, or in general a UTC timestamp.

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

Get the modfied julian day of a Date or DateTime object.

mjd()

Get the current modified julian day of right now.

years2mjd()

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

gaia_plx(gaia_id=12123)

Get a distribution (truncated Normal) of parallax distance in mas of a source with GAIA catalog id gaia_id.

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.

Visual{OrbitType}(..., plx=...)

This wraps another orbit to add the parallax distance field plx, thus allowing projected quantities to be calculated. It forwards everything else to the parent orbit.

For example, the KepOrbit type supports calculating x and y positions in AU. A Visual{KepOrbit} additionally supports calculating projected right ascension and declination offsets.

AbsoluteVisual{OrbitType}(..., ref_epoch=, ra=, dec=, plx=, rv=, pmra=, pmdec=)

This wraps another orbit object to add parallax, proper motion, and RV fields, at a given reference epoch.

Like a Visual{OrbitType} this allows for calculating projected quantities, eg. separation in milliarcseconds.

What this type additionally does is correct for the star's 3D motion through space (RV and proper motion) and differential light travel-time compared to a reference epoch when calculating various quantities. This becomes necessary when computing eg. RVs over a long time period.

ra : degrees dec : degrees parallax : mas pmra : mas/yr pmdec : mas/yr rv : m/s ref_epoch : years

TODO: account for viewing angle differences and differential light travel time between a planet and its host.

RadialVelocityOrbit(a, e, ω, tp, M)

Represents an orbit of a planet with only the information retrievable from radial velocity measurements. That is, without inclination, longitude of ascending node, or distance to the system.

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 (M_jup) that gives the absolute magnitude of the planet in that bandpass.

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, massMjup) -> temp_K using Sonora Bobcat cooling model grids.

octofit(
    [rng::Random.AbstractRNG],
    model::Octofitter.LogDensityModel
    target_accept::Number=0.8,
    ensemble::AbstractMCMC.AbstractMCMCEnsemble=MCMCSerial();
    adaptation,
    iterations,
    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.

using Pigeons
octofit_pigeons(model; nrounds, n_chains=16, n_chains_variational=16)

Use Pigeons.jl to sample from intractable posterior distributions. Pigeons must be loaded by the user.

using Pigeons
model = Octofitter.LogDensityModel(System, autodiff=:ForwardDiff, verbosity=4)
chain, pt = octofit_pigeons(model)

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...)

Initialize the model with optional fixed parameters provided as a named tuple. Fixed parameters will be held constant during optimization and sampling.

The fixed_params can include:

• System-level variables: (; plx=24.4, pmra=10.2, ...)
• Planet variables: (; planets=(; 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,
    planets=(;
        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, planet variables under planets, and observation variables under observations. Unlike initialize!, which accepts a partial named tuple and fits the rest, startingpoints! needs a value for every free variable in the model, since it is setting a complete starting point rather than pinning a few of them. 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,
    planets=(;
        b=(;
            a=5.0,
            e=0.3,
            i=0.6,
            ω=1.2,
            Ω=2.4,
            θ=0.8,
        )
    ),
))

chain = octofit(model)

See also initialize!.

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

Save an MCMCChain to a FITS file binary table.

Octofitter.loadchain("saved-chain.fits")

Load an MCMCChain from a FITS file binary table.

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.

Results are cached locally in _gaia_nss_{catalog}/ directories.

Arguments

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

nss_to_starting_point(nss_sol, model; planet_key=:b)

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

This function inspects the model's planet parameters and maps NSS values to whichever parameterization the user has chosen (Thiele-Innes or Campbell).

Mapped parameters

For Thiele-Innes models (planet has A, B, F, G):

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

For Campbell models (planet has 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, Ω, ω

Common parameters mapped in both cases:

• e (eccentricity)
• tp (periastron time in MJD, if the planet has tp)

Arguments

• nss_sol: Named tuple from query_nss
• model: A LogDensityModel
• planet_key: Symbol identifying which planet to set (default :b)

Returns

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

initialize_from_nss!(model; gaia_id, planet_key=: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 starting points.

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

Example

model = Octofitter.LogDensityModel(sys)
chain = initialize_from_nss!(model; gaia_id=4295745059252873600, planet_key=:b)

Any additional keyword arguments are forwarded to initialize!.

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

Build a minimal Octofitter model and chain from an NSS orbital solution, suitable for comparison plotting with 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 one-planet LogDensityModel using a ThieleInnesOrbit parameterization.

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)

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 [Mjup]
• 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 planets
• 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) -> (; planets=(; 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 module note)
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) -> (; planets=(; 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 [Mjup]
• 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) -> (; planets=(; 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 [Mjup]
• separation::Float64 — semi-major axis [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 [Mjup]
• 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