harv

harv is a Python package for inferring Keplerian orbital parameters of binary-star and star–exoplanet systems from time series data. Built on JAX, NumPyro, and unxt for units-aware computation.

It's pronounced harvey.

⚠️ Warning! ⚠️

harv is in rapid development and is pre-alpha. Meaning the API is not stable or guaranteed! Once the first version is released, we will have some guarantees about backwards compatibility, but there are no guarantees for API stability with the current development versions of this package. Sorry!

Installation

Requires Python 3.12+. Install from GitHub:

pip install git+https://github.com/adrn/harv

Or, if using uv:

uv add git+https://github.com/adrn/harv

Quickstart

TODO: move to a quickstart page in documentation and link here.

Radial velocity modeling

from unxt import Q
import harv
import harv.models as hm

# Load or create RV data with explicit units - some sample RV data:
data = harv.RVData(
    time=Q([56000.0, 56100.0, 56250.0, 56400.0, 56600.0], "day"),
    rv=Q([12.3, -8.7, 5.1, -14.2, 10.8], "km/s"),
    rv_err=Q([1.2, 0.9, 1.1, 0.8, 1.0], "km/s"),
)

# Set up a prior with default structure (log-uniform in period, etc.)
prior = hm.StandardRV().default_prior(
    period_min=Q(10, "day"),
    period_max=Q(1000, "day"),
    sigma_K0=Q(30, "km/s"),   # RV semi-amplitude scale
    sigma_v0=Q(10, "km/s"),   # systemic velocity prior width
)

# Run the rejection sampler
sampler = harv.RejectionSampler(prior, harv.RVModel())
samples = sampler.run(data, n_prior_samples=1_000_000, seed=42)

# Inspect results — quantities carry units:
print(f"Accepted {samples.n_samples} posterior samples")
print(f"Period: {samples['period']}")
print(f"Eccentricity: {samples['eccentricity']}")

Gaia epoch astrometry

import numpyro.distributions as dist
import quaxed.numpy as jnp

astro_data = harv.GaiaAstrometryData(
    time=Q([958.110978, 994.910525, 995.086642, 1010.091395, 1076.918577], "day"),
    al_position=Q([147.066, 379.996, 378.656, 74.666, -293.923], "mas"),
    al_position_err=Q([0.370, 0.446, 0.428, 0.270, 0.247], "mas"),
    scan_angle=Q([-59.047, -5.114, -5.783, -68.579, -134.155], "deg"),
    parallax_factor=jnp.array([0.70828, -0.46657, -0.45946, 0.19659, -0.21379])
)

prior = hm.StandardGaiaAstrometry().default_prior(
    period_min=Q(50, "day"),
    period_max=Q(3000, "day"),
    sigma_a0=Q(1, "au"),   # astrometric semi-major axis prior scale
    sigma_parallax=Q(1.0, "mas"),  # parallax prior width
    sigma_pos=Q(1, "mas"),  # position offset prior width
    sigma_vtan=Q(30, "km/s"),  # tangential velocity prior width
)

sampler = harv.RejectionSampler(prior, harv.GaiaAstrometryModel())
samples = sampler.run(astro_data, n_prior_samples=1_000_000, seed=42)

MCMC continuation

When the rejection sampler returns a small number of samples, you can refine with NumPyro MCMC, started from the posterior samples:

from harv.samplers import NumpyroSampler

samples = sampler.run(n_prior_samples=1_000_000, max_posterior_samples=128)
mcmc_sampler = NumpyroSampler(model=model, prior=prior)
mcmc_samples = mcmc_sampler.run(samples, seed=0)

How it works

harv is built on the same tricks as The Joker (Price-Whelan et al. 2017) to make rejection sampling practical for Keplerian orbit inference. It uses a two-level parameterization that separates orbital parameters into:

1. Nonlinear parameters (period, eccentricity, argument of pericenter, phase) — sampled directly from the prior via rejection sampling, and
2. Linear parameters (RV semi-amplitude, systemic velocity, astrometric offsets) — these can be analytically marginalized given each set of nonlinear parameters.

This approach makes rejection sampling possible for high-dimensional parameter spaces by reducing the effective dimensionality of the sampling problem. The main utility of the rejection sampler is to map out the multi-modal posterior distribution of orbital parameters in the low signal-to-noise or low number of observations regime, which is challenging or impossible for MCMC methods. The samples returned by the rejection sampler are exact draws from the posterior distribution, and can be used to initialize MCMC samplers for further refinement, if necessary.

All objects are Equinox modules (valid JAX pytrees), so they work seamlessly with jax.jit, jax.vmap, and jax.grad.

Documentation

TODO: Full documentation will be available at harv.readthedocs.io.

License

harv is free software released under the MIT License. See LICENSE for details.