"""silx fitting engine helpers shared by the peak-fit backend and the widget.
This module wraps :mod:`silx.math.fit` (the engine behind
:class:`silx.gui.fit.FitWidget`) so the headless ewoks task and the interactive
Orange widget configure the *same* fit:
* :func:`make_fit_manager` builds a :class:`~silx.math.fit.fitmanager.FitManager`
carrying only the peak theories we expose (so every selectable shape yields a
meaningful descriptor) and, optionally, only the supported backgrounds.
* :func:`descriptors_from_fit_results` turns silx's raw ``fit_results`` into the
physically meaningful per-peak descriptors (center, height, area, FWHM and the
pseudo-Voigt mixing) the task and widget report.
* :func:`background_curve` evaluates the fitted background over the manager's
x-grid, splitting the total model into background and peak contributions.
* :func:`strip_theory`/:func:`snip_theory` are drop-in replacements for silx's
stock Strip/Snip peak-stripping backgrounds, whose process-wide caches crash
when the data length changes between fits (a region-of-interest change) and
tear under concurrent fits; :func:`make_fit_manager` installs them on every
manager.
silx peak theories report *Height* (the true peak maximum) directly; the
integrated *area* is derived from the height and width(s) with the closed-form
expressions in :func:`area_from_height` -- numerically integrating the model
curve would truncate the heavy Lorentzian/Voigt tails.
"""
from __future__ import annotations
from typing import TYPE_CHECKING
import numpy as np
from silx.math.fit import fittheories
from silx.math.fit.bgtheories import CONFIG, configure, estimate_snip, estimate_strip
from silx.math.fit.filters import savitsky_golay, snip1d, strip
from silx.math.fit.fitmanager import FitManager
from silx.math.fit.fittheory import FitTheory
if TYPE_CHECKING:
from ewoksxas.fit.data import ArrayF64
# silx peak theories we expose, keyed so a fit always yields the descriptors we
# report. Height-based theories (not the "Area ..." variants) so the height is
# read directly and the area is the single derived quantity.
CURATED_THEORIES: tuple[str, ...] = (
"Gaussians",
"Lorentz",
"Pseudo-Voigt Line",
"Split Gaussian",
"Split Lorentz",
"Split Pseudo-Voigt",
)
# silx background theories we expose (the polynomial backgrounds are dropped:
# they rarely make sense for spectroscopic peaks and would clutter the GUI).
# Strip and Snip are silx's iterative peak-stripping backgrounds, suited to
# multi-peak spectra.
SUPPORTED_BACKGROUNDS: tuple[str, ...] = (
"No Background",
"Constant",
"Linear",
"Strip",
"Snip",
)
# Area-under-a-unit-height profile of full width at half maximum 1, per family:
# Gaussian integral of exp(-4 ln2 x^2 / w^2) = w * sqrt(pi / (4 ln2))
# Lorentzian integral of 1 / (1 + 4 x^2 / w^2) = w * pi / 2
_GAUSSIAN_AREA_COEFF = float(np.sqrt(np.pi / (4.0 * np.log(2.0))))
_LORENTZIAN_AREA_COEFF = float(np.pi / 2.0)
[docs]
def make_fit_manager(curated_backgrounds: bool = False) -> FitManager:
"""Return a FitManager carrying the curated peak theories.
The standard silx backgrounds (No Background, Constant, Linear, Strip, Snip,
polynomials) are loaded by :class:`FitManager` itself; this registers only
the peak theories we map to descriptors, and swaps Strip/Snip for the safely
cached replacements (:func:`strip_theory`/:func:`snip_theory`). Shared by the
headless backend and the interactive FitWidget so both run the same fit.
Args:
curated_backgrounds: Keep only :data:`SUPPORTED_BACKGROUNDS`, dropping
the polynomial backgrounds (used by the widget to keep the GUI
selection meaningful for spectroscopy).
"""
manager = FitManager()
# Replace silx's stock Strip/Snip backgrounds: their module-global caches
# crash when the data length changes between fits (a region-of-interest
# change) and are shared across managers and threads. Same names, so the
# FitWidget's background dialog and the combo order are unchanged (see the
# strip_theory/snip_theory factories below).
manager.addbgtheory("Strip", theory=strip_theory())
manager.addbgtheory("Snip", theory=snip_theory())
for name in CURATED_THEORIES:
manager.addtheory(name, fittheories.THEORY[name])
if curated_backgrounds:
for name in list(manager.bgtheories):
if name not in SUPPORTED_BACKGROUNDS:
del manager.bgtheories[name]
manager.settheory(CURATED_THEORIES[0])
manager.setbackground("No Background")
return manager
[docs]
def area_from_height( # noqa: PLR0913
family: str,
height: float,
*,
fwhm: float | None = None,
fwhm_low: float | None = None,
fwhm_high: float | None = None,
eta: float = 0.0,
) -> float:
"""Return the analytic integrated area of a peak from its height and width.
For asymmetric (split) peaks each half contributes half of a symmetric peak
of the corresponding side width. Pseudo-Voigt areas blend the Gaussian and
Lorentzian coefficients by the mixing fraction ``eta``.
Args:
family: Profile family: ``"gaussian"``, ``"lorentzian"`` or ``"voigt"``.
height: Peak maximum (silx ``Height``).
fwhm: Full width at half maximum (symmetric peaks).
fwhm_low: Low-side FWHM (asymmetric peaks).
fwhm_high: High-side FWHM (asymmetric peaks).
eta: Pseudo-Voigt Lorentzian fraction (Voigt family only).
"""
if family == "lorentzian":
coeff = _LORENTZIAN_AREA_COEFF
elif family == "voigt":
coeff = eta * _LORENTZIAN_AREA_COEFF + (1.0 - eta) * _GAUSSIAN_AREA_COEFF
else:
coeff = _GAUSSIAN_AREA_COEFF
if fwhm is not None:
return float(height * fwhm * coeff)
return float(0.5 * height * (fwhm_low + fwhm_high) * coeff)
def _value_sigma(results: dict[str, dict], name: str) -> tuple[float, float]:
"""Return the (value, sigma) of a fitted parameter; sigma NaN if missing."""
entry = results[name]
value = float(entry["fitresult"])
sigma = entry.get("sigma")
sigma_f = float("nan") if sigma is None else float(sigma)
if not np.isfinite(sigma_f):
sigma_f = float("nan")
return value, sigma_f
def _relative(value: float, sigma: float) -> float:
"""Relative uncertainty ``sigma / value`` (0 when undefined)."""
if not np.isfinite(sigma) or value == 0.0:
return 0.0
return sigma / value
[docs]
def descriptors_from_fit_results(
fit_results: list[dict], family: str, asymmetric: bool
) -> list[dict[str, float]]:
"""Return the per-peak descriptors from silx ``fit_results``.
Each peak (silx group ``>= 1``; group ``0`` holds the background parameters)
yields a descriptor dict with ``center``, ``center_stderr``, ``height``,
``area``, ``area_stderr`` and the shape-dependent width(s): a single ``fwhm``
for symmetric peaks, per-side ``fwhm_low``/``fwhm_high`` for asymmetric
ones, plus ``eta`` for the Voigt family. The peaks are returned in silx
order (left to right ordering is done by the caller).
Args:
fit_results: silx ``FitManager.fit_results`` after ``runfit``.
family: Profile family of the fitted shape.
asymmetric: Whether the shape reports per-side widths.
"""
by_name = {entry["name"]: entry for entry in fit_results}
groups = sorted({entry["group"] for entry in fit_results if entry["group"] >= 1})
peaks: list[dict[str, float]] = []
for group in groups:
center, center_sigma = _value_sigma(by_name, f"Position{group}")
height, height_sigma = _value_sigma(by_name, f"Height{group}")
eta = _value_sigma(by_name, f"Eta{group}")[0] if family == "voigt" else 0.0
descriptors: dict[str, float] = {
"center": center,
"center_stderr": center_sigma,
"height": height,
}
if asymmetric:
fwhm_low, low_sigma = _value_sigma(by_name, f"LowFWHM{group}")
fwhm_high, high_sigma = _value_sigma(by_name, f"HighFWHM{group}")
descriptors["fwhm_low"] = fwhm_low
descriptors["fwhm_high"] = fwhm_high
area = area_from_height(
family, height, fwhm_low=fwhm_low, fwhm_high=fwhm_high, eta=eta
)
width_rel = _relative(fwhm_low + fwhm_high, np.hypot(low_sigma, high_sigma))
else:
fwhm, fwhm_sigma = _value_sigma(by_name, f"FWHM{group}")
descriptors["fwhm"] = fwhm
area = area_from_height(family, height, fwhm=fwhm, eta=eta)
width_rel = _relative(fwhm, fwhm_sigma)
descriptors["area"] = area
area_rel = np.hypot(_relative(height, height_sigma), width_rel)
descriptors["area_stderr"] = (
float(area * area_rel) if area_rel else float("nan")
)
if family == "voigt":
descriptors["eta"] = eta
peaks.append(descriptors)
return peaks
[docs]
def background_curve(manager: FitManager) -> ArrayF64:
"""Evaluate the fitted background of a manager over its x-grid.
silx background functions take ``(x, ydata, *bg_params)``; the group-0
entries of ``fit_results`` hold the background parameters (fitted for
Constant/Linear, fixed filter settings for Strip/Snip).
"""
bgfun = manager.bgtheories[manager.selectedbg].function
bg_values = [
entry["fitresult"] for entry in manager.fit_results if entry["group"] == 0
]
return np.asarray(bgfun(manager.xdata, manager.ydata, *bg_values), dtype=np.float64)
# ---------------------------------------------------------------------------
# Safe Strip/Snip background theories (installed by ``make_fit_manager``).
#
# silx's stock Strip/Snip background functions cache the last evaluation in
# module-level globals shared by every FitManager in the process. Two defects
# follow (still present upstream as of June 2026):
#
# * ``snip_bg`` compares the cached y-array with ``numpy.sum(old == y) == len(y)``,
# which raises a broadcast ``ValueError`` whenever the data length changes
# between fits with an unchanged width -- exactly what re-fitting after a
# region-of-interest change does (``FitManager.setdata`` slices to xmin/xmax).
# * Both functions publish their cache in two steps (data first, background
# later), so concurrent fits -- the widget's interactive FitManager in the GUI
# thread and the batch task's managers in a worker thread -- can read a torn,
# wrong-length cache.
#
# The factories below build behaviour-identical theories whose cache lives in a
# closure (one per FitManager, no cross-manager state), compares data with the
# shape-safe ``numpy.array_equal`` and is published as a single tuple only after
# the background is computed. ``Strip`` additionally keys its cache on
# ``StripThresholdFactor``, which the stock cache ignores. The silx parameter
# names, ``estimate_strip``/``estimate_snip``, ``configure`` and the shared
# ``CONFIG`` dict are reused, so the FitWidget's background dialog and the fit
# configuration flow are unchanged. A cache is required at all: ``leastsq``
# evaluates the background function hundreds of times per fit with identical
# arguments, and the Strip filter alone defaults to 5000 iterations.
# ---------------------------------------------------------------------------
# One cached evaluation: (copy of the y data, settings key, background).
# Published as a single tuple so a reader never observes a half-updated cache.
_Cache = tuple[np.ndarray, tuple, np.ndarray] | None
def _anchor_indices(
x: ArrayF64, anchors_flag: bool, anchors_list: tuple[float, ...]
) -> list[int] | None:
"""Return the indices in ``x`` of the configured anchor abscissae.
Vendored from silx's private ``bgtheories._convert_anchors_to_indices`` so
the replacement theories do not depend on silx internals: each anchor maps
to the first index where ``x >= anchor`` (anchors at or before ``x[0]`` are
skipped). Returns None when anchors are disabled or none falls inside.
"""
if not anchors_flag or not anchors_list:
return None
indices = []
for anchor_x in anchors_list:
if anchor_x <= x[0]:
continue
candidates = np.nonzero(x >= anchor_x)[0]
if len(candidates):
indices.append(int(candidates.min()))
return indices or None
[docs]
def strip_theory() -> FitTheory:
"""Return a fresh Strip background theory with a per-manager cache."""
cache: _Cache = None
def strip_bg(x: ArrayF64, y0: ArrayF64, width: float, niter: float) -> ArrayF64:
nonlocal cache
# Snapshot the shared CONFIG once so a concurrent edit (e.g. from the
# FitWidget's background dialog) cannot tear a single evaluation.
smoothing_flag = bool(CONFIG["SmoothingFlag"])
smoothing_width = int(CONFIG["SmoothingWidth"])
anchors_flag = bool(CONFIG["AnchorsFlag"])
anchors_list = tuple(CONFIG["AnchorsList"] or ())
factor = CONFIG["StripThresholdFactor"]
key = (
width,
niter,
smoothing_flag,
smoothing_width,
anchors_flag,
anchors_list,
factor,
)
snapshot = cache
if (
snapshot is not None
and snapshot[1] == key
and np.array_equal(snapshot[0], y0)
):
return snapshot[2]
y1 = savitsky_golay(y0, smoothing_width) if smoothing_flag else y0
background = strip(
y1,
w=width,
niterations=niter,
factor=factor,
anchors=_anchor_indices(x, anchors_flag, anchors_list),
)
cache = (np.array(y0, copy=True), key, background)
return background
return FitTheory(
description="Compute background using a strip filter\n"
"Parameters 'StripWidth', 'StripIterations'",
function=strip_bg,
parameters=["StripWidth", "StripIterations"],
estimate=estimate_strip,
configure=configure,
is_background=True,
)
[docs]
def snip_theory() -> FitTheory:
"""Return a fresh Snip background theory with a per-manager cache."""
cache: _Cache = None
def snip_bg(x: ArrayF64, y0: ArrayF64, width: float) -> ArrayF64:
nonlocal cache
smoothing_flag = bool(CONFIG["SmoothingFlag"])
smoothing_width = int(CONFIG["SmoothingWidth"])
anchors_flag = bool(CONFIG["AnchorsFlag"])
anchors_list = tuple(CONFIG["AnchorsList"] or ())
key = (width, smoothing_flag, smoothing_width, anchors_flag, anchors_list)
snapshot = cache
if (
snapshot is not None
and snapshot[1] == key
and np.array_equal(snapshot[0], y0)
):
return snapshot[2]
y1 = savitsky_golay(y0, smoothing_width) if smoothing_flag else y0
# The snip filter runs per anchor-delimited segment, as in silx.
anchors = _anchor_indices(x, anchors_flag, anchors_list)
if not anchors:
anchors = [0, len(y1) - 1]
background = np.zeros_like(y1)
previous = 0
for anchor in anchors:
if previous < anchor < len(y1):
background[previous:anchor] = snip1d(y1[previous:anchor], width)
previous = anchor
if previous < len(y1):
background[previous:] = snip1d(y1[previous:], width)
cache = (np.array(y0, copy=True), key, background)
return background
return FitTheory(
description="Compute background using a snip filter\nParameter 'SnipWidth'",
function=snip_bg,
parameters=["SnipWidth"],
estimate=estimate_snip,
configure=configure,
is_background=True,
)