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
        jitter ~ LogUniform(0.1, 100.0)  # RV jitter (m/s)
    end
)

# Example with 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",
    gaussian_process = θ_obs -> GP(θ_obs.gp_η₁^2 * SqExponentialKernel() ∘ ScaleTransform(1/θ_obs.gp_η₂)),
    variables=@variables begin
        jitter ~ LogUniform(0.1, 100.0)     # RV jitter (m/s)
        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.

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

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)