Metadata-Version: 2.4
Name: kinextract
Version: 1.0.0
Summary: Non-parametric binned-LOSVD spectral fitter for galaxy kinematics
Author-email: "Thomas K. Waters" <waterstk@umich.edu>
License: MIT License
        
        Copyright (c) 2026 Thomas K. Waters
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
        
Project-URL: Homepage, https://github.com/thomas-k-waters/kinextract
Project-URL: Bug Tracker, https://github.com/thomas-k-waters/kinextract/issues
Project-URL: Documentation, https://kinextract.readthedocs.io
Keywords: astronomy,spectroscopy,kinematics,LOSVD,galaxies
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Astronomy
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Requires-Dist: matplotlib>=3.7
Requires-Dist: jax[cpu]>=0.4
Requires-Dist: threadpoolctl>=3.0
Requires-Dist: astropy>=5.0
Provides-Extra: eso
Requires-Dist: astroquery>=0.4.6; extra == "eso"
Provides-Extra: nuts
Requires-Dist: numpyro>=0.13; extra == "nuts"
Provides-Extra: docs
Requires-Dist: sphinx>=7.0; extra == "docs"
Requires-Dist: sphinx-rtd-theme>=2.0; extra == "docs"
Requires-Dist: myst-parser>=2.0; extra == "docs"
Requires-Dist: sphinx-copybutton>=0.5; extra == "docs"
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: pytest-cov>=4.0; extra == "dev"
Requires-Dist: hypothesis>=6.0; extra == "dev"
Requires-Dist: ruff>=0.4; extra == "dev"
Requires-Dist: black>=24.0; extra == "dev"
Dynamic: license-file

# kinextract

A non-parametric LOSVD (line-of-sight velocity distribution) extractor
for galaxy spectroscopy.

`kinextract` recovers the LOSVD directly from an absorption-line
galaxy spectrum as a smoothly-positive, low-dimensional log-density
spline -- without assuming a Gauss-Hermite (or any other fixed)
functional form for the underlying kinematics. Positivity is structural
(built into the representation via an `exp` link), not an afterthought
enforced after the fact, and the number of free parameters is chosen to
match how much velocity information a spectrum actually contains, not
fixed at an arbitrarily fine grid.

## Status

**v1.0.0** -- a complete, from-scratch rebuild of the fitting engine,
uncertainty quantification, and real-data pipeline, independently
validated end to end before being merged back in (27-stage validation
campaign, `dev_notes/VALIDATION_CAMPAIGN.md`). The default uncertainty
engine (Laplace + AIC-weighted) was validated via a production-scale
coverage sweep across sigma 50-600 km/s, matching or beating parametric
bootstrap's own coverage in every regime tested. Marked `Beta` (not
`Stable`) because it's still actively developed and because several
specific, real limitations are tracked openly below rather than hidden --
not because the validated parts are shaky.

The core fitting, uncertainty (Laplace and bootstrap), multi-template,
cleaning, and continuum machinery is implemented and validated on
synthetic data, and has also been run end-to-end on real VLT/MUSE and
HST/STIS spectra (see `examples/notebooks/05_noise_correlation.ipynb`,
`06_other_instruments.ipynb`, and `07_full_galaxy_analysis.ipynb` for a
full-galaxy analysis). **Real-data-specific finding worth knowing before
trusting a real fit**: real chi2_red on real IFU data commonly runs well
above 1 (systematics-dominated residuals -- template mismatch, unmodeled
continuum shape, correlated per-pixel noise -- not just photon noise),
which affects how much regularization gets selected and how well a
sigma-clipping/outlier-rejection scheme behaves; see
`05_noise_correlation.ipynb` and `07_full_galaxy_analysis.ipynb`'s own
Limitations section for the current, honest picture. Treat every
synthetic validation result in this package as solid, and every real-data
number as requiring the same care those notebooks show.

## Example notebooks

`examples/notebooks/` walks from simple synthetic data to real data, in
order:

- **`01_basic_mock_fit.ipynb`** -- foundations: the spectral model
  equation, the LOSVD's B-spline/`exp`-link parameterization, and the
  core LOSVD + template-mixture recovery machinery on a synthetic mock
  built from real E-MILES templates *you choose* (edit `TEMPLATE_RECIPE`:
  any number of templates, any ages/metallicities/weights -- no separate
  single-template-vs-mixture code path), no external files needed. Stops
  at the point estimate.
