"""Peak fitting of 1D spectroscopic data on top of silx's fitting engine.
This module is a thin headless wrapper around :mod:`silx.math.fit` -- the engine
behind :class:`silx.gui.fit.FitWidget` -- fitting one *or several* peaks of a
single shape (Gaussian, Lorentzian or pseudo-Voigt, symmetric or
split/asymmetric) together with a silx background (constant, linear, or the
Strip/Snip peak-stripping filters), and reporting the physically meaningful peak
descriptors (center, height, area, FWHM and the pseudo-Voigt mixing) for each
fitted peak.
The fit is exactly silx's native flow (``estimate`` -- which auto-detects the
peak count via silx's peak search -- followed by ``runfit``), so a headless
batch fit reproduces what the interactive
:class:`~silx.gui.fit.FitWidget` computes for the same configuration. Surplus
peaks the fit drives to zero are pruned. The
:class:`~silx.math.fit.fitmanager.FitManager` and the curated theory/background
set come from :mod:`ewoksxas.fit.silx_engine`, shared with the Orange widget.
"""
from __future__ import annotations
import copy
import logging
from enum import Enum
from functools import cached_property
from typing import TYPE_CHECKING
import numpy as np
from numpy.linalg import LinAlgError
from ewoksxas.fit.data import XYData, r2, rmse
from ewoksxas.fit.silx_engine import (
background_curve,
descriptors_from_fit_results,
make_fit_manager,
)
if TYPE_CHECKING:
from silx.math.fit.fitmanager import FitManager
from ewoksxas.fit.data import ArrayF64, DataDict
__all__ = [
"Baseline",
"MultiPeakFitResult",
"PeakModel",
"PeakShape",
"peak_descriptors",
]
logger = logging.getLogger(__name__)
[docs]
class PeakShape(Enum):
"""A peak shape and the silx fit theory that realises it.
Symmetric and asymmetric variants are separate members. Each member carries,
in order, the silx theory name, the profile *family* (used to derive the area
from the height and to decide whether an ``eta`` mixing fraction is reported)
and whether the shape is asymmetric (reports per-side FWHMs). The lowercase
``label`` (used in config and the GUI) is derived from the member name.
silx provides pseudo-Voigt (not a true Voigt) and models asymmetry with two
half-widths ("Split ..."), so ``VOIGT`` maps to a pseudo-Voigt and the skewed
shapes map to split profiles.
"""
# (silx theory name, family, asymmetric)
GAUSSIAN = ("Gaussians", "gaussian", False)
SKEWED_GAUSSIAN = ("Split Gaussian", "gaussian", True)
LORENTZIAN = ("Lorentz", "lorentzian", False)
SPLIT_LORENTZIAN = ("Split Lorentz", "lorentzian", True)
VOIGT = ("Pseudo-Voigt Line", "voigt", False)
SKEWED_VOIGT = ("Split Pseudo-Voigt", "voigt", True)
def __init__(self, theory_name: str, family: str, asymmetric: bool) -> None:
self.label = self.name.lower().replace("_", " ")
self.theory_name = theory_name
self.family = family
self.asymmetric = asymmetric
[docs]
@classmethod
def from_label(cls, value: str | PeakShape) -> PeakShape:
"""Resolve a label string (or pass a member through) to a PeakShape.
Raises:
ValueError: If the label does not match any shape.
"""
if isinstance(value, cls):
return value
for shape in cls:
if shape.label == value:
return shape
raise ValueError(f"Unknown shape: {value!r}.")
[docs]
@classmethod
def from_theory_name(cls, theory_name: str) -> PeakShape:
"""Resolve a silx theory name (e.g. ``"Gaussians"``) to a PeakShape.
Raises:
ValueError: If no shape maps to the theory name.
"""
for shape in cls:
if shape.theory_name == theory_name:
return shape
raise ValueError(f"Unknown silx theory: {theory_name!r}.")
@classmethod
def _missing_(cls, value: object) -> PeakShape | None:
# Backward-compat: lmfit-era pickles store (model_class, family, skew_param)
if isinstance(value, tuple) and len(value) == 3 and isinstance(value[0], type):
_, family, skew_param = value
asymmetric = skew_param is not None
for shape in cls:
if shape.family == family and shape.asymmetric == asymmetric:
return shape
return None
[docs]
class Baseline(Enum):
"""A baseline option and the silx background theory it maps to.
The background is fitted by silx together with the peaks. ``STRIP`` and
``SNIP`` are silx's iterative peak-stripping backgrounds, suited to
multi-peak spectra; their filter settings come from the fit configuration
(see :class:`PeakModel`).
"""
# (label, silx background theory name)
NONE = ("none", "No Background")
OFFSET = ("offset", "Constant")
LINEAR = ("linear", "Linear")
STRIP = ("strip", "Strip")
SNIP = ("snip", "Snip")
def __init__(self, label: str, theory_name: str) -> None:
self.label = label
self.theory_name = theory_name
[docs]
@classmethod
def from_label(cls, value: str | Baseline) -> Baseline:
"""Resolve a label string (or pass a member through) to a Baseline.
Raises:
ValueError: If the label does not match any baseline.
"""
if isinstance(value, cls):
return value
for baseline in cls:
if baseline.label == value:
return baseline
raise ValueError(f"Unknown baseline: {value!r}.")
[docs]
@classmethod
def from_theory_name(cls, theory_name: str) -> Baseline:
"""Resolve a silx background name (e.g. ``"Constant"``) to a Baseline.
Raises:
ValueError: If no baseline maps to the background name.
"""
for baseline in cls:
if baseline.theory_name == theory_name:
return baseline
raise ValueError(f"Unknown silx background: {theory_name!r}.")
@classmethod
def _missing_(cls, value: object) -> Baseline | None:
# Backward-compat: lmfit-era pickles store (label, model_or_None).
# Shirley was removed in the silx migration -> degrade to no background.
if isinstance(value, tuple) and len(value) == 2 and isinstance(value[0], str):
label = value[0]
for baseline in cls:
if baseline.label == label:
return baseline
if label == "shirley":
return cls.NONE
return None
# Per-peak descriptor keys common to every shape, in GUI fit-table order.
_COMMON_DESCRIPTORS = (
"center",
"center_stderr",
"height",
"area",
"area_stderr",
)
[docs]
def peak_descriptors(shape: PeakShape) -> tuple[str, ...]:
"""Return the per-peak descriptor keys emitted for a given peak shape.
Symmetric shapes report a single ``fwhm``; asymmetric (skewed/split) shapes
report the per-side ``fwhm_low``/``fwhm_high`` instead. Voigt-family shapes
additionally report the pseudo-Voigt mixing fraction ``eta``. This is the
single place that defines the per-shape column set, shared with the task.
"""
widths = ("fwhm_low", "fwhm_high") if shape.asymmetric else ("fwhm",)
mixing = ("eta",) if shape.family == "voigt" else ()
return _COMMON_DESCRIPTORS + widths + mixing
# A peak whose height falls below this fraction of the tallest peak is a
# vanishing (degenerate) artefact of over-fitting and disqualifies its fit.
_DEGENERATE_HEIGHT_FRAC = 1e-3
[docs]
class MultiPeakFitResult:
"""The result of an N-peak fit, exposing per-peak physical descriptors.
The bare scalar properties (``center``, ``area``, ...) report the first
(leftmost) peak so that single-peak callers keep working unchanged; the full
per-peak descriptors are available through :attr:`peaks`.
"""
def __init__(self, manager: FitManager, shape: PeakShape, roi: XYData) -> None:
"""Wrap a silx fit and expose the derived peak descriptors.
Args:
manager: The silx ``FitManager`` after ``estimate``/``runfit`` on the
region-of-interest data (background and peaks fitted together).
shape: The fitted peak shape.
roi: The fitted region-of-interest data.
"""
self.shape: PeakShape = shape
self._manager = manager
model_y = np.asarray(manager.gendata(), dtype=np.float64)
baseline_y = background_curve(manager)
self.model = XYData(x=roi.x, y=model_y)
self.baseline = XYData(x=roi.x, y=baseline_y)
self.peak = XYData(x=roi.x, y=model_y - baseline_y)
self._observed = roi.y
# Number of varied parameters, for the degrees of freedom.
self._n_varied = sum(
1
for entry in manager.fit_results
if entry["code"] not in ("FIXED", "IGNORE")
)
# Order peaks left to right so peak 0 is the leftmost (stable indexing).
peaks = descriptors_from_fit_results(
manager.fit_results, shape.family, shape.asymmetric
)
peaks.sort(key=lambda descriptors: descriptors["center"])
self._peaks: list[dict[str, float]] = peaks
@property
def peaks(self) -> list[dict[str, float]]:
"""Per-peak descriptor dictionaries, ordered by ascending center."""
return self._peaks
@property
def n_peaks(self) -> int:
"""Number of fitted peaks."""
return len(self._peaks)
@property
def center(self) -> float:
"""Center x0 of the first (leftmost) peak."""
return self._peaks[0]["center"]
@property
def area(self) -> float:
"""Integrated area of the first peak."""
return self._peaks[0]["area"]
@property
def height(self) -> float:
"""Maximum height of the first peak above the baseline."""
return self._peaks[0]["height"]
@property
def fwhm(self) -> float:
"""Total full width at half maximum of the first (symmetric) peak."""
return self._peaks[0].get("fwhm", float("nan"))
@property
def fwhm_low(self) -> float:
"""Low-side full width at half maximum of the first (asymmetric) peak."""
return self._peaks[0].get("fwhm_low", float("nan"))
@property
def fwhm_high(self) -> float:
"""High-side full width at half maximum of the first (asymmetric) peak."""
return self._peaks[0].get("fwhm_high", float("nan"))
@property
def eta(self) -> float:
"""Pseudo-Voigt mixing fraction of the first peak (Voigt shapes only)."""
return self._peaks[0].get("eta", float("nan"))
[docs]
def without_negligible_peaks(self) -> MultiPeakFitResult:
"""Return a result with vanishing peaks dropped (for automatic counting).
silx may over-detect peaks; a surplus peak's height and area are driven
towards zero by the fit. Those are removed so the output reports only the
peaks the data supports. A zeroed peak's removal barely perturbs the
model curve, so no re-fit is needed; at least the tallest peak is kept,
and ``self`` is returned when nothing is pruned.
"""
if self.n_peaks <= 1:
return self
heights = [peak["height"] for peak in self._peaks]
areas = [peak["area"] for peak in self._peaks]
tallest = max(heights)
largest_area = max(areas)
survivors = [
peak
for peak, height, area in zip(self._peaks, heights, areas, strict=True)
if height >= _DEGENERATE_HEIGHT_FRAC * tallest
and area >= _DEGENERATE_HEIGHT_FRAC * largest_area
]
if len(survivors) == self.n_peaks:
return self
if not survivors:
survivors = [self._peaks[int(np.argmax(heights))]]
pruned = copy.copy(self)
pruned._peaks = survivors
return pruned
@cached_property
def rmse(self) -> float:
"""Root-mean-square error of the total model against the data."""
return rmse(self._observed, self.model.y)
@cached_property
def r2(self) -> float:
"""Coefficient of determination of the total model."""
return r2(self._observed, self.model.y)
@property
def dof(self) -> int:
"""Degrees of freedom of the fit (data points minus varied params)."""
return int(self._observed.size - self._n_varied)
[docs]
def as_dict(self) -> dict[str, float]:
"""Return the first peak's descriptors and the fit quality as a dict."""
return {**self._peaks[0], "r2": self.r2, "rmse": self.rmse}
[docs]
def summary(self) -> str:
"""Return a formatted report of the fit."""
return repr(self)
def __repr__(self) -> str:
ilen = len("FWHM (high):")
first = self._peaks[0]
lines = [
"-" * 72,
f"Peak Fit Results ({self.shape.label}, {self.n_peaks} peak(s))",
"-" * 72,
f"{'Center:':>{ilen}} {first['center']:.6g}",
f"{'Height:':>{ilen}} {first['height']:.6g}",
f"{'Area:':>{ilen}} {first['area']:.6g}",
]
# Only the descriptors emitted for this shape are shown.
width_labels = {
"fwhm": "FWHM:",
"fwhm_low": "FWHM (low):",
"fwhm_high": "FWHM (high):",
"eta": "eta:",
}
for key, label in width_labels.items():
if key in first:
lines.append(f"{label:>{ilen}} {first[key]:.6g}")
lines.append(f"{'R2:':>{ilen}} {self.r2:.6g}")
lines.append(f"{'RMSE:':>{ilen}} {self.rmse:.6g}")
return "\n".join(lines) + "\n"
# Default silx peak-search sensitivity (matches silx's own default).
_DEFAULT_SENSITIVITY = 2.5
# Relative size of the deterministic jitter used to break the optimiser
# degeneracy that perfectly-noiseless data can trigger in silx's least squares.
_JITTER_REL = 1e-7
# Peak-search FWHM (in points) used by the last-resort retry when silx's
# automatic FWHM estimation leads to a singular system.
_FALLBACK_FWHM_POINTS = 20
def _silent(data: dict | None = None) -> None:
"""No-op callback so silx's estimate/runfit error path does not crash.
silx calls the (optional) progress callback unconditionally in its
``LinAlgError`` handler; passing this no-op lets the real error propagate.
"""
[docs]
class PeakModel:
"""Fit one or several peaks of a single shape, with an optional baseline."""
def __init__( # noqa: PLR0913, PLR0917
self,
data: XYData,
shape: str | PeakShape = PeakShape.GAUSSIAN,
baseline: str | Baseline = Baseline.NONE,
roi: tuple[float, float] | None = None,
sensitivity: float = _DEFAULT_SENSITIVITY,
fwhm_points: int | None = None,
auto_fwhm: bool = True,
strip_width: int | None = None,
strip_iterations: int | None = None,
snip_width: int | None = None,
) -> None:
"""Initialise a peak fit for a single x/y observation.
Args:
data: The x and y data to fit.
shape: Peak shape, as a :class:`PeakShape` or its label, e.g.
"gaussian", "skewed gaussian", "lorentzian", "split lorentzian",
"voigt" or "skewed voigt".
baseline: Baseline, as a :class:`Baseline` or its label: "none",
"offset", "linear", "strip" or "snip".
roi: Optional ``(xmin, xmax)`` region of interest; data outside it is
excluded from the fit.
sensitivity: silx peak-search sensitivity (detection threshold,
higher detects fewer peaks).
fwhm_points: silx peak-search FWHM in points (used when
``auto_fwhm`` is False).
auto_fwhm: When True (default) silx estimates the search FWHM from the
data.
strip_width: Strip-background filter width (silx ``StripWidth``).
strip_iterations: Strip-background iterations (``StripIterations``).
snip_width: Snip-background filter width (silx ``SnipWidth``).
Raises:
ValueError: If the shape or baseline is unknown, or the region of
interest contains fewer than three points.
"""
self.shape: PeakShape = PeakShape.from_label(shape)
self.baseline: Baseline = Baseline.from_label(baseline)
self.sensitivity: float = float(sensitivity)
self.fwhm_points: int | None = None if fwhm_points is None else int(fwhm_points)
self.auto_fwhm: bool = bool(auto_fwhm)
self.strip_width = None if strip_width is None else int(strip_width)
self.strip_iterations = (
None if strip_iterations is None else int(strip_iterations)
)
self.snip_width = None if snip_width is None else int(snip_width)
xmin, xmax = roi if roi is not None else (-np.inf, np.inf)
self.data: XYData = data
self.xmin: float = max(float(data.x.min()), xmin)
self.xmax: float = min(float(data.x.max()), xmax)
self.roi: XYData = self._set_roi()
if self.roi.x.size < 3:
raise ValueError(f"ROI has too few points to fit: {self.roi.x.size}.")
def _set_roi(self) -> XYData:
"""Return the data restricted to the region of interest."""
if self.xmin == self.data.x.min() and self.xmax == self.data.x.max():
return self.data
mask = (self.data.x >= self.xmin) & (self.data.x <= self.xmax)
return XYData(x=self.data.x[mask], y=self.data.y[mask])
def _make_manager(self) -> FitManager:
"""Build a silx FitManager configured for this shape and baseline."""
manager = make_fit_manager()
manager.settheory(self.shape.theory_name)
manager.setbackground(self.baseline.theory_name)
config: dict[str, object] = {
"Sensitivity": self.sensitivity,
"AutoFwhm": self.auto_fwhm,
"ForcePeakPresence": True,
"PositiveHeightAreaFlag": True,
"PositiveFwhmFlag": True,
}
if self.fwhm_points is not None:
config["FwhmPoints"] = self.fwhm_points
if self.strip_width is not None:
config["StripWidth"] = self.strip_width
if self.strip_iterations is not None:
config["StripIterations"] = self.strip_iterations
if self.snip_width is not None:
config["SnipWidth"] = self.snip_width
manager.configure(**config)
return manager
def _fit(self, y: ArrayF64, **extra_config: object) -> FitManager:
"""Run silx's native estimate + runfit on ``y`` over the ROI grid."""
manager = self._make_manager()
if extra_config:
manager.configure(**extra_config)
manager.setdata(x=self.roi.x, y=np.ascontiguousarray(y, dtype=np.float64))
manager.estimate(callback=_silent)
manager.runfit(callback=_silent)
return manager
def _robust_fit(self) -> FitManager:
"""Fit with a retry ladder for silx's singular-system corner cases.
Perfectly-noiseless or very smooth data can make silx's estimation or
least squares singular. Retries escalate: a tiny deterministic jitter
(breaks the degeneracy without affecting the recovered parameters),
then a fixed peak-search FWHM instead of silx's automatic estimate.
"""
y = self.roi.y
rng = np.random.default_rng(0)
scale = _JITTER_REL * max(1.0, float(np.max(np.abs(y))))
jittered = y + rng.standard_normal(y.size) * scale
fixed_fwhm = {"AutoFwhm": False, "FwhmPoints": _FALLBACK_FWHM_POINTS}
attempts = ((y, {}), (jittered, {}), (y, fixed_fwhm), (jittered, fixed_fwhm))
last_error: LinAlgError | None = None
for values, extra_config in attempts:
try:
return self._fit(values, **extra_config)
except LinAlgError as error: # noqa: PERF203
last_error = error
if last_error is not None:
raise last_error
raise AssertionError() from last_error
[docs]
def fit(self) -> MultiPeakFitResult:
"""Fit the peaks and background to the region of interest.
Runs silx's native flow -- ``estimate`` (auto peak count) then
``runfit`` with the selected background theory -- and prunes peaks the
fit drove to zero.
Returns:
The fit wrapped in a :class:`MultiPeakFitResult`.
"""
manager = self._robust_fit()
result = MultiPeakFitResult(manager=manager, shape=self.shape, roi=self.roi)
return result.without_negligible_peaks()
[docs]
@classmethod
def from_dict(cls, data: DataDict, parameters: dict | None = None) -> PeakModel:
"""Construct a PeakModel from plain dictionaries.
Args:
data: Mapping with ``"x"`` and ``"y"`` arrays.
parameters: Optional settings: ``shape``, ``baseline``, ``xmin``,
``xmax``, ``sensitivity``, ``fwhm_points``, ``auto_fwhm``,
``strip_width``, ``strip_iterations`` and ``snip_width``.
"""
parameters = parameters or {}
xmin = parameters.get("xmin")
xmax = parameters.get("xmax")
roi = None
if xmin is not None or xmax is not None:
roi = (
xmin if xmin is not None else -np.inf,
xmax if xmax is not None else np.inf,
)
return cls(
data=XYData(x=data["x"], y=data["y"]),
shape=parameters.get("shape", "gaussian"),
baseline=parameters.get("baseline", "none"),
roi=roi,
sensitivity=parameters.get("sensitivity", _DEFAULT_SENSITIVITY),
fwhm_points=parameters.get("fwhm_points"),
auto_fwhm=parameters.get("auto_fwhm", True),
strip_width=parameters.get("strip_width"),
strip_iterations=parameters.get("strip_iterations"),
snip_width=parameters.get("snip_width"),
)