"""Gauss-Hermite characterization of the non-parametric LOSVD.
This module provides tools to *summarize* a non-parametric line-of-sight
velocity distribution (LOSVD) — a histogram ``b`` over ``nl`` velocity bins
on a grid ``xl`` (km/s), recovered by the main spectral-fitting pipeline —
with a compact parametric Gauss-Hermite (GH) model. Following the van der
Marel & Franx (1993) convention,
LOSVD(v) = amp / (sqrt(2*pi) * sigma) * exp(-w^2/2)
* (1 + h3*H3(w) + h4*H4(w) + ...), w = (v - V) / sigma,
this GH fit is *not* the primary fit performed against the galaxy spectrum;
it is a post-hoc least-squares fit of the parametric model to the already-
recovered non-parametric LOSVD, used to compress it into physically
interpretable numbers (V, sigma, h3, h4, ...) for reporting and comparison
across spectra. Two numerically distinct Hermite-polynomial normalizations
are implemented and clearly distinguished: a legacy "fortran-like" 4-moment
convention (``fh3_fortran_like``, ``fh4_fortran_like``,
``gauss_hermite_losvd_model``, ``fit_losvd_gauss_hermite``) matching the
original Fortran pipeline, and the van der Marel & Franx (1993) polynomial
normalization used by pPXF and by the higher-order fitter
(``gauss_hermite_losvd_model_ho``, ``fit_losvd_gauss_hermite_higher``). The
two conventions give numerically different h3/h4 values for the same LOSVD
shape, so they must not be mixed.
"""
from __future__ import annotations
import math
import numpy as np
from scipy.optimize import least_squares
# =============================================================================
# Section 15 - LOSVD analysis
# =============================================================================
[docs]
def fh3_fortran_like(w: np.ndarray) -> np.ndarray:
"""Third-order Gauss-Hermite basis polynomial, legacy Fortran normalization.
Evaluates the h3 basis function used by the original Fortran pipeline
(pallmc.f-style convention), which differs by a numerical normalization
factor from the van der Marel & Franx (1993) polynomial used elsewhere
in this package (see :func:`gauss_hermite_losvd_model_ho`). Used inside
:func:`gauss_hermite_losvd_model` to build the GH LOSVD model and by
:func:`fit_losvd_gauss_hermite` to recover h3 in the fortran convention.
Parameters
----------
w : ndarray
Standardized velocity, ``w = (v - V) / sigma`` (dimensionless).
Returns
-------
ndarray
Value of the fortran-convention H3 basis polynomial at each ``w``.
"""
return (2 * w ** 3 - 3 * w) / np.sqrt(3.0)
[docs]
def fh4_fortran_like(w: np.ndarray) -> np.ndarray:
"""Fourth-order Gauss-Hermite basis polynomial, legacy Fortran normalization.
Evaluates the h4 basis function used by the original Fortran pipeline
(pallmc.f-style convention). See :func:`fh3_fortran_like` for the
companion h3 term and a note on the normalization convention; this
fortran convention differs numerically from the van der Marel & Franx
(1993) polynomials used by :func:`gauss_hermite_losvd_model_ho`.
Parameters
----------
w : ndarray
Standardized velocity, ``w = (v - V) / sigma`` (dimensionless).
Returns
-------
ndarray
Value of the fortran-convention H4 basis polynomial at each ``w``.
"""
return (4 * w ** 4 - 12 * w ** 2 + 3) / np.sqrt(24.0)
def _gh_poly_vdm(n: int, w: np.ndarray) -> np.ndarray:
"""Normalized Gauss-Hermite polynomial of order n (van der Marel & Franx 1993).
H_n(w) = He_n(w) / sqrt(n!) where He_n are probabilist Hermite polynomials
generated by He_{k+1} = w*He_k - k*He_{k-1}. Standard pPXF/GH literature convention.
"""
w = np.asarray(w, float)
if n == 0:
return np.ones_like(w)
if n == 1:
return w.copy()
he_prev = np.ones_like(w)
he_curr = w.copy()
for k in range(1, n):
he_next = w * he_curr - k * he_prev
he_prev = he_curr
he_curr = he_next
return he_curr / np.sqrt(float(math.factorial(n)))
[docs]
def gauss_hermite_losvd_model_ho(
v: np.ndarray, amp: float, vel: float, sig: float, *hn: float,
) -> np.ndarray:
"""Higher-order Gauss-Hermite LOSVD model, van der Marel & Franx (1993) convention.
Evaluates ``amp * exp(-w^2/2) / (sqrt(2*pi)*sig) * (1 + sum_i hn[i] * H_{i+3}(w))``
where ``w = (v - vel) / sig`` and ``H_n`` are the vdM1993-normalized
Gauss-Hermite polynomials (see :func:`_gh_poly_vdm`). Unlike
:func:`gauss_hermite_losvd_model`, this supports an arbitrary number of
higher-order moments (h3, h4, h5, h6, ...) passed positionally via
``*hn``, enabling a more flexible characterization of non-Gaussian LOSVD
wings. Uses the same polynomial normalization as pPXF, which is
numerically distinct from the legacy fortran convention used by
:func:`gauss_hermite_losvd_model` — do not compare h3/h4 values between
the two model functions directly.
Parameters
----------
v : ndarray
Velocities (km/s) at which to evaluate the model.
amp : float
Overall normalization (total flux/area under the LOSVD).
vel : float
Mean velocity V (km/s).
sig : float
Velocity dispersion sigma (km/s). Must be positive and finite; if
not, the model returns an array of NaN with the same shape as ``v``.
*hn : float
Higher-order Gauss-Hermite moments h3, h4, h5, ... (dimensionless),
in ascending order starting at h3.
Returns
-------
ndarray
Model LOSVD amplitude evaluated at each ``v``, same shape as ``v``.
All-NaN if ``sig`` is non-positive or non-finite.
"""
sig = float(sig)
if sig <= 0 or not np.isfinite(sig):
return np.full_like(v, np.nan, float)
w = (np.asarray(v, float) - vel) / sig
poly_sum = np.ones_like(w, float)
for i, h in enumerate(hn):
poly_sum = poly_sum + float(h) * _gh_poly_vdm(i + 3, w)
return amp * np.exp(-0.5 * w ** 2) / (np.sqrt(2 * np.pi) * sig) * poly_sum
[docs]
def fit_losvd_gauss_hermite_higher(
v: np.ndarray, y: np.ndarray, max_order: int = 6,
) -> dict:
"""Fit a higher-order Gauss-Hermite model (h3..h_max_order) to a LOSVD.
Least-squares fits :func:`gauss_hermite_losvd_model_ho` to the sampled,
non-parametric LOSVD ``y(v)``, using van der Marel & Franx (1993)
normalized Hermite polynomials (vdM1993) — the same convention as pPXF.
Including higher-order terms beyond h3/h4 lets the fit absorb
non-Gaussian wings that would otherwise bias a truncated 4-moment fit
(:func:`fit_losvd_gauss_hermite`), typically giving more accurate
recovered V and sigma at the cost of extra nuisance parameters. Initial
guesses for amplitude, center, and width are derived from
:func:`getfwhm_fortran_like`; all higher-order ``h`` terms start at zero.
Bounded least-squares (`scipy.optimize.least_squares`) is used with log
parametrization of sigma to keep it positive.
Parameters
----------
v : ndarray
Velocity grid (km/s) on which the LOSVD ``y`` is sampled.
y : ndarray
Non-parametric LOSVD amplitude at each ``v`` (same shape as ``v``).
max_order : int, optional
Highest Gauss-Hermite moment order to fit (minimum enforced value
is 4, i.e. at least h3 and h4 are always fit). Default is 6, fitting
h3, h4, h5, h6.
Returns
-------
dict
Dictionary with keys:
- ``vherm``, ``sherm`` : fitted mean velocity and dispersion (km/s).
- ``amp`` : fitted overall amplitude/normalization.
- ``h3``, ``h4``, ... up to ``h{max_order}`` : fitted higher-order
Gauss-Hermite moments (dimensionless, vdM1993 convention — these
differ numerically from the fortran-convention h3/h4 returned by
:func:`fit_losvd_gauss_hermite`).
- ``model`` : the best-fit model evaluated on the full input ``v``
grid (NaN where the fit could not be performed).
- ``fit_success`` : bool, whether the least-squares solver converged.
- ``fit_message`` : str, solver status message (or failure reason).
- ``max_order`` : int, the (possibly clamped) max order actually used.
- ``convention`` : str, always ``"vdM1993"`` for this function.
If fewer than 5 finite ``(v, y)`` points are available, all numeric
keys are NaN and ``fit_success`` is False.
"""
v, y = np.asarray(v, float), np.asarray(y, float)
ok = np.isfinite(v) & np.isfinite(y)
vfit, yfit = v[ok], y[ok]
max_order = max(max_order, 4)
n_hm = max_order - 2 # h3..h_{max_order}
ho_keys = [f"h{n}" for n in range(3, max_order + 1)]
nan_result: dict = {k: np.nan for k in ["vherm", "sherm", "amp"] + ho_keys}
nan_result.update({"model": np.full_like(v, np.nan), "fit_success": False,
"fit_message": "Too few points", "max_order": max_order,
"convention": "vdM1993"})
if vfit.size < 5:
return nan_result
f = getfwhm_fortran_like(vfit, yfit)
fwhm_sigma = (
f["fwhm"] / 2.35 if np.isfinite(f["fwhm"]) and f["fwhm"] > 0
else max(float(np.std(vfit)), 1.0)
)
amp0 = f["ymax"] * np.sqrt(2 * np.pi) * fwhm_sigma
if not np.isfinite(amp0) or amp0 <= 0:
amp0 = max(float(np.sum(np.maximum(yfit, 0))), 1e-6)
vel0 = f["vmax"] if np.isfinite(f["vmax"]) else float(vfit[np.argmax(yfit)])
sig0 = max(fwhm_sigma, 1.0)
vspan = max(float(np.ptp(vfit)), float(np.std(vfit)), 1.0)
sigy = 0.01
p0 = [amp0, vel0, np.log(sig0)] + [0.0] * n_hm
lo = [0.0, vfit.min() - 2 * vspan, np.log(1e-6)] + [-0.5] * n_hm
hi = [np.inf, vfit.max() + 2 * vspan, np.log(1e6)] + [0.5] * n_hm
def _resid(p: np.ndarray) -> np.ndarray:
return (
gauss_hermite_losvd_model_ho(vfit, p[0], p[1], np.exp(p[2]), *p[3:]) - yfit
) / sigy
try:
res = least_squares(_resid, p0, bounds=(lo, hi), max_nfev=10000)
except Exception as exc:
nan_result["fit_message"] = str(exc)
return nan_result
amp_fit = float(res.x[0])
vel_fit = float(res.x[1])
sig_fit = float(np.exp(res.x[2]))
hn_fit = res.x[3:]
result: dict = {
"vherm": vel_fit, "sherm": sig_fit, "amp": amp_fit,
"convention": "vdM1993",
"max_order": max_order,
"model": gauss_hermite_losvd_model_ho(v, amp_fit, vel_fit, sig_fit, *hn_fit),
"fit_success": bool(res.success), "fit_message": str(res.message),
}
for i, h in enumerate(hn_fit, start=3):
result[f"h{i}"] = float(h)
return result
[docs]
def gauss_hermite_losvd_model(
v: np.ndarray, amp: float, vel: float, sig: float, h3: float, h4: float,
) -> np.ndarray:
"""Four-moment Gauss-Hermite LOSVD model, legacy Fortran normalization.
Evaluates the classic van der Marel & Franx (1993) LOSVD parametrization
truncated at h4:
``amp/sqrt(2*pi)/sig * exp(-w^2/2) * (1 + h3*H3(w) + h4*H4(w))`` with
``w = (v - vel) / sig``, using the legacy "fortran-like" H3/H4 basis
polynomials (:func:`fh3_fortran_like`, :func:`fh4_fortran_like`) rather
than the vdM1993 polynomial normalization used by
:func:`gauss_hermite_losvd_model_ho`. This is the model fit by
:func:`fit_losvd_gauss_hermite` and is kept numerically consistent with
the original Fortran pipeline (pallmc.f) for backward compatibility;
h3/h4 values from this function are **not** directly comparable to
h3/h4 from the ``_ho``/vdM1993 functions.
Parameters
----------
v : ndarray
Velocities (km/s) at which to evaluate the model.
amp : float
Overall normalization (total flux/area under the LOSVD).
vel : float
Mean velocity V (km/s).
sig : float
Velocity dispersion sigma (km/s). Must be positive and finite; if
not, the model returns an array of NaN with the same shape as ``v``.
h3 : float
Third Gauss-Hermite moment (dimensionless, fortran convention) —
skewness-like asymmetry of the LOSVD.
h4 : float
Fourth Gauss-Hermite moment (dimensionless, fortran convention) —
kurtosis-like symmetric deviation from a Gaussian.
Returns
-------
ndarray
Model LOSVD amplitude evaluated at each ``v``, same shape as ``v``.
All-NaN if ``sig`` is non-positive or non-finite.
"""
sig = float(sig)
if sig <= 0 or not np.isfinite(sig):
return np.full_like(v, np.nan, float)
w = (np.asarray(v, float) - vel) / sig
return (
amp * np.exp(-0.5 * w ** 2) / np.sqrt(2 * np.pi) / sig
* (1 + h3 * fh3_fortran_like(w) + h4 * fh4_fortran_like(w))
)
[docs]
def getfwhm_fortran_like(x: np.ndarray, y: np.ndarray, frac: float = 0.5) -> dict:
"""Estimate the full-width-at-fraction-of-max and flux moments of a curve.
Locates the peak of ``y(x)`` and linearly interpolates the crossing
points where ``y`` falls to ``frac`` (default: half) of its peak value
on either side, giving a robust width estimate that does not require a
parametric fit. Also returns the simple flux-weighted first and second
moments (v1, v2) of ``y`` as a model-free estimate of the centroid and
width. Used to construct good initial guesses (amplitude, center,
sigma) for the Gauss-Hermite least-squares fits
(:func:`fit_losvd_gauss_hermite`, :func:`fit_losvd_gauss_hermite_higher`),
matching the legacy Fortran pipeline's FWHM-based initialization scheme.
Parameters
----------
x : ndarray
Independent variable, typically the LOSVD velocity grid (km/s).
y : ndarray
Dependent variable (e.g. LOSVD amplitude) sampled at each ``x``.
Assumed single-peaked; only the interval on either side of the
global maximum is searched for the half-max crossings.
frac : float, optional
Fraction of the peak value at which to measure the width (0.5 for
standard FWHM). Default is 0.5.
Returns
-------
dict
Dictionary with keys:
- ``fwhm`` : full width at ``frac`` of maximum (km/s if ``x`` is
velocity).
- ``vmax`` : location of the peak (same units as ``x``).
- ``ymax`` : peak value of ``y``.
- ``x_left``, ``x_right`` : interpolated crossing points of the
``frac``-of-max level on the left/right of the peak.
- ``v1`` : flux-weighted first moment (centroid) of ``y`` over ``x``.
- ``v2`` : flux-weighted second moment (RMS width) of ``y`` about
``v1``.
All values are NaN if fewer than 3 points are supplied.
"""
x, y = np.asarray(x, float), np.asarray(y, float)
if x.size < 3:
return {k: np.nan for k in ["fwhm", "vmax", "ymax", "x_left", "x_right", "v1", "v2"]}
imax = int(np.argmax(y))
ymax = float(y[imax]); vmax = float(x[imax]); yhalf = ymax * frac
x1 = float(x[0])
for i in range(imax):
if yhalf >= y[i] and yhalf < y[i + 1]:
x1 = x[i] + (yhalf - y[i]) / (y[i + 1] - y[i]) * (x[i + 1] - x[i])
x2 = float(x[-1])
for i in range(imax, len(x) - 1):
if yhalf >= y[i + 1] and yhalf < y[i]:
x2 = x[i + 1] + (yhalf - y[i + 1]) / (y[i] - y[i + 1]) * (x[i] - x[i + 1])
sumy = float(np.sum(y))
v1 = float(np.sum(x * y) / sumy) if sumy else np.nan
v2 = float(np.sqrt(np.sum(y * (x - v1) ** 2) / sumy)) if sumy else np.nan
return {
"fwhm": float(x2 - x1), "vmax": vmax, "ymax": ymax,
"x_left": x1, "x_right": x2, "v1": v1, "v2": v2,
}
[docs]
def fit_losvd_gauss_hermite(v: np.ndarray, y: np.ndarray, fit_h3h4: bool = True) -> dict:
"""Fit the classic 4-moment Gauss-Hermite model to a non-parametric LOSVD.
This is the standard, publication-facing GH characterization of the
recovered non-parametric LOSVD: a bounded least-squares fit of
:func:`gauss_hermite_losvd_model` (legacy fortran H3/H4 normalization,
matching the original Fortran pipeline pallmc.f) to ``y(v)``, giving the
mean velocity V, dispersion sigma, and the skewness/kurtosis-like
moments h3 and h4. Initial guesses are derived from
:func:`getfwhm_fortran_like`, and sigma is fit in log-space to enforce
positivity. This is the fitter used by :func:`kinextract.plotting.plot_losvd`
and by :class:`kinextract.errors.LOSVDErrorEstimator` to characterize both
the MAP LOSVD and each bootstrap replicate. For higher-order (h5, h6, ...)
characterization using the vdM1993 convention instead, see
:func:`fit_losvd_gauss_hermite_higher`.
Parameters
----------
v : ndarray
Velocity grid (km/s) on which the LOSVD ``y`` is sampled.
y : ndarray
Non-parametric LOSVD amplitude at each ``v`` (same shape as ``v``).
fit_h3h4 : bool, optional
If True (default), fit h3 and h4 along with amplitude, V, and sigma.
If False, h3 and h4 are fixed to zero and only a pure Gaussian
(amp, V, sigma) is fit — useful for comparing against a
strictly-Gaussian LOSVD assumption.
Returns
-------
dict
Dictionary with keys:
- ``vherm``, ``sherm`` : fitted mean velocity and dispersion (km/s).
- ``h3``, ``h4`` : fitted Gauss-Hermite moments (dimensionless,
fortran convention — see :func:`fh3_fortran_like`); zero if
``fit_h3h4=False``.
- ``amp`` : fitted overall amplitude/normalization.
- ``fwhm_sigma`` : sigma estimate derived from the FWHM initial
guess (km/s), useful as a sanity check against ``sherm``.
- ``v1``, ``v2`` : flux-weighted centroid and RMS width of the raw
LOSVD (km/s), model-free moments from :func:`getfwhm_fortran_like`.
- ``fwhm``, ``vmax``, ``ymax``, ``x_left``, ``x_right`` : raw
half-max diagnostics from :func:`getfwhm_fortran_like`.
- ``model`` : the best-fit GH model evaluated on the full input
``v`` grid (NaN where the fit could not be performed).
- ``fit_success`` : bool, whether the least-squares solver converged.
- ``fit_message`` : str, solver status message (or failure reason).
If fewer than 5 finite ``(v, y)`` points are available, all numeric
keys are NaN and ``fit_success`` is False.
"""
v, y = np.asarray(v, float), np.asarray(y, float)
ok = np.isfinite(v) & np.isfinite(y)
vfit, yfit = v[ok], y[ok]
nan_result = {k: np.nan for k in [
"vherm", "sherm", "fwhm_sigma", "h3", "h4",
"v1", "v2", "fwhm", "vmax", "ymax", "x_left", "x_right", "amp",
]}
if vfit.size < 5:
return {**nan_result, "model": np.full_like(v, np.nan), "fit_success": False,
"fit_message": "Too few points"}
f = getfwhm_fortran_like(vfit, yfit)
fwhm_sigma = (
f["fwhm"] / 2.35 if np.isfinite(f["fwhm"]) and f["fwhm"] > 0
else max(float(np.std(vfit)), 1.0)
)
amp0 = f["ymax"] * np.sqrt(2 * np.pi) * fwhm_sigma
if not np.isfinite(amp0) or amp0 <= 0:
try:
amp0 = np.trapezoid(np.maximum(yfit, 0), vfit)
except AttributeError:
amp0 = np.trapz(np.maximum(yfit, 0), vfit)
if not np.isfinite(amp0) or amp0 <= 0:
amp0 = max(float(np.sum(yfit)), 1e-6)
vel0 = f["vmax"] if np.isfinite(f["vmax"]) else float(vfit[np.argmax(yfit)])
sig0 = max(
fwhm_sigma if np.isfinite(fwhm_sigma) and fwhm_sigma > 0 else float(np.std(vfit)),
1.0,
)
vspan = max(float(np.ptp(vfit)), float(np.std(vfit)), 1.0)
sigy = 0.01
if fit_h3h4:
p0 = [amp0, vel0, np.log(sig0), 0.0, 0.0]
lo = [0.0, vfit.min() - 2 * vspan, np.log(1e-6), -1.0, -1.0]
hi = [np.inf, vfit.max() + 2 * vspan, np.log(1e6), 1.0, 1.0]
def resid(p):
return (
gauss_hermite_losvd_model(vfit, p[0], p[1], np.exp(p[2]), p[3], p[4]) - yfit
) / sigy
res = least_squares(resid, p0, bounds=(lo, hi), max_nfev=10000)
amp, vel, log_sig, h3, h4 = res.x
sig = np.exp(log_sig)
else:
p0 = [amp0, vel0, np.log(sig0)]
lo = [0.0, vfit.min() - 2 * vspan, np.log(1e-6)]
hi = [np.inf, vfit.max() + 2 * vspan, np.log(1e6)]
def resid(p):
return (
gauss_hermite_losvd_model(vfit, p[0], p[1], np.exp(p[2]), 0.0, 0.0) - yfit
) / sigy
res = least_squares(resid, p0, bounds=(lo, hi), max_nfev=10000)
amp, vel, log_sig = res.x
sig = np.exp(log_sig)
h3 = h4 = 0.0
return {
"vherm": float(vel), "sherm": float(sig), "fwhm_sigma": float(fwhm_sigma),
"h3": float(h3), "h4": float(h4),
"v1": float(f["v1"]), "v2": float(f["v2"]),
"fwhm": float(f["fwhm"]), "vmax": float(f["vmax"]), "ymax": float(f["ymax"]),
"x_left": float(f["x_left"]), "x_right": float(f["x_right"]), "amp": float(amp),
"model": gauss_hermite_losvd_model(v, amp, vel, sig, h3, h4),
"fit_success": bool(res.success), "fit_message": str(res.message),
}