- **`02_realistic_mock_fit.ipynb`** -- picks up where 01 left off: a
  genuinely degenerate 4-template mixture, `penalty_order`/`lam_wing`
  regularization trade-offs, a Cramer-Rao-bound noise-floor analysis for
  `h4`, this package's two uncertainty engines (Laplace, the default, and
  bootstrap as a cross-check), and continuum-cofitting -- including why
  a recovered continuum's *shape* is trustworthy but its absolute *level*
  is not identifiable (a real, root-caused property of the model, not a
  bug).
- **`03_template_mismatch.ipynb`** -- fits the *same* mock spectrum with
  a template library that actually contains the true stellar population,
  and with one that systematically doesn't, and compares the recovered
  LOSVD directly -- motivates `04_template_calibration.ipynb`.
- **`04_template_calibration.ipynb`** -- run this first for any real MUSE
  target. Gets a MUSE datacube (a real ESO archive download, a synthetic
  Sersic-profile mock with known ground truth, or shipped pre-computed
  products), bins it with the same PA-aware radial x angular scheme
  `07_full_galaxy_analysis.ipynb` uses for its real kinematic maps, fits
  every bin against the full 636-template E-MILES grid, and iteratively
  reduces to a representative template subset (converges when the
  *union* of surviving templates across bins stabilizes) -- saved for
  reuse in `07_full_galaxy_analysis.ipynb`.
- **`05_noise_correlation.ipynb`** -- measures real MUSE pixel-noise
  correlation (AR(1) lag-1 correlation length) on a real spectrum,
  connects it to the Bacon et al. (2017) MUSE variance-underestimate
  factor, shows correlated-noise whitening (`FitConfig.rho`,
  `measure_and_whiten`) applied to a real fit and its uncertainty band,
  and confirms the whole mechanism on a synthetic mock with known,
  injected correlated noise.
- **`06_other_instruments.ipynb`** -- adapting this package to a
  different instrument (HST/STIS): follows the same 01 -> 02 -> real-fit
  progression, adapted to STIS -- a basic mock and a continuum-cofitting
  mock built entirely from real STIS-observed stellar templates (no
  LSF-matching step needed -- templates and data share the same
  instrument), then a real HST/STIS NGC 5102 spectrum fit end-to-end --
  demonstrates this package's workflow is not MUSE-specific and shows how
  to adapt it to your own instrument.
- **`07_full_galaxy_analysis.ipynb`** -- the capstone: a full real-galaxy
  analysis start to finish -- PA-aware dartboard binning, per-bin
  masking/cleaning/fitting/Laplace uncertainty in parallel, moment maps,
  and radial kinematic profiles along the major axis.

## Features

- **LOSVD recovery** (`kinextract.losvd`, `kinextract.config`,
  `kinextract.multi_template`): one fitting engine (`fit_losvd`)
  recovers the LOSVD -- and, for more than one template, the mixture
  weights jointly via alternating minimization (exact NNLS for the
  weights, the same multi-start MAP fit for the spline coefficients) --
  for any `n_templates >= 1`. See `losvd.py`'s and `multi_template.py`'s
  module docstrings for the method, and
  `examples/notebooks/01_basic_mock_fit.ipynb` for a worked example with
  a configurable template recipe (any number of templates, one mechanism),
  including a template-conditioning diagnostic. Includes an opt-in wing ramp on the
  roughness penalty (`FitConfig.lam_wing`) -- see Known limitations.
