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),
)
PlanetRelAstromLikelihood(data)

Represents a likelihood function of relative astometry 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.

StarAbsoluteRVLikelihood(
    (;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:
StarAbsoluteRVLikelihood(
    (;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.

PlanetRelativeRVLikelihood(
    (;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:
PlanetRelativeRVLikelihood(
    (;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),
)
PhotometryLikelihood(
    data,
    name="INSTRUMENT",
    variables=@variables begin
        flux ~ Uniform(0, 10)
    end
)

A likelihood 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 PhotometryLikelihood 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.

HGCALikelihood(;
    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.

plotchains(
    chain, planet_key;
    N=1500,
    ii = rand(1:size(chain,1)*size(chain,3), N),
    color=length(model.system.planets) == 0 || !haskey(chain, string(planet_key)*"_a") ? nothing : string(planet_key)*"_a",
    colorbartitle=color,
    clims=nothing,
    cmap=:plasma,
    alpha=30/length(ii),
    attime=nothing,
    kwargs...,
)

Draw samples from a posterior chain for a given planet given by name planet_key and visualize them in some way. Use kind to control what plot is made. A few options: :astrometry, :radvel, :trueanom, :meananom, :eccanom, :x, :y, :z, (:x, :y), :raoff, :decoff, :pmra, :pmdec, :accra, :accdec, :radvel, :posangle, :projectedseparation. See PlanetOrbits documentation for more details.

Inputs:

• chain The chain to draw from
• planet_key Planet name in the model (symbol)
• N=1500 Number of samples to draw for the plot
• kind=nothing Specify what kind of plot to make.
• ii=... Specific row numbers to use, if you want to e.g. plot the same 100 samples in a few different plots
• color="planetkeya" Column name to to map colors to. Semi-major axis by default but can be any column or an arbitrary array.
• colorbartitle=color Name for colourbar
• clims=nothing Tuple of colour limits (min and max)
• cmap=:plasma Colormap
• alpha=... Transparency of the lines

projectpositions(model, chains.planets[1], mjd("2020-02-02"))

Given the posterior for a particular planet in the model and a modified julian date(s), return ra and dec offsets in mas for each sampling in the posterior.

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)

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

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.