diff --git a/.github/workflows/joss_pdf.yml b/.github/workflows/joss_pdf.yml new file mode 100644 index 0000000..b78decb --- /dev/null +++ b/.github/workflows/joss_pdf.yml @@ -0,0 +1,33 @@ +name: JOSS Draft PDF + +on: + push: + paths: + - joss_data/** + - .github/workflows/joss_pdf.yml + pull_request: + paths: + - joss_data/** + - .github/workflows/joss_pdf.yml + +jobs: + paper: + runs-on: ubuntu-latest + name: Build JOSS paper PDF + steps: + - name: Checkout repository + uses: actions/checkout@v4 + + - name: Build draft PDF + uses: openjournals/openjournals-draft-action@master + with: + journal: joss + # Path to the paper within the repo + paper-path: joss_data/paper.md + + - name: Upload compiled PDF artifact + uses: actions/upload-artifact@v4 + with: + name: joss-paper + # Output path where Pandoc writes the compiled PDF + path: joss_data/paper.pdf \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml index a1867e8..f8aca64 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -131,6 +131,11 @@ name = "gauss_prelim_orbit" harness = false required-features = ["jpl-download"] +[[bench]] +name = "iod_batch_diff" +harness = false +required-features = ["jpl-download"] + [profile.bench] opt-level = 3 lto = "thin" diff --git a/benches/gauss_prelim_orbit.rs b/benches/gauss_prelim_orbit.rs index 7621ca1..3c244d9 100644 --- a/benches/gauss_prelim_orbit.rs +++ b/benches/gauss_prelim_orbit.rs @@ -1,13 +1,34 @@ -//! Benchmarks for GaussObs::prelim_orbit (single-threaded) +//! Benchmarks for `GaussObs::prelim_orbit` (single-threaded). //! -//! Exemples d'exécution : -//! cargo bench --bench gauss_prelim_orbit --features jpl-download -//! cargo bench gauss_prelim_orbit -- gauss_prelim_orbit/single_call -//! cargo bench gauss_prelim_orbit -- gauss_prelim_orbit/batch_100 -//! cargo bench gauss_prelim_orbit -- gauss_prelim_orbit/noisy_batch_100 +//! # Overview //! -//! Astuce : vous pouvez aussi lancer en ligne de commande avec -//! RAYON_NUM_THREADS=1 cargo bench --bench gauss_prelim_orbit +//! This benchmark suite exercises the Gauss preliminary orbit determination +//! routine on a fixed synthetic fixture (`GaussObs`), in different usage +//! patterns: +//! +//! * Single call (micro-benchmark of the core algorithm). +//! * Fixed-size batches (e.g. 100 calls) to isolate algorithmic cost. +//! * Noisy realizations (Monte Carlo–like scenario). +//! * Scaling with batch size to study throughput and amortized overhead. +//! +//! All benchmarks are **single-threaded** by construction to avoid any +//! interference from Rayon or downstream parallelism. +//! +//! # Example commands +//! +//! ```bash +//! cargo bench --bench gauss_prelim_orbit --features jpl-download +//! cargo bench gauss_prelim_orbit -- gauss_prelim_orbit/single_call +//! cargo bench gauss_prelim_orbit -- gauss_prelim_orbit/batch_100 +//! cargo bench gauss_prelim_orbit -- gauss_prelim_orbit/noisy_single_call +//! cargo bench gauss_prelim_orbit -- gauss_prelim_orbit_batch_scaling/batch_size_256 +//! ``` +//! +//! Tip: you can also enforce single-threading from the command line with: +//! +//! ```bash +//! RAYON_NUM_THREADS=1 cargo bench --bench gauss_prelim_orbit +//! ``` #![cfg_attr(not(feature = "jpl-download"), allow(dead_code))] @@ -21,21 +42,64 @@ use outfit::outfit::Outfit; use outfit::outfit_errors::OutfitError; use rand::SeedableRng; -/// Build global Outfit state (ephemerides, frames, etc.) -/// Note: keep this outside the hot loops. +/// Build the global Outfit state (ephemerides, reference frames, etc.). +/// +/// This helper constructs an [`Outfit`] environment using the DE440 ephemeris +/// and the FCCT14 astrometric error model, matching the production setup. +/// +/// Arguments +/// ----------------- +/// * None. +/// +/// Return +/// ---------- +/// * A `Result` containing the initialized [`Outfit`] environment, or an +/// [`OutfitError`] if initialization fails. +/// +/// See also +/// ------------ +/// * [`make_fixture_gaussobs`] – Provides a deterministic Gauss fixture. +/// * [`bench_prelim_orbit`] – Uses this state for all prelim-orbit benchmarks. +/// * [`bench_prelim_orbit_vs_batch_size`] – Reuses the same state for +/// batch-scaling experiments. fn build_state() -> Result { - // English in-code comments per user preference: // Using FCCT14 as in production; adjust if you want to compare error models. use outfit::error_models::ErrorModel; Outfit::new("horizon:DE440", ErrorModel::FCCT14) } -/// Deterministic GaussObs fixture (angles in radians, times in MJD TT) +/// Create a deterministic `GaussObs` test fixture (angles in radians, times in MJD TT). +/// +/// The fixture encodes: +/// * Three observation indices (0, 1, 2), +/// * Right ascensions and declinations (radians), +/// * Observation times (MJD TT), +/// * Heliocentric positions of the observer at each epoch (AU). +/// +/// Arguments +/// ----------------- +/// * None. +/// +/// Return +/// ---------- +/// * A [`GaussObs`] instance with fixed values suitable for reproducible +/// benchmarking. +/// +/// See also +/// ------------ +/// * [`build_state`] – Returns the global Outfit environment. +/// * [`bench_prelim_orbit`] – Benchmarks `prelim_orbit` on this fixture. +/// * [`bench_prelim_orbit_vs_batch_size`] – Reuses the same fixture to explore +/// performance scaling with batch size. fn make_fixture_gaussobs() -> GaussObs { let idx_obs = Vector3::new(0, 1, 2); - let ra = Vector3::new(1.6894680985108945, 1.6898614520910629, 1.7526450904422723); + let ra = Vector3::new( + 1.689_468_098_510_894_5, + 1.689_861_452_091_062_9, + 1.752_645_090_442_272_3, + ); let dec = Vector3::new( - 1.0825984522657437, + 1.082_598_452_265_743_7, 0.943_679_018_934_623_1, 0.827_517_321_571_201_4, ); @@ -61,14 +125,32 @@ fn make_fixture_gaussobs() -> GaussObs { } /// Force Rayon (if used by downstream code) to a single worker thread. -/// Must be called before any Rayon pool is created by dependencies. +/// +/// This function sets the `RAYON_NUM_THREADS` environment variable to `"1"` +/// **before** any Rayon thread pool is created, ensuring that all parallel +/// code in dependencies is effectively single-threaded. +/// +/// Arguments +/// ----------------- +/// * None. +/// +/// Return +/// ---------- +/// * Nothing. The function is used for its side-effect on the process +/// environment. +/// +/// See also +/// ------------ +/// * [`bench_prelim_orbit`] – Calls this before any benchmark runs. +/// * [`bench_prelim_orbit_vs_batch_size`] – Also enforces single-threaded +/// execution for scaling benchmarks. fn force_single_thread_pool() { - // Pure env-var approach: works even if we don't depend on rayon directly. - // Must be set very early, before any rayon usage. + // Pure env-var approach: works even if we do not depend on Rayon directly. + // Must be set very early, before any Rayon usage. std::env::set_var("RAYON_NUM_THREADS", "1"); - // If you prefer hard enforcement and have rayon as a dev-dependency, - // you could uncomment this block (but it's optional): + // If you prefer hard enforcement and have Rayon as a dev-dependency, + // you could uncomment this block (but it is optional): // // #[cfg(any())] // { @@ -79,6 +161,34 @@ fn force_single_thread_pool() { // } } +/// Benchmark various usage patterns of `GaussObs::prelim_orbit` on a fixed fixture. +/// +/// This benchmark group (`gauss_prelim_orbit`) contains: +/// +/// 1. **single_call** – A single call to `prelim_orbit`, micro-benchmarking +/// the solver itself. +/// 2. **batch_100** – Reuse the same input and call `prelim_orbit` 100 times +/// per iteration, focusing on algorithmic cost. +/// 3. **noisy_single_call** – Generate noisy realizations of the same geometry +/// (Monte Carlo–style) and call `prelim_orbit` on each noisy instance. +/// +/// Single-threaded execution is enforced to avoid interference from any +/// Rayon-based parallelism in downstream code. +/// +/// Arguments +/// ----------------- +/// * `c` - Criterion handle used to register the benchmarks. +/// +/// Return +/// ---------- +/// * Nothing. The function registers benchmark cases in the given +/// [`Criterion`] context. +/// +/// See also +/// ------------ +/// * [`bench_prelim_orbit_vs_batch_size`] – Explores scaling vs. batch size. +/// * [`make_fixture_gaussobs`] – Provides the deterministic Gauss fixture. +/// * [`build_state`] – Constructs the shared Outfit environment. fn bench_prelim_orbit(c: &mut Criterion) { // Ensure single-threaded execution before anything else touches Rayon. force_single_thread_pool(); @@ -102,7 +212,7 @@ fn bench_prelim_orbit(c: &mut Criterion) { }) }); - // 2) Batch 100: reuse same input, measure algorithmic cost only + // 2) Batch 100: reuse same input, measure algorithmic cost only. group.bench_function("batch_100", |b| { b.iter_batched( || gauss.clone(), @@ -116,21 +226,20 @@ fn bench_prelim_orbit(c: &mut Criterion) { ) }); - // 3) Noisy: generate one noisy realization then run prelim_orbit + // 3) Noisy: generate noisy realizations then run prelim_orbit on each. group.bench_function("noisy_single_call", |b| { b.iter_batched( || gauss.clone(), |g| { - // 0.3" ≈ 1.454e-6 rad; use your production noise scale if different. + // 0.3" ≈ 1.454e-6 rad; we use ~1.5e-6 rad as a representative astrometric noise. let sigma_rad = 1.5e-6_f64; let mut rng = rand::rngs::StdRng::seed_from_u64(42); - let noisy = g.generate_noisy_realizations( - &(Vector3::zeros().add_scalar(1.0) * sigma_rad), - &(Vector3::zeros().add_scalar(1.0) * sigma_rad), - 100, - 1.0, - &mut rng, - ); + + // Per-angle standard deviations (same for RA and Dec here). + let sigma_vec = Vector3::zeros().add_scalar(sigma_rad); + + let noisy = + g.generate_noisy_realizations(&sigma_vec, &sigma_vec, 100, 1.0, &mut rng); for gg in noisy { let res = gg.prelim_orbit(&state, &IODParams::default()); @@ -144,5 +253,68 @@ fn bench_prelim_orbit(c: &mut Criterion) { group.finish(); } -criterion_group!(gauss_benches, bench_prelim_orbit); +/// Benchmark the scaling of `prelim_orbit` with the number of calls per batch. +/// +/// This group (`gauss_prelim_orbit_batch_scaling`) explores how the cost per +/// orbit-fit behaves as you increase the batch size: +/// +/// * For small batches, overheads (Criterion, function calls, etc.) dominate. +/// * For large batches, the raw cost of the preliminary orbit algorithm +/// becomes more apparent, yielding a stable "time per fit". +/// +/// This is particularly useful to derive a "time per orbit-fit" curve for +/// performance notes or publications. +/// +/// Arguments +/// ----------------- +/// * `c` - Criterion handle used to register the benchmarks. +/// +/// Return +/// ---------- +/// * Nothing. The function registers the batch-scaling benchmarks. +/// +/// See also +/// ------------ +/// * [`bench_prelim_orbit`] – Baseline single-call and fixed-batch benchmarks. +/// * [`make_fixture_gaussobs`] – Fixed geometry used for all fits. +/// * [`build_state`] – Shared Outfit environment reused across batches. +fn bench_prelim_orbit_vs_batch_size(c: &mut Criterion) { + // Ensure single-threaded execution before anything else touches Rayon. + force_single_thread_pool(); + + let mut group = c.benchmark_group("gauss_prelim_orbit_batch_scaling"); + + // Build global state once (outside hot loops). + let state = build_state().expect("Outfit state"); + + // Deterministic fixture. + let gauss = make_fixture_gaussobs(); + + // Different batch sizes to explore scaling behaviour. + let batch_sizes = [1_usize, 4, 16, 64, 128, 256, 512]; + + for &batch in &batch_sizes { + let name = format!("batch_size_{batch}"); + group.bench_function(&name, |b| { + b.iter_batched( + || gauss.clone(), + |g| { + for _ in 0..batch { + let res = g.prelim_orbit(&state, &IODParams::default()); + black_box(&res); + } + }, + BatchSize::SmallInput, + ) + }); + } + + group.finish(); +} + +criterion_group!( + gauss_benches, + bench_prelim_orbit, + bench_prelim_orbit_vs_batch_size +); criterion_main!(gauss_benches); diff --git a/benches/iod_batch_diff.rs b/benches/iod_batch_diff.rs new file mode 100644 index 0000000..8646be1 --- /dev/null +++ b/benches/iod_batch_diff.rs @@ -0,0 +1,284 @@ +//! Benchmark Gauss Initial Orbit Determination (IOD) on survey-like +//! trajectories loaded from a Parquet file. +//! +//! # Overview +//! +//! This benchmark loads a Parquet dataset, builds a `TrajectorySet`, and +//! evaluates the performance of Gauss IOD across +//! increasing batch sizes of *distinct* trajectories (survey-like scenario). +//! +//! The goal is to estimate the average cost per orbit-fit when processing +//! large collections of independent tracks, similar to real survey pipelines. +//! +//! # Usage +//! ```bash +//! cargo bench --bench iod_batch_diff --features jpl-download +//! ``` +//! +//! The `jpl-download` feature is required so that the DE440 ephemeris can be +//! fetched if not already cached locally. + +#![cfg_attr(not(feature = "jpl-download"), allow(dead_code))] + +use std::cell::RefCell; + +use camino::Utf8Path; +use criterion::{black_box, criterion_group, criterion_main, BatchSize, Criterion}; +use outfit::constants::ObjectNumber; +use outfit::error_models::ErrorModel; +use outfit::initial_orbit_determination::gauss_result::GaussResult; +use outfit::initial_orbit_determination::IODParams; +use outfit::observations::observations_ext::ObservationIOD; +use outfit::outfit::Outfit; +use outfit::outfit_errors::OutfitError; +use outfit::trajectories::trajectory_file::TrajectoryFile; +use outfit::TrajectorySet; +use rand::rngs::StdRng; +use rand::SeedableRng; + +/// Run Gauss Initial Orbit Determination on a single trajectory. +/// +/// This function retrieves the observation list for a given trajectory, +/// constructs reasonable IOD parameters (noise realizations, triplet limits, +/// etc.) and calls [`ObservationIOD::estimate_best_orbit`]. +/// +/// # Arguments +/// * `env_state` – Mutable Outfit runtime environment (ephemerides, settings…) +/// * `traj_set` – The global trajectory set loaded from the Parquet file +/// * `traj_number` – The ID of the trajectory to solve +/// +/// # Returns +/// A tuple `(GaussResult, rms)` if a valid orbit is found, or a wrapped +/// [`OutfitError`] otherwise. +/// +/// Note that failures are **expected** for some trajectories in survey-like +/// datasets (insufficient geometry, poor conditioning, etc.). +fn run_iod( + env_state: &mut Outfit, + traj_set: &mut TrajectorySet, + traj_number: &ObjectNumber, +) -> Result<(GaussResult, f64), OutfitError> { + // Retrieve the mutable list of observations for this trajectory. + let obs = traj_set.get_mut(traj_number).expect("trajectory not found"); + + // Use a deterministic RNG so the benchmark remains reproducible. + let mut rng = StdRng::seed_from_u64(42); + + // Build the parameters controlling the Gauss IOD solver. + let params = IODParams::builder() + .n_noise_realizations(10) + .noise_scale(1.1) + .max_obs_for_triplets(obs.len()) + .max_triplets(30) + .build()?; + + // Perform the actual orbit estimation. + obs.estimate_best_orbit(env_state, &ErrorModel::FCCT14, &mut rng, ¶ms) +} + +/// Prepare the Outfit environment and the full `TrajectorySet` from the +/// Fink-like Parquet test file. +/// +/// This function: +/// 1. Loads DE440 via Outfit, +/// 2. Initializes the observer (ZTF Palomar I41), +/// 3. Reads all trajectories from the Parquet file, +/// 4. Returns the environment, the trajectory set, and all available IDs. +/// +/// This preparation is done **once** and reused across benchmark iterations. +fn prepare_env_parquet() -> (Outfit, TrajectorySet, Vec) { + // Create the Outfit environment using DE440. + let mut env_state = + Outfit::new("horizon:DE440", ErrorModel::FCCT14).expect("Outfit init (DE440)"); + + let path_file = std::env::var("OUTFIT_IOD_BATCH_DIFF_PARQUET") + .unwrap_or_else(|_| "tests/data/test_from_fink.parquet".to_string()); + let path_file = Utf8Path::new(&path_file); + + // Using the Palomar/ZTF observer (MPC code I41). + let ztf_observer = env_state.get_observer_from_mpc_code(&"I41".into()); + + // Load all trajectories from the Parquet file. + // + // The last two parameters are the RA/Dec clipping radii used when building + // tangent-plane residuals. They are intentionally loose in this benchmark. + let traj_set = + TrajectorySet::new_from_parquet(&mut env_state, path_file, ztf_observer, 0.5, 0.5, None) + .expect("load TrajectorySet from Parquet"); + + // Collect all trajectory IDs. Assumes the underlying representation exposes `.keys()`. + let ids: Vec = traj_set.keys().cloned().collect(); + + (env_state, traj_set, ids) +} + +/// Benchmark the scaling of Gauss IOD across increasingly large batches of +/// **distinct trajectories**. +/// +/// Unlike micro-benchmarks which refit the same orbit repeatedly, this test +/// mimics a survey scenario: different objects, heterogeneous geometries, +/// and varying observation counts. +/// +/// This provides a more realistic estimation of the throughput one can expect +/// in daily survey operations. +fn bench_gauss_iod_parquet_batch_scaling(c: &mut Criterion) { + // Prepare environment and dataset once; reuse across all benchmark runs. + let (env0, traj0, all_ids) = prepare_env_parquet(); + let env = RefCell::new(env0); + let traj = RefCell::new(traj0); + + // Batch sizes (numbers of distinct trajectories per iteration). + // These remain moderate so the benchmark does not take too long. + let batch_sizes = [ + 1_usize, 4, 16, 64, 256, 512, 1024, 2048, 4096, 8192, 16384, //, 32768, 65536, 131072, + ]; + + let mut group = c.benchmark_group("gauss_iod_parquet_batch_scaling"); + group.sample_size(10); + + for &batch in &batch_sizes { + if batch > all_ids.len() { + continue; + } + + let name = format!("parquet_batch_{batch}"); + group.bench_function(&name, |b| { + b.iter_batched( + || (), // No heavy per-iteration setup needed + |_| { + let mut env_borrow = env.borrow_mut(); + let mut traj_borrow = traj.borrow_mut(); + + // Counters to ensure useful work is not optimized away. + let mut ok_count: usize = 0; + let mut err_count: usize = 0; + + // Fit orbits for the first `batch` trajectory IDs. + for id in all_ids.iter().take(batch) { + match run_iod(&mut env_borrow, &mut traj_borrow, id) { + Ok(result) => { + black_box(result); + ok_count += 1; + } + Err(err) => { + // Failures may occur depending on observation geometry. + black_box(&err); + err_count += 1; + } + } + } + + // Black-box the counters to prevent loop elimination. + black_box((ok_count, err_count)); + }, + BatchSize::SmallInput, + ) + }); + } + + group.finish(); +} + +#[cfg(feature = "parallel")] +mod iod_bench_parallel { + + use super::*; + use outfit::TrajectoryFit; + + /// Benchmark Gauss IOD using the parallel batch API over subsets of trajectories. + /// + /// For each `batch_size`, we build a `TrajectorySet` restricted to the first + /// `batch_size` objects (according to the ID list returned by `prepare_env_parquet`), + /// and run `estimate_all_orbits_in_batches_parallel` on that subset. + /// + /// This way, the x-axis ("batch size") represents the *number of trajectories* + /// processed, just like in the sequential benchmark. + pub fn bench_gauss_iod_parquet_parallel_batches(c: &mut Criterion) { + // Prepare environment, full dataset, and the full list of IDs once. + let (env0, traj0, all_ids) = prepare_env_parquet(); + let env = RefCell::new(env0); + + // Same batch sizes as the sequential benchmark so that comparisons are easy. + let batch_sizes = [ + 1_usize, 4, 16, 64, 256, 512, 1024, 2048, 4096, 8192, + 16384, // 32768, 65536, 131072, + ]; + + let mut group = c.benchmark_group("gauss_iod_parquet_parallel_batches"); + group.sample_size(10); + + for &batch_size in &batch_sizes { + // Skip if we don't have enough trajectories. + if batch_size > all_ids.len() { + continue; + } + + let name = format!("parquet_parallel_batch_{batch_size}"); + + // ---- Build a subset TrajectorySet with exactly `batch_size` trajectories ---- + // + // This is done once, *outside* Criterion's inner loop so that we are + // measuring the IOD cost only (not the cost of building the subset). + let selected_ids: Vec = + all_ids.iter().take(batch_size).cloned().collect(); + + let selected_set: std::collections::HashSet = + selected_ids.iter().cloned().collect(); + + // Clone the full set and retain only the selected trajectories. + let mut traj_subset = traj0.clone(); + traj_subset.retain(|obj_id, _obs| selected_set.contains(obj_id)); + + let traj = RefCell::new(traj_subset); + + // ---- Benchmark on the restricted subset ---- + group.bench_function(&name, |b| { + b.iter_batched( + || { + // Per-iteration setup: deterministic RNG + IOD parameters. + let rng = StdRng::seed_from_u64(42); + + let chunk_size = batch_size.min(64); + let params = IODParams::builder() + .n_noise_realizations(10) + .noise_scale(1.1) + .max_triplets(30) + .batch_size(chunk_size) + .build() + .expect("build IODParams"); + + (rng, params) + }, + |(mut rng, params)| { + let env_borrow = env.borrow(); + let mut traj_borrow = traj.borrow_mut(); + + let full_result = traj_borrow.estimate_all_orbits_in_batches_parallel( + &env_borrow, + &mut rng, + ¶ms, + ); + + // Black-box the result to avoid dead-code elimination. + black_box(full_result); + }, + BatchSize::SmallInput, + ) + }); + } + + group.finish(); + } +} + +#[cfg(feature = "parallel")] +criterion_group!( + gauss_parquet_benches, + bench_gauss_iod_parquet_batch_scaling, + iod_bench_parallel::bench_gauss_iod_parquet_parallel_batches, +); + +#[cfg(not(feature = "parallel"))] +criterion_group!(gauss_parquet_benches, bench_gauss_iod_parquet_batch_scaling,); + +criterion_main!(gauss_parquet_benches); diff --git a/benches/outfit_gauss_iod.rs b/benches/outfit_gauss_iod.rs index d9acdc2..2795862 100644 --- a/benches/outfit_gauss_iod.rs +++ b/benches/outfit_gauss_iod.rs @@ -18,7 +18,29 @@ mod benches_impl { use rand::rngs::StdRng; use rand::SeedableRng; - /// Run Gauss IOD on a single trajectory and return the best orbit + RMS. + /// Run Gauss Initial Orbit Determination on a single trajectory. + /// + /// This function: + /// - retrieves the mutable observation list for the given trajectory, + /// - builds a set of IOD parameters (noise realizations, triplet limits, …), + /// - calls [`ObservationIOD::estimate_best_orbit`] and returns the result. + /// + /// Arguments + /// ----------------- + /// * `env_state` - Mutable Outfit environment (DE ephemerides, settings, etc.). + /// * `traj_set` - Global trajectory set containing the observations. + /// * `traj_number` - ID of the trajectory to solve. + /// + /// Return + /// ---------- + /// * A `Result` containing `(GaussResult, rms)` if a valid orbit is found, + /// or an [`OutfitError`] if the IOD fails. + /// + /// See also + /// ------------ + /// * [`prepare_env`] - Helper to initialize the environment and trajectories. + /// * [`bench_gauss_iod_e2e`] - End-to-end benchmark on three trajectories. + /// * [`bench_gauss_iod_batch_scaling`] - Batch-scaling benchmark on the same set. fn run_iod( env_state: &mut Outfit, traj_set: &mut TrajectorySet, @@ -37,7 +59,34 @@ mod benches_impl { obs.estimate_best_orbit(env_state, &ErrorModel::FCCT14, &mut rng, ¶ms) } - /// Prepare environment with DE440 and 3 trajectories. + /// Prepare the Outfit environment (DE440) and a small set of test trajectories. + /// + /// This helper: + /// - instantiates an [`Outfit`] environment using the DE440 ephemeris, + /// - loads three 80-column observation files into a [`TrajectorySet`], + /// - returns the environment, the trajectory set, and the three IDs. + /// + /// The three objects are: + /// - `K09R05F` (corresponding to `2015AB.obs`), + /// - `8467`, + /// - `33803`. + /// + /// Arguments + /// ----------------- + /// * None. + /// + /// Return + /// ---------- + /// * A tuple `(env, set, ids)` where: + /// - `env` is the initialized [`Outfit`] environment, + /// - `set` is the [`TrajectorySet`] populated with the three objects, + /// - `ids` is the fixed array of [`ObjectNumber`] identifiers. + /// + /// See also + /// ------------ + /// * [`run_iod`] - Core IOD routine using this prepared environment. + /// * [`bench_gauss_iod_core`] - Reuses `env` but reloads trajectories each iteration. + /// * [`bench_gauss_iod_batch_scaling`] - Uses these same objects in round-robin. fn prepare_env() -> (Outfit, TrajectorySet, [ObjectNumber; 3]) { let mut env = Outfit::new("horizon:DE440", ErrorModel::FCCT14).expect("Outfit init"); @@ -55,7 +104,19 @@ mod benches_impl { (env, set, ids) } - /// End-to-end benchmark. + /// End-to-end benchmark on three real trajectories (full pipeline). + /// + /// This benchmark measures the cost of running Gauss IOD on all three + /// trajectories, including environment and trajectory set initialization. + /// + /// Arguments + /// ----------------- + /// * `c` - Criterion handle used to register the benchmark. + /// + /// See also + /// ------------ + /// * [`bench_gauss_iod_core`] - Reuses the environment and reloads data per iteration. + /// * [`bench_gauss_iod_batch_scaling`] - Scales the number of orbit fits per batch. pub fn bench_gauss_iod_e2e(c: &mut Criterion) { c.bench_function("gauss_iod_e2e_all", |b| { b.iter_batched( @@ -71,7 +132,25 @@ mod benches_impl { }); } - /// Core benchmark. + /// Core Gauss IOD benchmark: reuse environment, reload one trajectory per iteration. + /// + /// This benchmark keeps the [`Outfit`] environment in a shared [`RefCell`], + /// and repeatedly: + /// - rebuilds a `TrajectorySet` from `2015AB.obs`, + /// - runs Gauss IOD only on `K09R05F`. + /// + /// This isolates the cost of: + /// - loading a single trajectory from 80-column format, + /// - performing one orbit determination for this trajectory. + /// + /// Arguments + /// ----------------- + /// * `c` - Criterion handle used to register the benchmark. + /// + /// See also + /// ------------ + /// * [`bench_gauss_iod_e2e`] - End-to-end benchmark over three trajectories. + /// * [`bench_gauss_iod_batch_scaling`] - Varies the number of fits per batch. pub fn bench_gauss_iod_core(c: &mut Criterion) { let (env0, _set_once, _ids) = prepare_env(); let env = RefCell::new(env0); @@ -93,6 +172,42 @@ mod benches_impl { ) }); } + + /// Batch-scaling benchmark: N orbit fits per batch using the same 3 trajectories. + /// + /// This benchmark is useful to build "time per orbit-fit vs batch size" + /// plots. It reuses the three real objects (K09R05F, 8467, 33803) and + /// cycles through them in round-robin fashion for each batch. + /// + /// Arguments + /// ----------------- + /// * `c` - Criterion handle used to register the benchmark. + /// + /// See also + /// ------------ + /// * [`bench_gauss_iod_e2e`] - Baseline for full end-to-end performance. + /// * [`bench_gauss_iod_core`] - Focuses on a single trajectory reload + fit. + pub fn bench_gauss_iod_batch_scaling(c: &mut Criterion) { + // Test several batch sizes: reuse the same three objects in round-robin. + let batch_sizes = [1_usize, 2, 3, 6, 12, 24, 48, 96]; + + for &batch in &batch_sizes { + let name = format!("gauss_iod_batch_{batch}"); + c.bench_function(&name, |b| { + b.iter_batched( + prepare_env, + |(mut env, mut set, ids)| { + for i in 0..batch { + let id = &ids[i % ids.len()]; + let res = run_iod(&mut env, &mut set, id).expect("IOD"); + black_box(res); + } + }, + BatchSize::SmallInput, + ) + }); + } + } } #[cfg(feature = "jpl-download")] @@ -106,13 +221,15 @@ criterion_group!( .warm_up_time(std::time::Duration::from_secs(5)) .measurement_time(std::time::Duration::from_secs(25)) .with_plots(); - targets = benches_impl::bench_gauss_iod_e2e, benches_impl::bench_gauss_iod_core + targets = + benches_impl::bench_gauss_iod_e2e, + benches_impl::bench_gauss_iod_core, + benches_impl::bench_gauss_iod_batch_scaling ); #[cfg(feature = "jpl-download")] criterion_main!(benches); -// Fallback quand la feature est absente : fournit une main au crate. #[cfg(not(feature = "jpl-download"))] fn main() { eprintln!( diff --git a/joss_data/paper.bib b/joss_data/paper.bib new file mode 100644 index 0000000..160b764 --- /dev/null +++ b/joss_data/paper.bib @@ -0,0 +1,289 @@ +@misc{rubinobservatorylsstsolarsystemsciencecollaboration2020scientificimpactverac, + title = {The Scientific Impact of the Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST) for Solar System Science}, + author = { LSST-SSSC and R. Lynne Jones and Michelle T. Bannister and Bryce T. Bolin and Colin Orion Chandler and Steven R. Chesley and Siegfried Eggl and Sarah Greenstreet and Timothy R. Holt and Henry H. Hsieh and Zeljko Ivezić and Mario Jurić and Michael S. P. Kelley and Matthew M. Knight and Renu Malhotra and William J. Oldroyd and Gal Sarid and Megan E. Schwamb and Colin Snodgrass and Michael Solontoi and David E. Trilling}, + year = {2020}, + eprint = {2009.07653}, + archiveprefix = {arXiv}, + primaryclass = {astro-ph.IM}, + url = {https://arxiv.org/abs/2009.07653} +} + +@article{LSST_from_science_to_design, + author = {{Ivezi{\'c}}, {\v{Z}}eljko and {Kahn}, Steven M. and {Tyson}, J. Anthony and {Abel}, Bob and {Acosta}, Emily and {Allsman}, Robyn and {Alonso}, David and {AlSayyad}, Yusra and {Anderson}, Scott F. and {Andrew}, John and {Angel}, James Roger P. and {Angeli}, George Z. and {Ansari}, Reza and {Antilogus}, Pierre and {Araujo}, Constanza and {Armstrong}, Robert and {Arndt}, Kirk T. and {Astier}, Pierre and {Aubourg}, {\'E}ric and {Auza}, Nicole and {Axelrod}, Tim S. and {Bard}, Deborah J. and {Barr}, Jeff D. and {Barrau}, Aurelian and {Bartlett}, James G. and {Bauer}, Amanda E. and {Bauman}, Brian J. and {Baumont}, Sylvain and {Bechtol}, Ellen and {Bechtol}, Keith and {Becker}, Andrew C. and {Becla}, Jacek and {Beldica}, Cristina and {Bellavia}, Steve and {Bianco}, Federica B. and {Biswas}, Rahul and {Blanc}, Guillaume and {Blazek}, Jonathan and {Blandford}, Roger D. and {Bloom}, Josh S. and {Bogart}, Joanne and {Bond}, Tim W. and {Booth}, Michael T. and {Borgland}, Anders W. and {Borne}, Kirk and {Bosch}, James F. and {Boutigny}, Dominique and {Brackett}, Craig A. and {Bradshaw}, Andrew and {Brandt}, William Nielsen and {Brown}, Michael E. and {Bullock}, James S. and {Burchat}, Patricia and {Burke}, David L. and {Cagnoli}, Gianpietro and {Calabrese}, Daniel and {Callahan}, Shawn and {Callen}, Alice L. and {Carlin}, Jeffrey L. and {Carlson}, Erin L. and {Chandrasekharan}, Srinivasan and {Charles-Emerson}, Glenaver and {Chesley}, Steve and {Cheu}, Elliott C. and {Chiang}, Hsin-Fang and {Chiang}, James and {Chirino}, Carol and {Chow}, Derek and {Ciardi}, David R. and {Claver}, Charles F. and {Cohen-Tanugi}, Johann and {Cockrum}, Joseph J. and {Coles}, Rebecca and {Connolly}, Andrew J. and {Cook}, Kem H. and {Cooray}, Asantha and {Covey}, Kevin R. and {Cribbs}, Chris and {Cui}, Wei and {Cutri}, Roc and {Daly}, Philip N. and {Daniel}, Scott F. and {Daruich}, Felipe and {Daubard}, Guillaume and {Daues}, Greg and {Dawson}, William and {Delgado}, Francisco and {Dellapenna}, Alfred and {de Peyster}, Robert and {de Val-Borro}, Miguel and {Digel}, Seth W. and {Doherty}, Peter and {Dubois}, Richard and {Dubois-Felsmann}, Gregory P. and {Durech}, Josef and {Economou}, Frossie and {Eifler}, Tim and {Eracleous}, Michael and {Emmons}, Benjamin L. and {Fausti Neto}, Angelo and {Ferguson}, Henry and {Figueroa}, Enrique and {Fisher-Levine}, Merlin and {Focke}, Warren and {Foss}, Michael D. and {Frank}, James and {Freemon}, Michael D. and {Gangler}, Emmanuel and {Gawiser}, Eric and {Geary}, John C. and {Gee}, Perry and {Geha}, Marla and {Gessner}, Charles J.~B. and {Gibson}, Robert R. and {Gilmore}, D. Kirk and {Glanzman}, Thomas and {Glick}, William and {Goldina}, Tatiana and {Goldstein}, Daniel A. and {Goodenow}, Iain and {Graham}, Melissa L. and {Gressler}, William J. and {Gris}, Philippe and {Guy}, Leanne P. and {Guyonnet}, Augustin and {Haller}, Gunther and {Harris}, Ron and {Hascall}, Patrick A. and {Haupt}, Justine and {Hernandez}, Fabio and {Herrmann}, Sven and {Hileman}, Edward and {Hoblitt}, Joshua and {Hodgson}, John A. and {Hogan}, Craig and {Howard}, James D. and {Huang}, Dajun and {Huffer}, Michael E. and {Ingraham}, Patrick and {Innes}, Walter R. and {Jacoby}, Suzanne H. and {Jain}, Bhuvnesh and {Jammes}, Fabrice and {Jee}, M. James and {Jenness}, Tim and {Jernigan}, Garrett and {Jevremovi{\'c}}, Darko and {Johns}, Kenneth and {Johnson}, Anthony S. and {Johnson}, Margaret W.~G. and {Jones}, R. Lynne and {Juramy-Gilles}, Claire and {Juri{\'c}}, Mario and {Kalirai}, Jason S. and {Kallivayalil}, Nitya J. and {Kalmbach}, Bryce and {Kantor}, Jeffrey P. and {Karst}, Pierre and {Kasliwal}, Mansi M. and {Kelly}, Heather and {Kessler}, Richard and {Kinnison}, Veronica and {Kirkby}, David and {Knox}, Lloyd and {Kotov}, Ivan V. and {Krabbendam}, Victor L. and {Krughoff}, K. Simon and {Kub{\'a}nek}, Petr and {Kuczewski}, John and {Kulkarni}, Shri and {Ku}, John and {Kurita}, Nadine R. and {Lage}, Craig S. and {Lambert}, Ron and {Lange}, Travis and {Langton}, J. Brian and {Le Guillou}, Laurent and {Levine}, Deborah and {Liang}, Ming and {Lim}, Kian-Tat and {Lintott}, Chris J. and {Long}, Kevin E. and {Lopez}, Margaux and {Lotz}, Paul J. and {Lupton}, Robert H. and {Lust}, Nate B. and {MacArthur}, Lauren A. and {Mahabal}, Ashish and {Mandelbaum}, Rachel and {Markiewicz}, Thomas W. and {Marsh}, Darren S. and {Marshall}, Philip J. and {Marshall}, Stuart and {May}, Morgan and {McKercher}, Robert and {McQueen}, Michelle and {Meyers}, Joshua and {Migliore}, Myriam and {Miller}, Michelle and {Mills}, David J.}, + title = {{LSST: From Science Drivers to Reference Design and Anticipated Data Products}}, + journal = {\apj}, + keywords = {astrometry, cosmology: observations, Galaxy: general, methods: observational, stars: general, surveys, Astrophysics}, + year = 2019, + month = mar, + volume = {873}, + number = {2}, + eid = {111}, + pages = {111}, + doi = {10.3847/1538-4357/ab042c}, + archiveprefix = {arXiv}, + eprint = {0805.2366}, + primaryclass = {astro-ph}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2019ApJ...873..111I}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@article{Fink_2020, + title = {fink, a new generation of broker for the LSST community}, + volume = {501}, + issn = {1365-2966}, + url = {http://dx.doi.org/10.1093/mnras/staa3602}, + doi = {10.1093/mnras/staa3602}, + number = {3}, + journal = {Monthly Notices of the Royal Astronomical Society}, + publisher = {Oxford University Press (OUP)}, + author = {Möller, Anais and Peloton, Julien and Ishida, Emille E O and Arnault, Chris and Bachelet, Etienne and Blaineau, Tristan and Boutigny, Dominique and Chauhan, Abhishek and Gangler, Emmanuel and Hernandez, Fabio and Hrivnac, Julius and Leoni, Marco and Leroy, Nicolas and Moniez, Marc and Pateyron, Sacha and Ramparison, Adrien and Turpin, Damien and Ansari, Réza and Allam Jr, Tarek and Bajat, Armelle and Biswas, Biswajit and Boucaud, Alexandre and Bregeon, Johan and Campagne, Jean-Eric and Cohen-Tanugi, Johann and Coleiro, Alexis and Dornic, Damien and Fouchez, Dominique and Godet, Olivier and Gris, Philippe and Karpov, Sergey and Nebot Gomez-Moran, Ada and Neveu, Jérémy and Plaszczynski, Stephane and Savchenko, Volodymyr and Webb, Natalie}, + year = {2020}, + month = nov, + pages = {3272–3288} +} + +@software{OrbFit, + author = {{Orbfit Consortium}}, + title = "{OrbFit: Software to Determine Orbits of Asteroids}", + howpublished = {Astrophysics Source Code Library, record ascl:1106.015}, + year = 2011, + month = jun, + eid = {ascl:1106.015}, +archivePrefix = {ascl}, + eprint = {1106.015}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2011ascl.soft06015O}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} + +@misc{FindOrb, + author = {Gray, Bill}, + title = {Find\_Orb: Orbit Determination Software}, + howpublished = {\url{https://www.projectpluto.com/find_orb.htm}}, + note = {Accessed: 2025-02-12}, + year = {2000} +} + +@article{Berres_2022, +doi = {10.3847/2515-5172/ac8e08}, +url = {https://doi.org/10.3847/2515-5172/ac8e08}, +year = {2022}, +month = {sep}, +publisher = {The American Astronomical Society}, +volume = {6}, +number = {9}, +pages = {174}, +author = {Berres, Aidan and Greenstreet, Sarah and Moeyens, Joachim and Jurić, Mario and Lu, Edward}, +title = {Orb_It: A Validation Package for Orbit Integrators}, +journal = {Research Notes of the AAS}, +abstract = {We present a testing platform and package that can validate the accuracy of an orbit integrator. Orb_It is written in Python and is capable of validating the most heavily used orbital integrators such as Find_Orb, OpenOrb, and OrbFit. Orb_It performs an “end-to-end” test of each integrator to check for internal consistency of the integrators functions. For a given objects state vector and epoch, Orb_It uses a user-specified integrator to generate the objects ephemeris for a specified time range and observer, uses those ephemerides as observations to fit an orbit, and then compares that fitted orbit to the original input state vector propagated to the fitted orbits epoch. This allows users to understand the internal accuracy of each integrator. This testing suite and package can be found on GitHub with a demo on Binder (https://github.com/B612-Asteroid-Institute/orb_it).} +} + + +@misc{MPC80col, + title = {Minor Planet Center Optical Observations: 80-column format}, + author = {Minor Planet Center}, + year = {2025}, + howpublished = {Web page}, + url = {https://minorplanetcenter.net/iau/info/OpticalObs.html}, + note = {Accessed: 2025-12-05} +} + +@misc{ADES, + title = {Astrometric Data Exchange Standard (ADES)}, + author = {IAU Minor Planet Center}, + year = {2025}, + howpublished = {Web page}, + url = {https://minorplanetcenter.net/iau/info/ADES.html}, + note = {Accessed: 2025-12-05} +} + +@article{DE440, + author = {Park, R. S. and Folkner, W. M. and Jacobson, R. A. and others}, + title = {The {JPL} Planetary and Lunar Ephemerides {DE440} and {DE441}}, + journal = {The Astronomical Journal}, + year = {2021}, + volume = {161}, + number = {3}, + pages = {105}, + doi = {10.3847/1538-3881/abd414} +} + +@misc{NAIF, + author = {Acton, C. H. and {NAIF/NASA JPL}}, + title = {Navigation and Ancillary Information Facility (SPICE)}, + year = {2025}, + howpublished = {Web page}, + url = {https://naif.jpl.nasa.gov/naif/}, + note = {Accessed: 2025-12-05} +} + +@article{Arrow, + author = {P. E. {Apache Arrow Project}}, + title = {Apache Arrow: A Cross-Language Development Platform for In-Memory Data}, + journal = {Data Engineering}, + year = {2017}, + note = {Project documentation: https://arrow.apache.org/} +} + +@misc{Parquet, + author = {{Apache Parquet Project}}, + title = {Apache Parquet}, + year = {2025}, + howpublished = {Project page}, + url = {https://parquet.apache.org/}, + note = {Accessed: 2025-12-05} +} + +@misc{rust, + author = {{Rust Project}}, + title = {The Rust Programming Language}, + year = {2025}, + howpublished = {Official website}, + url = {https://www.rust-lang.org/}, + note = {Accessed: 2025-12-05} +} + +@book{rust_book, + author = {Klabnik, Steve and Nichols, Carol}, + title = {The Rust Programming Language}, + year = {2018}, + publisher = {No Starch Press}, + isbn = {978-1-59327-828-1}, + howpublished = {Online book}, + url = {https://doc.rust-lang.org/book/}, + note = {Accessed: 2026-04-01} +} + +@misc{nalgebra, + author = {Crozet, S. and Dimforge Contributors}, + title = {nalgebra: Linear algebra library for Rust}, + year = {2025}, + howpublished = {Crate documentation}, + url = {https://crates.io/crates/nalgebra}, + note = {Accessed: 2025-12-05} +} + +@misc{hifitime, + author = {Christopher Rabotin}, + title = {hifitime: High-Fidelity Time library for Rust}, + year = {2025}, + howpublished = {Crate documentation}, + url = {https://crates.io/crates/hifitime}, + note = {Accessed: 2025-12-05} +} + +@misc{rayon, + author = {N. Matsakis and {Rayon Contributors}}, + title = {Rayon: Data Parallelism in Rust}, + year = {2025}, + howpublished = {Crate documentation}, + url = {https://crates.io/crates/rayon}, + note = {Accessed: 2025-12-05} +} + +@misc{Aberth, + author = {ickk and contributors}, + title = {aberth: Aberth–Ehrlich polynomial root-finding in Rust}, + year = {2025}, + howpublished = {GitHub repository}, + url = {https://github.com/ickk/aberth}, + note = {Accessed: 2025-12-05} +} + +@misc{indicatif, + author = {F. {Indicatif Contributors}}, + title = {indicatif: Progress bars for Rust}, + howpublished = {Crate documentation}, + url = {https://crates.io/crates/indicatif}, + note = {Accessed: 2025-12-05} +} + +@misc{comfytable, + author = {{comfy-table Contributors}}, + title = {comfy-table: Table formatting for Rust}, + howpublished = {Crate documentation}, + url = {https://crates.io/crates/comfy-table}, + note = {Accessed: 2025-12-05} +} + +@misc{SemVer, + author = {Preston-Werner, T.}, + title = {Semantic Versioning 2.0.0}, + year = {2025}, + howpublished = {Specification}, + url = {https://semver.org/}, + note = {Accessed: 2025-12-05} +} + +@misc{KeepAChangelog, + author = {Oliveira, O. and others}, + title = {Keep a Changelog}, + year = {2025}, + howpublished = {Project page}, + url = {https://keepachangelog.com/en/1.0.0/}, + note = {Accessed: 2025-12-05} +} + +@software{pyo3, + author = {David Hewitt and Nikolay Kim and the pyo3 Admins}, + title = {{pyo3}: Rust bindings for Python}, + url = {https://crates.io/crates/pyo3}, + version = {0.27.2}, + description = {Rust bindings to the Python interpreter, enabling creation of native Python modules or embedding Python in Rust.}, + date = {2025}, + note = {Owners: Admins (github:pyo3:admins), David Hewitt, Nikolay Kim; Accessed: 2026-01-08}, +} + +@software{maturin, + author = {PyO3 Admins and messense and konsti}, + title = {{maturin}: Build and publish Rust-based Python packages}, + url = {https://crates.io/crates/maturin}, + version = {1.11.0}, + description = {Command-line tool to build, package, and publish Rust crates as Python packages (wheels) with support for pyo3 and other bindings.}, + date = {2025}, + note = {Owners: Admins (github:pyo3:admins), messense, konsti; Accessed: 2026-01-08}, +} + +@Inbook{Williams1997, +author="Williams, Gareth V.", +title="V{\"a}is{\"a}l{\"a} orbit", +bookTitle="Encyclopedia of Planetary Science", +year="1997", +publisher="Springer Netherlands", +address="Dordrecht", +pages="875--876", +isbn="978-1-4020-4520-2", +doi="10.1007/1-4020-4520-4_432", +url="https://doi.org/10.1007/1-4020-4520-4_432" +} + +@article{Le_Montagner_2023, + title={Enabling discoveries of Solar System objects in large alert data streams}, + volume={680}, + ISSN={1432-0746}, + url={http://dx.doi.org/10.1051/0004-6361/202346905}, + DOI={10.1051/0004-6361/202346905}, + journal={Astronomy & Astrophysics}, + publisher={EDP Sciences}, + author={Le Montagner, R. and Peloton, J. and Carry, B. and Desmars, J. and Hestroffer, D. and Mendez, R. A. and Perlbarg, A. C. and Thuillot, W.}, + year={2023}, + month=dec, pages={A17} } + +@article{10.1145/2934664, +author = {Zaharia, Matei and Xin, Reynold S. and Wendell, Patrick and Das, Tathagata and Armbrust, Michael and Dave, Ankur and Meng, Xiangrui and Rosen, Josh and Venkataraman, Shivaram and Franklin, Michael J. and Ghodsi, Ali and Gonzalez, Joseph and Shenker, Scott and Stoica, Ion}, +title = {Apache Spark: a unified engine for big data processing}, +year = {2016}, +issue_date = {November 2016}, +publisher = {Association for Computing Machinery}, +address = {New York, NY, USA}, +volume = {59}, +number = {11}, +issn = {0001-0782}, +url = {https://doi.org/10.1145/2934664}, +doi = {10.1145/2934664}, +abstract = {This open source computing framework unifies streaming, batch, and interactive big data workloads to unlock new applications.}, +journal = {Commun. ACM}, +month = oct, +pages = {56–65}, +numpages = {10} +} \ No newline at end of file diff --git a/joss_data/paper.md b/joss_data/paper.md new file mode 100644 index 0000000..1b72ad4 --- /dev/null +++ b/joss_data/paper.md @@ -0,0 +1,100 @@ +--- +title: "Outfit: a fast, safe Rust library for astrometric observations and preliminary orbit determination" +tags: + - Rust + - astronomy + - orbit determination + - astrometry + - ephemerides +authors: + - name: Roman Le Montagner + orcid: 0000-0002-6099-8939 + affiliation: "1" + - name: Hadrien Grasland + orcid: 0000-0003-3379-7100 + affiliation: "1" + - name: Julien Peloton + orcid: 0000-0002-8560-4449 + affiliation: "1" +affiliations: + - index: 1 + name: Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France +date: 5 December 2025 +bibliography: paper.bib +--- + +# Summary + +The NSF-DOE Vera C. Rubin Observatory [@LSST_from_science_to_design] conducts the Legacy Survey of Space and Time (LSST), which produces an unprecedented optical alert stream exceeding ten million alerts per night. A substantial fraction of these alerts corresponds to moving Solar System objects. LSST is expected to increase the number of known asteroids by at least an order of magnitude—from roughly 1.3 million today to several million new discoveries over the survey lifetime [@rubinobservatorylsstsolarsystemsciencecollaboration2020scientificimpactverac]. Handling this volume requires fast and reliable tools to read astrometric observations and estimate preliminary orbits automatically. + +Alert brokers such as Fink [@Fink_2020] receive the LSST alert stream ahead of all other users and classify alerts in real time. Their position at the earliest stage of data dissemination places them on the front line for recognising and flagging previously unknown asteroids. Early identification of such objects improves alert classification, enables rapid computation of ephemerides for follow-up, and supports linking future detections to newly derived orbits. It also allows brokers to filter out moving objects when searching for fast optical transients, an important capability for multi-messenger astronomy. + +**Outfit[^outfit]** is a software library that performs these calculations automatically. It reads standardised observation records, applies classical orbit-determination algorithms, and returns preliminary orbital solutions that can be used for object identification, follow-up planning, and alert classification. Outfit is designed to be embedded directly into automated data-processing pipelines, enabling real-time analysis of large-scale survey data. + +[^outfit]: + +# Statement of need + +Preliminary orbit determination is a foundational step in Solar System science. Researchers, survey teams, alert brokers, and even advanced amateur observers rely on software that can ingest astrometric measurements and rapidly estimate orbits for newly detected objects. + +Modern surveys operate at unprecedented scales, and their data-processing infrastructures require software components that are composable, reliable, and easy to integrate programmatically. Developers building next-generation pipelines need open-source components that are composable, reproducible, and easy to embed within larger scientific software stacks—whether for nightly survey operations, real-time alert classification, or educational and citizen-science applications. + +Outfit addresses this need by providing a modern library interface for preliminary orbit determination. It allows users to process observations in memory and estimate orbits programmatically, making classical orbit-determination methods accessible and usable within contemporary survey pipelines and alert-broker infrastructures. This design is particularly important for real-time systems such as Fink [@Fink_2020], which receive the LSST alert stream ahead of all other users and must classify alerts within seconds. Early identification of moving objects improves alert classification, enables rapid computation of ephemerides for follow-up, supports linking future detections to newly derived orbits, and allows brokers to filter out asteroids when searching for fast optical transients—an important capability for multi-messenger astronomy. + +# State of the field + +Within the astronomical community, tools such as OrbFit [@OrbFit] and Find_Orb [@FindOrb] have long played an essential role by providing scientifically trusted implementations of classical orbit-determination algorithms. OrbFit, developed in Fortran, is a mature and rigorously validated tool that has been widely adopted in observatories and research institutions. Find_Orb, implemented in C, is similarly well established and provides comparable functionality. Both tools are distributed as standalone executables that read input files and write output files, a design pattern that has served the community well for decades. Recent validation efforts such as Orb_It [@Berres_2022] provide end-to-end testing frameworks for these tools, confirming their internal consistency and accuracy. + +However, these tools were not designed for programmatic integration into automated pipelines. Invoking them from within a larger application requires launching external processes, writing temporary input files, parsing textual output, and managing inter-process communication—all sources of complexity and runtime overhead that become significant bottlenecks at the scales required by modern surveys. Furthermore, neither tool exposes a library API that would allow direct in-memory processing. + +While OrbFit and Find_Orb remain the gold standard for scientific accuracy and validation, their executable-based design makes them ill-suited for the use cases Outfit targets. Rather than contribute additional features to these existing tools, we chose to build a new library from the ground up with a different architectural philosophy: an embeddable, memory-safe library with a clear API designed for automated workflows. This design choice represents a meaningful scholarly contribution, as it enables classical orbit-determination methods to be integrated into modern data-processing systems that would otherwise struggle to adopt executable-based tools. + +# Software design + +The architecture of Outfit reflects a deliberate set of trade-offs made to balance scientific accuracy, computational performance, and ease of integration. The most fundamental design decision was to implement Outfit as a library rather than an executable. This choice prioritises programmatic access over interactive use, enabling developers to embed orbit-determination capabilities directly into their applications without inter-process communication overhead. + +Outfit is implemented in Rust [@rust_book], a programming language that provides memory safety without garbage collection and zero-cost abstractions for performance-critical code. This choice was motivated by the need for both speed and reliability in a library intended for automated, high-throughput pipelines. Rust's ownership system eliminates entire classes of runtime errors (null pointer dereferences, use-after-free, data races) without imposing runtime overhead, making it well-suited for scientific computing applications that must be both fast and correct. + +The library accepts observations directly as in-memory data structures, avoiding file I/O in the critical path. For interoperability with existing workflows, Outfit also supports standard formats in astronomy including MPC 80-column [@MPC80col], ADES XML [@ADES], and columnar layouts such as Parquet/Arrow [@Arrow; @Parquet]. These parsers perform validation and time-scale normalisation on load, ensuring that downstream orbit-fitting code operates on consistently formatted data. + +Initial orbit determination (IOD) follows the classical Gauss method on triplets of optical observations. The solver performs iterative refinement with numerical and geometric acceptability checks, evaluating RMS residuals to assess fit quality. The eighth-degree polynomial arising in Gauss's formulation is solved using an Aberth–Ehrlich root-finding method (via the aberth crate [@Aberth]), which provides numerical stability and efficiency. Observatory geometry is handled explicitly via MPC codes or geodetic coordinates, converted into consistent geocentric and topocentric reference frames. Outfit supports multiple orbital element representations (Keplerian, equinoctial, and cometary), two-body propagation, and planetary ephemerides from JPL DE kernels (e.g., DE440) [@DE440; @NAIF]. Time handling uses the hifitime library [@hifitime] for consistent UTC/TT/TDB conversions. + +For survey-scale workloads, Outfit provides batched and optionally parallel processing of trajectories. Parallelism is implemented with Rayon [@rayon], allowing multi-core evaluation of independent orbit fits while preserving deterministic behaviour when required. This design lets users tune memory use and throughput effectively across a wide range of data volumes. Core numerical operations rely on nalgebra [@nalgebra] for efficient linear-algebra routines, leveraging a mature and well-tested foundation for vector and matrix operations. + +To maximise accessibility, Outfit is also exposed to Python via **pyOutfit[^pyoutfit]**, a binding built with PyO3 [@pyo3] and distributed with maturin [@maturin]. This provides a thin, idiomatic Python interface to the Rust core, enabling integration in notebooks and Python-based pipelines. The binding is available on PyPI and ships precompiled manylinux wheels for Python 3.10, 3.11, and 3.12 at the time of writing. Project documentation is published on GitHub Pages using MkDocs Material, and the public API is type-annotated via .pyi stub files to improve IDE completion and static type checking. + +[^pyoutfit]: + +These design choices matter for research applications because they enable orbit determination to be embedded directly into real-time alert classification systems (such as Fink), distributed computing frameworks (such as Apache Spark), and other automated pipelines where launching external executables would introduce unacceptable latency or complexity. By providing a library interface with low integration friction, Outfit makes classical orbit-determination methods accessible to a broader range of research workflows. + +# Research impact statement + +Outfit has demonstrated realized impact through integration into operational research infrastructure and scholarly publication. The library is actively used within the Fink broker [@Fink_2020], which processes the LSST alert stream in real time and serves as a critical component of the survey's early data-dissemination pipeline. This operational deployment validates Outfit's design for high-throughput, low-latency applications. + +The research need that motivated Outfit's development is documented in [@Le_Montagner_2023], which used OrbFit to fit tens of thousands of trajectory hypotheses by invoking the Fortran executable through Apache Spark across multiple compute nodes—a workflow that required dozens of CPU cores and took several minutes to complete. With Outfit's library interface, the same analysis runs in seconds on a single multi-core machine, demonstrating the performance gains enabled by in-memory processing and eliminating the overhead of inter-process communication and file I/O. + +To ensure numerical reliability, Outfit includes regression tests on real astrometric datasets (2015AB[^2015AB], 8467[^8467], 33803[^33803]) that validate orbit-fitting accuracy against known solutions returned by OrbFit. These tests run automatically in continuous integration and serve as a reproducible benchmark for the library's correctness. The library is distributed on crates.io (the Rust package registry) and PyPI (the Python Package Index), with precompiled binaries for multiple platforms, lowering the barrier to adoption. Documentation is available on docs.rs and GitHub Pages, and the repository includes runnable examples demonstrating end-to-end usage: + +```bash +cargo run --release --example parquet_to_orbit --features "jpl-download" +``` + +[^2015AB]: +[^8467]: +[^33803]: + +# AI usage disclosure + +GitHub Copilot's auto-completion feature was used in VS Code during implementation, and the GitHub Copilot coding agent assisted with writing testing infrastructure, documentation, and refactoring. AI-assisted language editing was used to improve wording and clarity in parts of the manuscript. All code and text generated or modified by AI tools was carefully reviewed and proofread by the authors, and all AI-generated code was extensively tested. + +# Limitations and future work + +Outfit currently focuses on optical astrometry and provides initial orbit determination through the classical Gauss method. It does not yet support radar or radio observables, hyperbolic orbit fitting, or full least-squares refinement, and therefore cannot produce formal covariances or robust uncertainty estimates. In addition, while planetary ephemerides are fully integrated, the library does not yet expose a high-level API for generating predicted ephemerides from estimated orbits. + +Future development will extend the IOD capabilities beyond Gauss, notably through the addition of a Väisälä solver [@Williams1997] and a full differential-corrections pipeline with covariance output. Support for hyperbolic trajectories will enable the analysis of interstellar candidates, and a unified API for orbit-based ephemeris prediction will facilitate downstream residual analysis. + +# Acknowledgements + +We acknowledge the maintainers of OrbFit for foundational work in preliminary orbit determination, the JPL Horizons/NAIF teams for ephemerides and kernels, and contributors. Development benefited from the Rust ecosystem and libraries cited. + +# References