- **Uncertainty quantification** (`kinextract.laplace_uncertainty`,
  `kinextract.uncertainty`): two engines. `laplace_losvd_uncertainty`
  (and its per-bin-parallel counterpart, `laplace_losvd_uncertainty_
  parallel`) is the default -- a Laplace (curvature) approximation to the
  posterior at each candidate regularization strength, combined via
  AIC-weighted posterior-predictive sampling; measured directly to give
  substantially better h4 coverage than bootstrap (39% vs. 4% against a
  68% target, production-scale coverage sweep at sigma 50-150 km/s) at a
  fraction of the runtime, though still below the target especially at
  low sigma -- see `dev_notes/uncertainty/FINDINGS.md`, which also
  breaks the residual gap down into a genuine point-estimate-bias
  component and a specific, identified sampling artifact. A separate
  sweep extending through sigma 250-600 km/s (real black-hole-mass
  target regime, e.g. NGC 4751) confirms the same pattern holds well
  beyond the original validation range: sigma coverage stays on-target
  for both engines, and Laplace's h4 advantage is if anything more
  decisive there (50-83% vs. 0-53%). `bootstrap_losvd`/`bootstrap_losvd_
  parallel` (parametric bootstrap, any `n_templates >= 1`, including
  mixture-weight uncertainty) remain available as a cross-check -- real
  ensemble properties (empirical skew/multimodality, basin-selection
  risk) Laplace structurally cannot give. An optional NUTS (NumPyro)
  reference mode (`kinextract.nuts_uncertainty`, the `nuts` extra) is
  available for spot-checking the Laplace approximation itself.
  `bootstrap_losvd` also supports an opt-in `calibrate_noise_scale`/
  `noise_scale` mechanism (not yet available on the Laplace engine) for
  when the assumed per-pixel noise itself is suspected wrong (validated
  on synthetic uniform noise mis-specification; **not yet safe on real
  data**, see Known limitations). See `laplace_uncertainty.py`'s and
  `uncertainty.py`'s own module docstrings for the full, nuanced coverage
  picture.
- **Continuum co-fitting** (`kinextract.continuum`): a shared,
  structurally-positive multiplicative continuum
  (`C(lambda) = exp(spline(lambda))`), jointly recovered with the LOSVD
  and template weights via 3-block alternating minimization. Cleanly
  validated on synthetic data, including the key "does the continuum
  contaminate the weight recovery" adversarial check -- see
  `continuum.py`'s own docstring, and `02_realistic_mock_fit.ipynb`/
  `06_other_instruments.ipynb` for worked demonstrations. The recovered
  continuum's *shape* is reliable; its absolute *mean* is not
  identifiable from the model alone (a real, root-caused mathematical
  property -- `C`'s mean and the template weights' overall scale are the
  same degree of freedom -- not a bug to fix, see `dev_notes/continuum/
  FINDINGS.md`'s 2026-08-01 entry). `penalty_order=3` with a fixed `lam`
  in the 4096-32768 range is the validated default whenever
  `continuum_knots` is set. **Tried on real data and found not to
  help** (see Status above) -- kept because the mechanism itself is
  solid, not because it's currently recommended for this dataset.
- **Outlier rejection** (`kinextract.forward_model.
  mask_wavelength_windows`, `kinextract.cleaning`): fixed-wavelength
  masking for known bad-pixel ranges, and iterative sigma-clipping
  (`fit_losvd_clean`/`fit_losvd_clean_parallel`) for contamination whose
  location isn't known in advance. The clip threshold can be measured
  against a robust (MAD-based) estimate of the residual's own scale
  (`robust_scale=True`) rather than the raw input noise vector -- matters
  on real data, where the input noise vector commonly underestimates the
  true residual scatter (see the chi2_red note above); confirmed to
  behave far better than the raw-threshold version on real MUSE bins that
  otherwise clip runaway. See `07_full_galaxy_analysis.ipynb` for this
  used as the default per-bin cleaning step, and its own Limitations
  section for what it still doesn't handle (continuum-*shape* mismatch,
  in particular).
- **LOSVD moment estimation** (`kinextract.gh_moments`): summarizes
  a recovered LOSVD as (V, sigma, h3, h4) via a direct-moment estimator
  that beats a full Gauss-Hermite refit on this package's own validation
  grid; the joint Gauss-Hermite fit is kept as a general-purpose
  alternative and for comparison against the wider dynamics literature.
- **Noise-correlation whitening** (`kinextract.noise_model`): AR(1)
  pixel-covariance whitening (`whiten`, `FitConfig.rho`), a lag-1 ACF
  correlation-length estimator (`estimate_correlation_length`), and a
  measure-and-whiten convenience function (`measure_and_whiten`) that
  does both in one call -- see `05_noise_correlation.ipynb` for the real
  measurement + consequence walkthrough this exists to support.

