[wei-lingfeng]

Summary

This merges mm_rework_lingfeng into dev. It supersedes #25 (which compared this branch against mm_rework) — everything from that comparison is folded into the sections below, organized by what actually changed and why, alongside a substantial follow-on pass of correctness, performance, memory, and documentation work done since.

API changes from mm_rework

mm_rework already has a motion_model.py and a Fixed/Linear/Acceleration/Parallax system of its own, so unlike the rest of this description, the API differences below are specifically against that branch's design (not dev, which predates motion models entirely and has neither side of any of this).

• default_motion_model (a single model name) → motion_models (a list): mm_rework's MosaicSelfRef.__init__ takes default_motion_model='Fixed', and fit_velocities separately takes its own default_motion_model='Linear' — one model name, used as a single global fallback. Both are replaced by one motion_models parameter (a list, e.g. motion_models=['Linear']): each star is now fit with the most complex model its own number of valid epochs actually supports, rather than every star sharing one default.
• motion_model_dict → fixed_params_dict, and t0 stops being a special case: mm_rework's run_fit(self, t, x, y, xe, ye, t0, ...) takes t0 as its own required positional argument everywhere, with a separate, loosely-typed motion_model_dict={} alongside it for anything else. Both are unified into one fixed_params_dict (holding t0 for Linear/Acceleration; t0, ra, dec, pa, obsLocation for Parallax), so every model's non-fit parameters flow through the same mechanism instead of t0 being threaded separately from everything else.
• New Empty motion model: mm_rework only has Fixed/Linear/Acceleration/Parallax. Empty (a trivial placeholder for a star with no usable data — always fill_value/inf) is new.
• Method rename: StarTable.fit_velocities → StarTable.fit_motion_models, reflecting the move from "always fits a velocity" to "fits whichever model applies." fit_velocities_all_detected, a separate helper for the "every star has every epoch" case, has no equivalent — the general per-star model selection already covers it.
• fit_velocities parameters removed or replaced: mask_val → mask_value; fixed_t0=False (a special-cased boolean/array just for t0) is gone now that t0 is an ordinary fixed_params_dict entry; show_progress/reassign_motion_model are gone (reassignment now happens automatically every call from each star's current epoch count — see keep_existing, new, to control whether stars not being refit keep their existing values instead of being overwritten with nan, described further down). New, with no mm_rework equivalent at all: keep_existing, seed, method (scipy method name), fill_value, art_star, and the multiprocessing controls processes/chunksize/mp_star_threshold.
• MotionModel subclass internals renamed (relevant only if calling code subclasses or touches MotionModel directly, rather than going through StarTable/align): n_pts_req (hardcoded per class, e.g. Parallax.n_pts_req = 4) → n_params (derived as ceil(n_fit_params / 2), giving Parallax.n_params = 3 — the same fix as the n_pts_req point already covered above); fitter_param_names → fit_param_names; optional_param_names (a bare list of names) → optional_fixed_params (a dict of name → actual default value); get_pos_at_time/get_batch_pos_at_time (separate per-star and per-batch methods) unified into one model() method that vectorizes automatically based on input shape; motion_model.get_weights → calc_sigma (see "Uncertainty model" above); new run_fit_batch on Empty/Fixed/Linear, with no mm_rework equivalent (see "Vectorized, multiprocessing-aware motion-model fitting" below).
• Also new since mm_rework, covered in their own sections below: inherit_n_detect, match_workers, mp_star_threshold.

Uncertainty model: error propagation, not weighted variance

The codebase's biggest conceptual change: uncertainty on an averaged/fit quantity is now computed by analytic error propagation (the classic "error on the mean," sigma = sqrt(1 / sum(1/sigma_i^2))) everywhere, instead of the sqrt of a weighted sample variance.

• Why: sample variance is ill-defined at 1 point (forcing a special case) and deceptively well-behaved at 2–3 points — a line fits 2 points exactly, a plane-ish model fits 3, so the "spread of residuals" collapses toward zero regardless of how uncertain the underlying measurements actually were. That's not a real uncertainty estimate, and patching it back up would mean flystar imposing its own peculiar rescaling assumptions; that's a decision for calling code, not this package. Error propagation has neither problem, and is what scipy.curve_fit's covariance matrix already computes internally for models that use it — so this also makes every motion model's reported uncertainty consistent with every other one, instead of Fixed computing something structurally different from Linear/Acceleration/Parallax.
• combine_lists switched to error propagation (also removing its now-unnecessary >1-epoch requirement for weighted averaging), and Fixed.run_fit's uncertainty was brought in line with the scipy-covariance convention every other model already used.
• motion_model.get_weights's x_wt = 1/xe**2 → sigma = 1/x_wt**0.5 round-trip accumulated floating-point error across repeated align() iterations; replaced with calc_sigma, which returns xe or sqrt(xe) directly usable by scipy, computing an inverse-variance weight only where one's actually still needed.
• Parallax.run_fit had two real bugs from not using that same path: sigma was never square-rooted before being passed to curve_fit, and absolute_sigma was silently ignored. Both fixed by routing it through calc_sigma + absolute_sigma like everything else. (Worked example proving curve_fit(..., absolute_sigma=True) and analytic error propagation agree bit-for-bit is preserved from Update mm_rework to match mm_rework_lingfeng for comparison #25, at the bottom of this description.)
• Related default flip: MosaicSelfRef/MosaicToRef's absolute_sigma default changed from False to True — required for scipy's covariance matrix to actually mean "error propagation" rather than a reduced-chi2-rescaled quantity.

Uncertainty columns: inf, never nan, for "no information"

A second, related convention, enforced consistently for the first time: a star's uncertainty is inf when there's no real information to compute one from — never nan. 1/inf**2 == 0, so an inf uncertainty correctly contributes zero weight anywhere it's combined with other stars' measurements; nan instead poisons every downstream computation it touches (x + nan == nan) and vanishes silently from filters, since nan == nan is False and nan < anything is False. Several real bugs were all instances of this convention being violated in different places:

• xe/ye/me defaulted to nan instead of inf for invalid entries in the newly-vectorized Fixed/Empty batch paths — corrected to inf.
• combine_lists' weighted "no usable uncertainty anywhere" fallback used to patch a fake weight into the raw error array, fit through that, and then had to remember to force the reported error back to inf afterward — fragile, and exactly the kind of place a fabricated finite value could silently leak through instead. Rewritten so inf falls out naturally: the weight sum is built only from real, known uncertainties, so std = sqrt(1/wgt_sum) is inf via ordinary 1/0 whenever none exist, with no override to remember.
• absolute_sigma=False's chi2 rescaling could turn a singular fit's correct inf error into nan — inf * sqrt(nan) == nan. Fixed by reapplying the singular/insufficient-data override unconditionally as the fit's final step, so nothing downstream can rescale it away.
• Two related, genuine NaN-propagation bugs (not the inf convention itself, but bugs the convention's violation exposed on real, imperfect data): a non-finite m0 could never be excluded by a magnitude-range filter (NaN comparisons are always False), letting bad reference stars corrupt the magnitude-zeropoint transform; and fit_motion_models' default t0 could silently come out NaN for some multi-epoch stars from mixing a masked weights array with plain (not masked) averaging. Both fixed.

Vectorized, multiprocessing-aware motion-model fitting

• fit_motion_models now fits a whole group of stars sharing a motion model in one vectorized pass instead of one star at a time, for every model that supports it: Fixed and Empty (trivial/closed-form) and Linear (closed-form 2×2 solve). On a real dataset, processes=3 dropped from 7.07s to 2.68s (matching serial's 2.66s) once all three were batched — a multiprocessing pool no longer needs to spin up at all just to run an O(1) fill for Empty stars.
• A numpy.ma indexing bottleneck (tens of millions of slow per-element numpy.ma.core.__getitem__ calls) made fit_motion_models ~38% faster to fix; one multiprocessing pool is now reused across the whole call instead of respawned per motion-model group.
• match() gained a workers parameter (default 1) threading scipy.spatial.KDTree.query_ball_point; align.py exposes it as match_workers (default 1 — production runs typically share a machine, so the speedup is opt-in). Verified this doesn't affect the delicate dm_min == dr_min tie-break, which depends on within-list neighbor order: workers=1 vs workers=-1 give identical, order-preserved neighbor lists, checked against dense/duplicate-point edge cases and a full end-to-end run with all 37 ref_table output columns byte-identical.
• Added a configurable mp_star_threshold (default 100,000): a multiprocessing pool for motion-model fitting (Acceleration/Parallax, or any model with bootstrap > 0 — the models with no vectorized path) is only spun up when the number of stars actually needing it meets this threshold, even if processes > 1 was requested. Measured break-even for that fixed pool-spawn/IPC overhead was between 20,000 and 100,000 stars on a 10-core machine.

Memory footprint at mosaic scale

• StarTable construction gained an opt-in copy=False and now builds all columns in one constructor call instead of many add_column() calls (~1.2s / +4.6GB → ~0.001s / ~0GB for ~29 columns at 1.4M rows).
• Adding new-star rows built a whole parallel table and vstack-ed it onto the growing ref_table, transiently holding old + new + concatenated data for every column at once — roughly doubling peak memory on every "add new stars" step. Now concatenates columns directly and drops old references immediately.
• fit_motion_models did an O(N_stars) data-prep step regardless of how few stars actually needed refitting; now sliced down to just the selected rows first. A redundant double-copy in combine_lists/fit_motion_models's array prep (an arange-based fancy-index copy stacked on top of an already-copying masked_invalid/deepcopy call) was removed. A per-star fixed-params dict is now built lazily, only for stars whose motion model actually needs it — skipped for ~84% of stars (those already handled by the vectorized batch path, above) in one benchmark.

Other correctness & behavior changes

• fit_star_idxs no longer wipes unrelated stars: refitting only a subset of stars used to overwrite every other star's existing result with nan; now keeps existing values by default and only touches the requested subset.
• Parallax's minimum-epoch requirement was n_pts_req=4; 3 (x, y) pairs are enough to solve its 5-parameter model, so lowered to 3 — a caller wanting to require 4+ can still do so via the motion_model_input column.
• inherit_n_detect: new MosaicSelfRef/MosaicToRef parameter (default True) so a star's n_detect reflects the total number of raw detections it represents across nested alignment layers, not just 1 per input starlist.
• False-positive StarList warning: every pickle round-trip (astropy passes columns positionally to __init__) incorrectly warned about missing required arguments.
• Stale motion-model params: a star that dropped to a simpler motion model (e.g. Linear → Fixed, since epoch counts aren't monotonic across align() iterations) kept stale params (vx/vy) from its old model indefinitely; now reset to fill_value/inf right after reclassification.
• Minor: parallax.parallax_in_direction's Time(mjd + 2400000.5, format='jd', ...) simplified to Time(mjd, format='mjd', ...); a misleading tqdm progress bar over a handful-of-motion-model-types loop (never per-star) removed.

Documentation

• Added .readthedocs.yaml. Getting a working build also required fixing three real, pre-existing bugs unrelated to Read the Docs specifically — they'd have broken a fresh pip install . for anyone on a current Python: pyproject.toml's [build-system] pinned cython==0.29.14 (leftover astropy-template boilerplate for a package with zero .pyx/C extensions), which imports the stdlib cgi module removed in Python 3.13, breaking install outright; the declared dependencies were missing scipy, matplotlib, tqdm, joblib, and pandas, all imported unconditionally at module level; and setup.cfg's github_project was still the template default (astropy/astropy), which conf.py uses unconditionally to build doc issue-links. Verified via two independent fresh-venv builds (editable and the exact non-editable pip install .[docs] Read the Docs runs) — sphinx-build succeeds, 92 warnings, all pre-existing docstring formatting nits unrelated to this change.
• Docs are currently configured to build from wei-lingfeng/flystar rather than this repo directly: importing MovingUniverseLab/flystar into Read the Docs needs repo-admin access (to add its webhook) not currently available on this account, and the flystar project name is already taken there by an existing (currently stale) project for this same package — re-pointing that one at this branch, rather than creating a competing project, needs maintainer access to it instead. The fork tracks this branch 1:1 in the meantime, so content is identical either way.

Validation

Each change was validated individually — synthetic fuzz tests against reference implementations (thousands to tens of thousands of cases per change), byte-for-bit output comparison on real multi-thousand and multi-million-star datasets, and dedicated new unit tests where behavior changed. Full test suite passes except two pre-existing, environment/data-dependent failures unrelated to any change here: test_masked_cols (missing test fixture file) and test_generic_match (test data contains no finite values). (The test_MosaicSelfRef_vel failure present at one point in mm_rework — see the magnitude-weighting fix above — is resolved on this branch.)

Worked example: scipy.curve_fit vs. analytic error propagation (referenced above)

import numpy as np
from scipy.optimize import curve_fit

# 1. Define dummy data and explicit measurement errors
y_data = np.array([10.2, 9.8, 10.5])
sigma_y = np.array([0.1, 0.4, 0.2])   # Different errors for each point
x_dummy = np.zeros_like(y_data)       # x is required but ignored by a constant function

# --- Method A: Analytical Error Propagation ---
weights = 1.0 / (sigma_y ** 2)
c_prop = np.sum(y_data * weights) / np.sum(weights)
sigma_c_prop = np.sqrt(1.0 / np.sum(weights))

# --- Method B: SciPy Curve Fit ---
def constant_model(x, c):
    return c

# absolute_sigma=True is MANDATORY to match pure error propagation
popt, pcov = curve_fit(constant_model, x_dummy, y_data, sigma=sigma_y, absolute_sigma=True)

c_fit = popt[0]
sigma_c_fit = np.sqrt(pcov[0, 0]) # Square root of the variance element

# --- Display Results ---
print(f"Propagation:  Value = {c_prop:.5f}, Error = {sigma_c_prop:.5f}")
print(f"Curve Fit:    Value = {c_fit:.5f}, Error = {sigma_c_fit:.5f}")

Output:

Propagation:  Value = 10.23810, Error = 0.08729
Curve Fit:    Value = 10.23810, Error = 0.08729