layup.orbitfit

Attributes

get_ias15_adaptive_mode
logger
ARCSEC_PER_HOUR_TO_RAD_PER_DAY
US_TO_DAYS
_DEFAULT_DELAY_UNC_DAYS
_DEFAULT_DOPPLER_UNC_HZ
REQUIRED_INPUT_OBSERVATIONS_COLUMN_NAMES
REQUIRED_INPUT_GUESS_COLUMN_NAMES
INPUT_FORMAT_READERS
_NONGRAV_BITS
_AUTO_NONGRAV_LADDER
_AUTO_DEFAULT_THRESHOLDS
_FINGERPRINT_COLUMNS
_PICKER_MIN_R_HELIO_AU
_PICKER_SCREEN_ITER_MAX
_PICKER_FULL_ITER_MAX
_PREFILTER_THRESHOLD_SIGMA
_PICKER_IAS15_ADAPTIVE_MODE
_assist_python_ephem_cache
_MIN_ARC_DAYS

Classes

NongravAutoThresholds

Decision thresholds for adaptive non-grav selection (fit_nongrav="auto",

Functions

_run_fit(assist_ephem, initial_guess, observations, engine)	Dispatch a single LM fit step to the configured engine.
_parse_nongrav(fit_nongrav)	Normalize the fit_nongrav argument to (mask, names).
_gravity_fit_acceptable(csq, ndof[, thresholds])	Whether a gravity-only fit is good enough that no non-grav is warranted.
_nongrav_warranted(csq_gravity, res_ng, names[, ...])	Whether adopting the (converged, well-conditioned) non-grav fit res_ng
_gofr_arg(nongrav_gr)	Normalize a non-grav g(r) argument to the [alpha, nm, nn, nk, r0] list the
_select_nongrav_auto(assist_ephem, res_grav, observations)	Adaptive non-grav selection for fit_nongrav="auto" (issue #357).
_get_result_dtypes(primary_id_column_name[, ...])	Helper function to create the result dtype with the correct primary ID column name.
_fmt_fingerprint_value(v)	Canonical, round-trip-stable string for one observation field value.
_obs_fingerprint(data, column_names)	Deterministic fingerprint of an object's fittable observation set.
_is_radar(d, column_names)	True if a row carries a populated radar observable (delay or Doppler).
_radar_observation(objID, d, epoch_jd, column_names)	Build a radar Observation from a row, converting JPL units to the
_append_observer_acceleration(data, observatory[, dt_sec])	Append barycentric observer acceleration columns (ax, ay, az; au/day^2).
_is_occultation(data)	Check if a given data point is an occultation measurement.
_use_star_astrometry(data)	Use occulting star's astrometry to replace the ra and dec values.
_split_by_index(input_list, indices)	Split an input list into a list of sublists, with break
_time_distance(ic0, ic1, jds)	Find the separation in time between two set of indices
_nearest_chunk(target, index_chunks, jds[, self_match])	Given an index chunk (target) and the set of all the index
_next_nearest(chunk_sequence, index_chunks, jds)	Given a sequence of index chunks and a collection of other index
_iterate_sequence(sequence, other_chunks, jds)	Given a sequence of index chunks and a collection of other index
_build_sequence(jds[, sep_dt])	Given a set of julian dates, in the order of the observations,
create_empty_result(id, dtypes)	Create an empty return object
_carry_forward_result(prior_row, dtypes)	Copy a prior fit result forward under the current output dtype (issue #419).
_partition_unchanged(data, initial_guess, ...)	Split raw observations into (to_fit, carried_forward) by obs fingerprint.
do_gauss_iod(observations, seq)	Backward-compat wrapper for the Gauss IOD.
_get_python_ephem(cache_dir)	Lazy-load and cache the Python-side assist.Ephem for the filter.
_pick_best_root(candidates, min_r_au)	Pick the best converged candidate from a list of LM results.
do_fit(observations, seq, cache_dir[, iod, engine, ...])	Carry out an orbit fit to a list of observations.
do_other_fit(iod)	This is a place holder function for future IOD implememtations that require
_warn_if_short_arc(jds, obj_id)	Emit a helpful warning when a failed fit is likely caused by too short an
_orbitfit(data, cache_dir, primary_id_column_name[, ...])	This function will contain all of the calls to the c++ code that will
orbitfit(data, cache_dir[, initial_guess, ...])	This is the function that you would call interactively. i.e. from a notebook
_observations_for_update(data, cache_dir[, ...])	Augment one object's observations with the observer barycentric state and
_update_mahalanobis(prior_fit, updated_fit)	Mahalanobis distance of the state update in prior standard deviations:
sequential_update(prior, new_data, cache_dir, *[, ...])	Sequential / information-filter update of a prior orbit fit (issue #419).
_obs_row_keys(data, column_names)	Per-observation identity keys over the fit-relevant columns (issue #419).
_append_only_new_obs(current, prior_obs)	Return the rows of current absent from prior_obs if the change is
_group_by_id(data, primary_id_column_name)
_fitresult_to_row(fit, obj_id, obs_hash, nobs_fit, dtypes)	Pack a FitResult (from the sequential update or its refit fallback) into one
incremental_orbitfit(data, cache_dir, prior_catalog, *)	Steady-state incremental fit over a batch of objects (issue #419 capstone).
orbitfit_cli(input, input_file_format, output_file_stem)	This is the function that is called from the command line
_is_valid_data(data)	Check if the input data contains all valid values.
is_numeric(obj)

Module Contents

get_ias15_adaptive_mode

logger

ARCSEC_PER_HOUR_TO_RAD_PER_DAY

US_TO_DAYS = 1.1574074074074074e-11

_DEFAULT_DELAY_UNC_DAYS = 1.1574074074074074e-11

_DEFAULT_DOPPLER_UNC_HZ = 1.0

REQUIRED_INPUT_OBSERVATIONS_COLUMN_NAMES

REQUIRED_INPUT_GUESS_COLUMN_NAMES = ['epochMJD_TDB', 'FORMAT']

INPUT_FORMAT_READERS

_run_fit(assist_ephem, initial_guess, observations, engine, iter_max=100)

Dispatch a single LM fit step to the configured engine.

Centralizing the dispatch here keeps do_fit’s IOD-then-fit pipeline parameterization-agnostic and lets us add new engines (e.g., a future distance-dispatched ‘auto’) with a single edit instead of threading the choice through every call site.

iter_max is the LM iteration budget used by the multi-root picker’s two-tier (cheap-screen then full) passes. The Cartesian engine honors it; the BK-native engine uses its own internal cap (it takes mu for the bound-orbit energy prior rather than an iteration budget), so iter_max is ignored on that path.

_NONGRAV_BITS

_parse_nongrav(fit_nongrav)

Normalize the fit_nongrav argument to (mask, names).

Accepts False/None (no non-grav fit), True (== ["A2"], the common asteroid Yarkovsky case), a string naming params (e.g. "A2", "A1A2A3", "A1,A3"), or an iterable of names. Returns the C++ bitmask (bits 1/2/4 for A1/A2/A3) and the selected names ordered A1, A2, A3.

_AUTO_NONGRAV_LADDER = (('A2',), ('A1', 'A2'), ('A1', 'A2', 'A3'))

class NongravAutoThresholds

Decision thresholds for adaptive non-grav selection (fit_nongrav="auto", issue #357). Pass a customized instance to orbitfit to tune how readily a non-gravitational model is adopted; the defaults reproduce the standard “introduce a non-grav only when gravity is unacceptable and the parameter is statistically warranted” behavior.

Parameters:

• accept_reduced_chi2 (float) – A gravity-only fit whose reduced chi-square is at or below this is kept as-is; non-grav models are tried only above it. Default 1.5.

• delta_chi2_per_param (float) – Minimum chi-square drop required per added non-grav parameter to adopt a model (9.0 ~ 3-sigma). Default 9.0.

• nsigma (float) – Each added non-grav parameter must exceed this many times its 1-sigma uncertainty to be adopted. Default 3.0.

accept_reduced_chi2: float = 1.5

delta_chi2_per_param: float = 9.0

nsigma: float = 3.0

_AUTO_DEFAULT_THRESHOLDS

_gravity_fit_acceptable(csq, ndof, thresholds=_AUTO_DEFAULT_THRESHOLDS)

Whether a gravity-only fit is good enough that no non-grav is warranted.

True when the reduced chi-square is at or below the acceptance threshold (or there are no degrees of freedom to judge it by).

_nongrav_warranted(csq_gravity, res_ng, names, thresholds=_AUTO_DEFAULT_THRESHOLDS)

Whether adopting the (converged, well-conditioned) non-grav fit res_ng for parameters names over the gravity-only fit is statistically warranted: a significant chi-square drop AND every added parameter individually significant.

_gofr_arg(nongrav_gr)

Normalize a non-grav g(r) argument to the [alpha, nm, nn, nk, r0] list the C++ fit expects (empty -> the default inverse-square law).

_select_nongrav_auto(assist_ephem, res_grav, observations, thresholds=_AUTO_DEFAULT_THRESHOLDS, gofr=None)

Adaptive non-grav selection for fit_nongrav="auto" (issue #357).

res_grav is the converged gravity-only fit. Returns the most parsimonious acceptable model: the gravity-only result unless a non-grav model is both well-conditioned (flag 0) and statistically warranted per thresholds. gofr selects the g(r) sublimation law (see orbitfit’s nongrav_gr).

_get_result_dtypes(primary_id_column_name: str, nongrav_names=(), per_arc=False)

Helper function to create the result dtype with the correct primary ID column name.

For each fitted non-gravitational parameter in nongrav_names (a subset of A1/A2/A3), two columns are appended – e.g. a2 and a2_unc (the value and its 1-sigma uncertainty, au/day^2). With no non-grav params the default 6-parameter output schema is unchanged.

When per_arc is set (the two-apparition comet-linkage fit), the non-grav columns above hold the earlier arc’s amplitudes and a second block a2_arc2/a2_arc2_unc is appended for the later arc. Both are opt-in, so the ordinary fit schema is unaffected.

_FINGERPRINT_COLUMNS = ('obsTime', 'ra', 'dec', 'stn', 'raStar', 'decStar', 'rmsRA', 'rmsDec', 'astCat', 'raRate',...

_fmt_fingerprint_value(v)

Canonical, round-trip-stable string for one observation field value.

Uses Python’s shortest round-trip repr for floats (deterministic across runs and platforms) and decodes bytes so a numpy S/O string column hashes identically however it was loaded.

_obs_fingerprint(data, column_names)

Deterministic fingerprint of an object’s fittable observation set.

Returns (nobs, hash16) where nobs is the row count and hash16 is a 16-hex-character SHA-1 digest over the fit-relevant columns (_FINGERPRINT_COLUMNS). The per-row strings are sorted before hashing, so the fingerprint is order-independent: the same physical observations produce the same hash regardless of row order (an LM fit over a set is itself order-independent). Computed on the raw observations before any in-place debiasing so it reflects what was reported, not the current bias model.

_is_radar(d, column_names)

True if a row carries a populated radar observable (delay or Doppler).

Radar rows have no ra/dec; they are dispatched to Observation.from_radar rather than the astrometry/streak factories.

_radar_observation(objID, d, epoch_jd, column_names)

Build a radar Observation from a row, converting JPL units to the fitter’s internal units.

delay (us, round-trip) -> days; Doppler (Hz) -> round-trip range-rate (au/day) via the per-observation transmit frequency freqTx. The barycentric observer position/velocity columns (x,y,z,vx,vy,vz) must already be present (added by orbitfit); the observer acceleration columns (ax,ay,az) drive the two-leg light-time model and default to zero when absent. 1-sigma uncertainties come from rmsDelay/rmsDoppler when present, else the module defaults.

_append_observer_acceleration(data, observatory, dt_sec=2.0)

Append barycentric observer acceleration columns (ax, ay, az; au/day^2).

Finite-differences the barycentric station velocity from obscodes_to_barycentric at the observation epoch +/- dt_sec. Used only by radar observations, whose two-leg light-time model extrapolates the station state back to the signal transmit time. Requires the et column (seconds past J2000, TDB) that orbitfit adds before computing the observer states.

_is_occultation(data)

Check if a given data point is an occultation measurement. An occultation measurement is indicated by the lack of ra and dec values.

Parameters:

data (numpy structured array) – The object data to check for occultation measurements.

Returns:

True if the data point is an occultation measurement, False otherwise.

Return type:

bool

_use_star_astrometry(data)

Use occulting star’s astrometry to replace the ra and dec values.

Notes

The units are a bit odd here. raStar and decStar are in degrees. However, deltaRA and deltaDec are in arcseconds. Thus will either convert to degrees or radians depending on the context.

For more details see the ADES description here: https://github.com/IAU-ADES/ADES-Master/blob/master/ADES_Description.pdf

Parameters:

data (numpy structured array) – The object data to replace the ra and dec values.

Returns:

data – The object data with the ra and dec values replaced by the star’s astrometry.

Return type:

numpy structured array

_split_by_index(input_list, indices)

Split an input list into a list of sublists, with break points at the provided indices.

Parameters:

• input_list (a list of, in this case, integers.)

• indices (a list of indices at which to split the input list.)

Return type:

list of lists.

_time_distance(ic0, ic1, jds)

Find the separation in time between two set of indices (index chunks), given the corresponding set of julian dates.

Parameters:

• ic0 (a list of indices (index chunk 0))

• ic1 (a list of indices (index chunk 1))

• jds (a collection of julian dates.)

Returns:

float

Return type:

the size of the time gap separating the two chunks.

_nearest_chunk(target, index_chunks, jds, self_match=False)

Given an index chunk (target) and the set of all the index chunks, find the index chunk that is nearest in time, that has the smallest separation in time.

Parameters:

• target (a list of indices (index chunk))

• index_chunks (a list of lists of indices)

• jds (a collection of julian dates.)

• self_match (allow the target to match to itself)

Returns:

• list (the nearest index chunk.)

• float (the separation in time between the target and the nearest chunk.)

_next_nearest(chunk_sequence, index_chunks, jds)

Given a sequence of index chunks and a collection of other index chunks, find the index chunk that is nearest in time to the sequence. This will be the next one to include in the fit.

Parameters:

• chunk_sequence (a list of lists of indices)

• index_chunks (a list of lists of indices)

• jds (a collection of julian dates.)

• self_match (allow the target to match to itself)

Returns:

list

Return type:

the nearest index chunk.

_iterate_sequence(sequence, other_chunks, jds)

Given a sequence of index chunks and a collection of other index chunks, iteratively build the sequence.

Parameters:

• sequence (a list of lists of indices (the start of the sequence))

• other_chunks (a list of lists of indices)

• jds (a collection of julian dates.)

Returns:

list of lists

Return type:

the resulting list of lists of indices

_build_sequence(jds, sep_dt=90.0)

Given a set of julian dates, in the order of the observations, split the observations into sets of indices that have gaps of at least sep_dt (days) between them. These are index chunks. Then develop a sequence of the index chunks that will be fit in order. The first chunk will have the longest time span. The next will be the closest in time to the first. The next will be the closest in time to the first two, etc.

Parameters:

• jds (a collection of julian dates in order of the observations (sorted).)

• sep_dt (float the mininum separation between chunks)

Returns:

• list of lists (the resulting sequence of lists of indices that)

• will be fit.

create_empty_result(id, dtypes)

Create an empty return object

Parameters:

• id (str) – The id of the object to provide an empty result for

• dtypes (np.array) – The list of datatypes for the structured array.

Returns:

Empty numpy structured array

Return type:

np.array

_carry_forward_result(prior_row, dtypes)

Copy a prior fit result forward under the current output dtype (issue #419).

Used by the skip-unchanged path: when an object’s observations are unchanged, its stored fit is emitted verbatim rather than recomputed. Fields shared with dtypes are copied by name (so the carried row is identical to the prior fit); any field present in dtypes but absent from the prior (e.g. a non-grav column the prior run did not fit) is left at its type default.

_partition_unchanged(data, initial_guess, primary_id_column_name, fit_nongrav)

Split raw observations into (to_fit, carried_forward) by obs fingerprint.

The steady-state pre-filter for orbitfit(skip_unchanged=True) (issue #419): an object whose fingerprint matches a converged prior fit over the same observation set is carried forward verbatim (cast to the current output dtype); every other object’s rows are returned for fitting. Runs before any ephemeris/observatory setup, so a skipped object costs only a fingerprint hash. The fingerprint is computed on the raw (pre-debias) observations, so it matches the one _orbitfit stores.

do_gauss_iod(observations, seq)

Backward-compat wrapper for the Gauss IOD.

Prefer layup.iod.get_iod("gauss") for new code; this shim exists so callers that imported do_gauss_iod directly continue to work.

_PICKER_MIN_R_HELIO_AU = 0.3

_PICKER_SCREEN_ITER_MAX = 80

_PICKER_FULL_ITER_MAX = 100

_PREFILTER_THRESHOLD_SIGMA = 1000.0

_PICKER_IAS15_ADAPTIVE_MODE = 2

_assist_python_ephem_cache: dict

_get_python_ephem(cache_dir)

Lazy-load and cache the Python-side assist.Ephem for the filter.

_pick_best_root(candidates, min_r_au)

Pick the best converged candidate from a list of LM results.

“Best” means smallest χ² among candidates that

1. report flag == 0 (LM converged), and

2. have heliocentric distance > min_r_au (physical orbit).

Returns None if no candidate satisfies (1); in that case the caller typically retries at a larger LM budget. If (1) is met but (2) isn’t, the smallest-χ² convergent root is still returned (better than nothing).

do_fit(observations, seq, cache_dir, iod='gauss', engine='cartesian', screen_iter_max: int = _PICKER_SCREEN_ITER_MAX, full_iter_max: int = _PICKER_FULL_ITER_MAX, min_r_helio_AU: float = _PICKER_MIN_R_HELIO_AU, prefilter_threshold_sigma: float = _PREFILTER_THRESHOLD_SIGMA, picker_ias15_adaptive_mode: int = _PICKER_IAS15_ADAPTIVE_MODE)

Carry out an orbit fit to a list of observations.

Pipeline:

1. IOD: produce one or more candidate seed orbits via the registered method named by iod (default: “gauss”). The registry lives in layup.iod; register new methods with iod.register_iod(name, callable).

2. Multi-root picker: run LM from every IOD candidate on the primary segment (seq[0]) at a cheap screen_iter_max budget. Pick the smallest-χ² converged candidate with heliocentric distance above min_r_helio_AU. If nothing converges at the cheap budget, retry at full_iter_max. Then refit on the full observation set.

Parameters:

• observations (list) – Time-ordered list of layup Observations.

• seq (list of lists) – Per-segment index lists; seq[0] is the primary segment.

• cache_dir (str) – Directory holding the ASSIST kernels.

• iod (str) – Name of the registered IOD method (default “gauss”), or “auto” to run Gauss and fall back to the BK 5-parameter linear IOD (run_bk_iod) on the primary segment when every Gauss root fails to seed a converged fit.

• engine (str) –

  Which LM fitter to dispatch to. Supported:

  • ’cartesian’ (default): the existing 6D Cartesian-state fit.

  • ’bk_native’: the universal Bernstein-Khushalani fit (run_bk_native_fit), with a fixed bound-orbit energy prior on gdot. Recovers the Cartesian state at the same epoch.

• screen_iter_max (int) – Two-tier LM iteration caps for the multi-root picker.

• full_iter_max (int) – Two-tier LM iteration caps for the multi-root picker.

• min_r_helio_AU (float) – Lower bound on heliocentric distance for accepted IOD roots.

Returns:

Best converged fit (flag == 0) when one exists, else a best-effort or sentinel FitResult with a non-zero flag.

Return type:

FitResult

do_other_fit(iod: str)

This is a place holder function for future IOD implememtations that require more significant data manipulation.

_MIN_ARC_DAYS = 1.0

_warn_if_short_arc(jds, obj_id)

Emit a helpful warning when a failed fit is likely caused by too short an observational arc (less than ~24 hours / a single night).

See issue #312: an orbit fit needs a baseline of more than 24 hours, so when a fit fails on a sub-day arc we tell the user the likely cause rather than leaving them with an opaque failure.

_orbitfit(data, cache_dir: str, primary_id_column_name: str, initial_guess=None, bias_dict: dict = None, sort_array: bool = True, weight_data=False, iod: str = 'gauss', engine: str = 'cartesian', fit_nongrav: bool = False, nongrav_auto_thresholds=None, nongrav_gr=None, per_arc: bool = False, skip_unchanged: bool = False)

This function will contain all of the calls to the c++ code that will calculate an orbit given a set of observations. Note that all observations should correspond to the same object.

This is function that is passed to the parallelizer.

Parameters:

• data (numpy structured array) – The object data to derive an orbit for

• cache_dir (str) – The directory where the required orbital files are stored

• primary_id_column_name (str) – The name of the primary identifier column for the objects.

• initial_guess (numpy structured array) – Optional guess data to use for the orbit fit. Default is None.

• bias_dict (dict) – A dictionary containing bias corrections for different catalogs.

• sort_array (bool) – Whether to sort the observations by obstime before processing. Default is True.

• weight_data (bool or str) – Astrometric weighting. False (default) leaves the built-in default uncertainty. True applies the Veres 2017 model (observation code, date, catalog, program). "supplied" uses the per-observation rmsRA / rmsDec columns (arcseconds) directly – e.g. ADES-reported uncertainties, or an external weighting model such as era-based historical weighting for old comet apparitions (a row with a NaN/nonpositive value falls back to the default).

• iod (str) – The IOD used to generate an initial guess orbit. Supports ‘gauss’ (default) and ‘auto’ (Gauss with BK-IOD fallback). Default is ‘gauss’.

orbitfit(data, cache_dir: str, initial_guess=None, num_workers=1, primary_id_column_name='provID', debias=False, weight_data=False, iod='gauss', engine='cartesian', fit_nongrav=False, nongrav_auto_thresholds=None, nongrav_gr=None, per_arc=False, skip_unchanged=False)

This is the function that you would call interactively. i.e. from a notebook

Parameters:

• data (numpy structured array) – The object data to derive an orbit for

• cache_dir (str) – The directory where the required orbital files are stored

• initial_guess (numpy structured array) – Optional initial guess data to use for the orbit fit. Default is None.

• num_workers (int) – The number of workers to use for parallel processing. Default is 1

• primary_id_column_name (str) – The name of the primary identifier column for the objects. Default is “provID”.

• debias (bool) – Whether to apply debiasing corrections to the observations. Default is False.

• weight_data (bool or str) – Astrometric weighting. False (default) leaves the built-in default uncertainty. True applies the Veres 2017 model (observation code, date, catalog, program). "supplied" uses the per-observation rmsRA / rmsDec columns (arcseconds) directly – e.g. ADES-reported uncertainties, or an external weighting model such as era-based historical weighting for old comet apparitions (a row with a NaN/nonpositive value falls back to the default).

• iod (str) – The IOD used to generate an initial guess orbit. Supports ‘gauss’ (default) and ‘auto’ (Gauss with BK-IOD fallback). Default is ‘gauss’.

• fit_nongrav (bool | str | iterable of str) – Which non-gravitational Marsden parameters to fit after the 6-parameter orbit converges. False (default) fits none; True fits A2 (the transverse Yarkovsky term, the common asteroid case); a string or iterable naming params – e.g. "A2", "A1A2A3", ["A1", "A3"] – selects a subset of A1 (radial), A2 (transverse), A3 (normal). "auto" selects the model adaptively per object (issue #357): the gravity-only fit is kept unless its reduced chi-square is unacceptable, in which case the most parsimonious non-grav model that is well-conditioned and statistically significant is adopted (the A1/A2/A3 columns are all present, with only the adopted params filled and the rest NaN). For each fitted param an a{n} value and a{n}_unc 1-sigma column (au/day^2) are added to the result. Cartesian engine only; params weakly constrained on short arcs are reported as NaN (issue #351).

• nongrav_auto_thresholds (NongravAutoThresholds, optional) – Decision thresholds used when fit_nongrav="auto" – how unacceptable the gravity-only fit must be before a non-grav is tried, and how large the chi-square drop and per-parameter significance must be to adopt one. Default (None) uses the standard thresholds; pass a customized NongravAutoThresholds to tune. Ignored unless fit_nongrav="auto".

• nongrav_gr (sequence of float, optional) – The non-gravitational g(r) sublimation law as [alpha, nm, nn, nk, r0] in ASSIST’s parameterization g(r) = alpha*(r/r0)^-nm*(1+(r/r0)^nn)^-nk. Default (None) is the asteroidal inverse-square law (r/r0)^-2 used by Yarkovsky A2 fits; pass a cometary law (e.g. Marsden water-ice) to fit a comet’s non-gravs. Applies to any non-grav fit (explicit or "auto").

• per_arc (bool, optional) – Fit piecewise-constant per-apparition non-grav amplitudes (comet linkage). The state and g(r) are shared, but observations before the fit epoch (the earlier arc) and after it (the later arc) each get their own [A1,A2,A3]. Requires an explicit fit_nongrav mask and an initial_guess whose epoch sits between the two apparitions. The output adds a{1,2,3}_arc2 columns for the later arc; the base a{1,2,3} columns then hold the earlier arc. Default False.

• skip_unchanged (bool) – Incremental / steady-state mode (issue #419). When True and initial_guess is a prior result catalog carrying the obs_hash fingerprint columns, any object whose observation set is byte-for-byte unchanged since that catalog is carried forward verbatim instead of re-fitting. Objects with changed obs are re-fit, warm-started from the prior state when available. Default False (every object is fit). The output always carries the obs_hash/nobs_fit columns so it can seed the next cycle.

_observations_for_update(data, cache_dir, weight_data=False, bias_dict=None)

Augment one object’s observations with the observer barycentric state and build the C++ Observation list, mirroring orbitfit()’s preprocessing.

Supports optical astrometry and streak (rate) rows – the observation kinds a steady-state catalog update sees. (Radar/occultation are not yet handled by the sequential path; the driver’s full-refit fallback covers them.)

_update_mahalanobis(prior_fit, updated_fit)

Mahalanobis distance of the state update in prior standard deviations: sqrt((x1 - x0)^T P0^-1 (x1 - x0)). This is the dimensionless size of the move the new observations induced, and the natural nonlinearity gate – a large move means the linearization of the summarized old observations (which is anchored at the prior mean x0) is no longer trustworthy.

sequential_update(prior, new_data, cache_dir, *, all_data=None, weight_data=False, debias_data=False, max_update_sigma=4.0, iter_max=100)

Sequential / information-filter update of a prior orbit fit (issue #419).

Refines prior using ONLY new_data (the newly reported observations of one object); the previously-fit observations enter through the prior’s state and 6x6 covariance, which becomes the information matrix Lambda0 = P0^-1 added to the normal equations. To the extent the prior is locally Gaussian this equals a full batch refit over all observations, at the cost of integrating only the new ones – the throughput win for steady-state catalog maintenance.

Parameters:

• prior (FitResult or numpy structured array (one row)) – The prior fit (state, 6x6 covariance, epoch). A result-catalog row is accepted and parsed via parse_fit_result.

• new_data (numpy structured array) – The new observations of this object (optical astrometry and/or streaks).

• cache_dir (str or None) – Kernel/ephemeris cache directory (None -> the default layup os_cache).

• all_data (numpy structured array, optional) – The full observation set (old + new). Required for the nonlinearity fallback: when the update is too large the driver refits over all_data.

• weight_data (bool) – Apply Veres (2017) weighting / MPC debiasing to the new observations, matching the corresponding orbitfit options.

• debias_data (bool) – Apply Veres (2017) weighting / MPC debiasing to the new observations, matching the corresponding orbitfit options.

• max_update_sigma (float) – Nonlinearity gate. If the update moves the state more than this many prior standard deviations (Mahalanobis), fall back to a full refit over all_data when provided, else flag the result (flag=8).

• iter_max (int) – LM iteration cap.

Returns:

The updated fit. method is "sequential_update" for an accepted information-filter update, or "orbit_fit" when the fallback refit ran.

Return type:

FitResult

_obs_row_keys(data, column_names)

Per-observation identity keys over the fit-relevant columns (issue #419).

Each key is the same per-row string that _obs_fingerprint hashes, so a row’s key equals its contribution to the object’s fingerprint. Used to diff a current observation set against the one a prior fit was built from.

_append_only_new_obs(current, prior_obs)

Return the rows of current absent from prior_obs if the change is append-only, else None.

Append-only means every observation the prior was fit from is still present (none removed or modified) – the case the sequential update handles exactly. If any prior observation is gone, the summarised old-obs information no longer matches the current set, so the object must be fully re-fit and this returns None.

_group_by_id(data, primary_id_column_name)

_fitresult_to_row(fit, obj_id, obs_hash, nobs_fit, dtypes)

Pack a FitResult (from the sequential update or its refit fallback) into one result-catalog row carrying the current-obs fingerprint, so the row is identical in shape to orbitfit output and can seed the next cycle.

incremental_orbitfit(data, cache_dir, prior_catalog, *, prior_obs=None, primary_id_column_name='provID', weight_data=False, debias=False, max_update_sigma=4.0, iod='gauss', engine='cartesian', num_workers=1)

Steady-state incremental fit over a batch of objects (issue #419 capstone).

Ties the three levers into one operational maintenance pass. For each object in data (the current observations), routes:

• skip – the observation set is unchanged since prior_catalog (matching fingerprint): carry the prior fit forward verbatim, no fit.

• sequential update – observations were only appended (prior_obs given and every prior observation is still present): update the prior with the new observations only, via sequential_update() (integrating just the new obs). Its nonlinearity gate falls back to a full refit when the update is too large.

• full refit – observations were removed or changed, or no per-object prior_obs is available: refit over all current obs, warm-started from the prior state when there is one, cold (IOD) when the object is new.

Parameters:

• data (numpy structured array) – Current observations for all objects (grouped by primary_id_column_name).

• cache_dir (str or None) – Kernel/ephemeris cache directory (None -> the default layup os_cache).

• prior_catalog (numpy structured array or None) – Prior fit results carrying state/cov/epoch and the obs_hash/nobs_fit fingerprint columns (i.e. produced by orbitfit/this driver). None -> every object is cold-fit.

• prior_obs (numpy structured array, optional) – The observations the prior catalog was fit from, grouped by id. Enables the sequential route (needs the per-object obs to diff). Without it, changed objects are fully refit.

• weight_data (bool) – Veres (2017) weighting / MPC debiasing, as in orbitfit.

• debias (bool) – Veres (2017) weighting / MPC debiasing, as in orbitfit.

• max_update_sigma (float) – Nonlinearity gate passed to sequential_update().

Returns:

The updated result catalog (one row per object, same schema as orbitfit output) and a routing tally {"skip", "sequential", "sequential_fallback", "full", "cold"}.

Return type:

(numpy structured array, dict)

orbitfit_cli(input: str, input_file_format: Literal['MPC80col', 'ADES_csv', 'ADES_psv', 'ADES_xml', 'ADES_hdf5'], output_file_stem: str, output_file_format: Literal['csv', 'hdf5'] = 'csv', chunk_size: int = 10000, num_workers: int = -1, cli_args: argparse.Namespace | None = None)

This is the function that is called from the command line

Parameters:

• input (str) – Path to the input data file.

• input_file_format (Literal[MPC80col, ADES_csv, ADES_psv, ADES_xml, ADES_hdf5]) – The format of the input data file.

• output_file_stem (str) – The stem of the output file.

• output_file_format (Literal[csv, hdf5] optional (default="csv")) – The format of the output file. Must be one of: “csv”, “hdf5”

• num_workers (int, optional (default=-1)) – The number of workers to use for parallel processing of the individual chunk. If -1, the number of workers will be set to the number of CPUs on the system. The default is 1 worker.

• cli_args (argparse.Namespace, optional (default=None)) – The argparse object that was created when running from the CLI.

_is_valid_data(data)

Check if the input data contains all valid values.

Parameters:

data (numpy structured array) – The object data to validate.

Returns:

True if the data is valid, False otherwise.

Return type:

bool

is_numeric(obj)