Not yet implemented: the template-continuum-trust hierarchy (library-
keyed presets for how much of a raw template's own continuum to trust --
see `continuum.py`'s docstring), wavelength-dependent LSF matching, real
MUSE/STIS pre-normalized data loading, an end-to-end spectrum-file I/O
and orchestration layer.

## Known limitations

- **`penalty_order=3` (third-difference roughness penalty) is the
  default**, replacing an earlier second-difference penalty. Its null
  space is any Gaussian (not a straight line in log-density), so heavy
  smoothing shrinks toward the best-fitting Gaussian rather than a
  slowly-decaying tail -- eliminates the wing-amplitude sigma inflation
  the earlier penalty suffered at low S/N. Real, structural cost in
  exchange: genuine non-Gaussian LOSVDs (nonzero h3/h4) at moderate
  sigma/S-N get their h3/h4 shrunk toward zero (`penalty_order=2` keeps a
  modest edge there instead -- see `02_realistic_mock_fit.ipynb`'s own
  regularization section). Neither is strictly better; which one you want
  depends on whether you're measuring moments or full-curve shape.
- **JAX runs in float64 by default** (opt out with `KINEXTRACT_NO_X64=1`).
  This is a correctness requirement, not a preference: under float32,
  L-BFGS-B can stall far from the optimum in the stiff heavy-penalty
  regime the third-difference penalty operates in on wide velocity grids,
  reporting `success=True` regardless.
- **Real chi2_red on real IFU data commonly runs well above 1** (double
  digits, not ~1) -- systematics-dominated residuals (template mismatch,
  unmodeled continuum shape, correlated per-pixel noise), not just photon
  noise. This affects more than the raw fit: it means a nominal "N sigma"
  threshold in any residual-based scheme (adaptive regularization
  selection, iterative cleaning) is measured against the wrong scale
  unless corrected for -- see the cleaning note below and
  `05_noise_correlation.ipynb`/`07_full_galaxy_analysis.ipynb` for the
  current, honest picture on real data.
- **Adaptive regularization selection (`lam_selection="aic_local"`, the
  default) can select heavy regularization even in the package's own
  validated velocity-grid convention, at low-to-moderate sigma.** Found
  while investigating a related warning in `FitConfig`: at sigma below
  roughly 30-50 km/s, `aic_local` can pick the top of its own candidate
  grid, producing real, sometimes severe sigma bias. Not yet root-caused
  or fixed -- flagged here rather than silently working around it.
- **Iterative sigma-clipping needs its threshold measured against the
  residual's own scale on real data, not the raw input noise vector.**
  `fit_losvd_clean`'s default (`robust_scale=False`) compares the residual
  directly to the input per-pixel noise -- correct when that noise vector
  is well-calibrated (chi2_red ~ 1), but on real data where chi2_red
  commonly runs well into the double digits (see above), a nominal
  "3 sigma" cut is actually far less than 1 true sigma, and can clip the
  majority of a spectrum on the first pass, collapsing the fit to an
  unconstrained, biased baseline. `robust_scale=True` instead measures the
  clip threshold against a MAD-based robust estimate of the *current*
  residual's own scale (recomputed each iteration) -- confirmed
  substantially better-behaved on real bins that broke the default, and
  used as the default cleaning step in `07_full_galaxy_analysis.ipynb`.
  It still needs real per-target tuning (too aggressive a threshold risks
  clipping into genuine absorption features), and has no mechanism for a
  continuum-*shape* mismatch (as opposed to a localized outlier) -- see
  `cleaning.py`'s own docstring and `07_full_galaxy_analysis.ipynb`'s
  Limitations section.
- **The regularization-bias correction (`bias_correction.py`,
  `apply_bias_correction`) is not yet trustworthy enough to default on,
  and is not yet available for the Laplace uncertainty engine at all**
  (calibrated only against `bootstrap_losvd`'s own point estimate). A
  4-gate acceptance test (per-cell bias, held-out coverage, worst-cell
  regression, regression tests) finds sigma coverage passes but h4
  coverage does not -- h4's point estimate is severely shrunk under
  regularization and the reported uncertainty never widens to reflect
  that. `apply_bias_correction` stays `False` by default until this is
  resolved; see `dev_notes/calibration/CALIBRATION_TRUST_STRATEGY.md` and
  `dev_notes/uncertainty/FINDINGS.md`.
- **h3 recovery is severely suppressed at low-to-moderate sigma.** For
  genuinely nonzero true h3, recovered h3 stays close to 0 (essentially
  fully suppressed, not just biased) across a wide low-sigma range, only
  starting to recover at higher sigma -- a calibration table
  characterizing this (`regularization_bias_emiles_h3wing.txt`) exists,
  but the shipped bias-correction machinery does not yet consume it (its
  file format doesn't carry `h3_mean`/`h3_std` columns yet -- see
  `dev_notes/calibration/CALIBRATION_TRUST_STRATEGY.md`'s 2026-08-01
  entry). No correction attempted yet.
- **Multi-template weight recovery is sensitive to template-library
  conditioning** -- a near-collinear template library (differing mainly
  in overall depth scale, not shape) gives genuinely unstable weight
  recovery across noise realizations, diagnosable via the condition
  number of the per-step NNLS design matrix. See
  `examples/notebooks/01_basic_mock_fit.ipynb` and
  `02_realistic_mock_fit.ipynb` (a real, physically-motivated near-
  collinear example), and `03_template_mismatch.ipynb` for the related
  question of what happens when the library is wrong, not just
  degenerate. Degrades gracefully (sparser/noisier weight recovery, never
  a crash or NaN); mitigated but not eliminated by
  template-conditioning (`04_template_calibration.ipynb`'s iterative
  reduction).
- **No continuum model by default.** The core fit expects
  continuum-normalized input; continuum-cofitting (`kinextract.continuum`)
  is available but opt-in, and adaptive lambda selection does not yet
  support it (falls back to the validated `penalty_order=3`/fixed-`lam`
  default with a warning).
- **Double-peaked/multi-component LOSVDs are not specifically handled.**
  The non-parametric spline recovery can represent them (unlike a fixed
  Gauss-Hermite parametrization), but no notebook in this sequence
  specifically validates that case.

## Installation

```
pip install -e ".[dev]"
```

Optional extras: `.[eso]` (ESO archive downloads), `.[nuts]` (the NUTS/
NumPyro reference uncertainty mode), `.[docs]` (building the docs).

## Testing

```
pytest tests/
```

200+ tests covering LOSVD recovery, moment estimation, multi-template
fitting, Laplace and bootstrap uncertainty (single- and multi-template,
serial and parallel), NUTS reference-mode sampling, noise-scale
calibration and AR(1) correlated-noise whitening, wavelength masking,
iterative cleaning (serial and parallel, including the robust-scale
option), continuum co-fitting, spatial binning and template reduction,
regularization-bias calibration, and feature-composition/input-validation
regressions. These are unit/sanity tests (correct shapes, right order of
magnitude, no crashes, and several real before/after demonstrations
against known synthetic ground truth) -- the full, wide-grid statistical
validation (bias and coverage across the entire physically-relevant
sigma/S-N/resolution range) is in `dev_notes/VALIDATION_CAMPAIGN.md` and
`dev_notes/validation_suite/`, referenced throughout the source
docstrings above.

`A_w`/`y_w` (and `A_stack_w`) must not contain NaN or inf -- mask bad
pixels out of the arrays before calling any fit function; a clear
`ValueError` is raised otherwise rather than a fit silently returning a
NaN-filled LOSVD.

A handful of tests are marked `@pytest.mark.slow` (genuinely
multi-minute -- e.g. a `lam_wing`-active fit refit inside every bootstrap
replicate) and are skipped by default, both via `-m "not slow"` (this
project's own CI path) and by a bundled pytest-astropy-ecosystem plugin
that requires `--run-slow` to execute any `slow`-marked test at all, even
under a plain `pytest tests/` with no `-m` filter. Pass `--run-slow` to
actually run them, e.g. before a release: `pytest tests/ --run-slow`.

## Quick start

One fitting engine, for any number of templates (a single template is
just `template_fluxes = [template_flux]`):

```python
from kinextract import FitConfig, build_template_stack
from kinextract.multi_template import fit_losvd

cfg = FitConfig()          # see FitConfig.describe() for every tunable parameter
v_grid = cfg.build_velocity_grid()
A_stack = build_template_stack(wave_obs, template_wave, template_fluxes, v_grid)  # (n_templates, n_pix, n_bins)
A_stack_w, y_w = A_stack / err[None, :, None], flux / err   # noise-whiten
result = fit_losvd(A_stack_w, y_w, cfg)
recovered_losvd = result["b"]
mixture_fractions = result["w"] / result["w"].sum()  # meaningful for n_templates > 1
```

Run `FitConfig.describe()` to see every tunable parameter grouped by
subsystem, or `FitConfig.describe("lam")` to filter by a substring. See
`examples/notebooks/` for complete, runnable, annotated examples
including uncertainty quantification and diagnostic plots.
