"""Shared fixtures synthesising ground-truth peaks for the peak-fit tests.
Peaks are built from silx's own profile functions (the engine behind the fit),
on a noiseless grid, so the fitted center, area, width(s) and mixing are known
exactly. ``peak_factory`` returns ``(y, truth)`` where ``truth`` holds the
descriptors the fit should recover. An independent lmfit-based Gaussian oracle
(``lmfit_gaussian``) guards against testing silx purely against itself.
"""
import numpy as np
import pytest
from lmfit.models import GaussianModel
from numpy.typing import NDArray
from silx.math.fit import functions
# Area under a unit-height profile of FWHM 1 (see ewoksxas.fit.silx_engine).
_GAUSS = float(np.sqrt(np.pi / (4.0 * np.log(2.0))))
_LORENTZ = float(np.pi / 2.0)
# Ground-truth parameters used to synthesise every test peak.
TRUE_CENTER = 2.0
TRUE_HEIGHT = 8.0
TRUE_FWHM = 4.0
TRUE_LOW_FWHM = 4.0 # asymmetric low-side FWHM
TRUE_HIGH_FWHM = 6.0 # asymmetric high-side FWHM
TRUE_ETA = 0.4 # pseudo-Voigt Lorentzian fraction
def _symmetric_area(family: str, height: float, fwhm: float, eta: float) -> float:
if family == "lorentzian":
coeff = _LORENTZ
elif family == "voigt":
coeff = eta * _LORENTZ + (1.0 - eta) * _GAUSS
else:
coeff = _GAUSS
return height * fwhm * coeff
def _truth(family: str, *, asymmetric: bool, eta: float = 0.0) -> dict:
truth: dict[str, float] = {"center": TRUE_CENTER, "height": TRUE_HEIGHT}
if asymmetric:
truth["fwhm_low"] = TRUE_LOW_FWHM
truth["fwhm_high"] = TRUE_HIGH_FWHM
truth["area"] = 0.5 * (
_symmetric_area(family, TRUE_HEIGHT, TRUE_LOW_FWHM, eta)
+ _symmetric_area(family, TRUE_HEIGHT, TRUE_HIGH_FWHM, eta)
)
else:
truth["fwhm"] = TRUE_FWHM
truth["area"] = _symmetric_area(family, TRUE_HEIGHT, TRUE_FWHM, eta)
if family == "voigt":
truth["eta"] = eta
return truth
[docs]
@pytest.fixture
def x() -> NDArray[np.float64]:
"""A fine, wide x-grid so even broad asymmetric tails decay at the edges."""
return np.linspace(-20.0, 20.0, 1601, dtype=np.float64)
[docs]
@pytest.fixture
def peak_factory(x):
"""Return a factory building a ground-truth peak and its descriptors.
The factory takes a peak-shape label and returns ``(y, truth)`` so each test
can assert the recovered descriptors against the values used to build them.
"""
def _build(shape: str = "gaussian"):
c, h = TRUE_CENTER, TRUE_HEIGHT
fw, lo, hi, eta = TRUE_FWHM, TRUE_LOW_FWHM, TRUE_HIGH_FWHM, TRUE_ETA
if shape == "gaussian":
y = functions.sum_gauss(x, h, c, fw)
return y, _truth("gaussian", asymmetric=False)
if shape == "lorentzian":
y = functions.sum_lorentz(x, h, c, fw)
return y, _truth("lorentzian", asymmetric=False)
if shape == "voigt":
y = functions.sum_pvoigt(x, h, c, fw, eta)
return y, _truth("voigt", asymmetric=False, eta=eta)
if shape == "skewed gaussian":
y = functions.sum_splitgauss(x, h, c, lo, hi)
return y, _truth("gaussian", asymmetric=True)
if shape == "split lorentzian":
y = functions.sum_splitlorentz(x, h, c, lo, hi)
return y, _truth("lorentzian", asymmetric=True)
if shape == "skewed voigt":
y = functions.sum_splitpvoigt(x, h, c, lo, hi, eta)
return y, _truth("voigt", asymmetric=True, eta=eta)
raise ValueError(f"Unknown shape: {shape!r}.")
return _build
[docs]
@pytest.fixture
def multi_peak_factory(x):
"""Return a factory building a sum of well-separated Gaussian peaks.
The factory takes ``(center, height)`` pairs and returns the summed
noiseless signal, so a test knows how many peaks and which centers/heights
it should recover.
"""
def _build(centers_heights, fwhm: float = TRUE_FWHM):
y = np.zeros_like(x)
for center, height in centers_heights:
y = y + functions.sum_gauss(x, height, center, fwhm)
return y
return _build
[docs]
@pytest.fixture
def lmfit_gaussian(x):
"""An independent (lmfit) Gaussian oracle: ``(y, center, area)``.
lmfit's Gaussian is mathematically identical to silx's, so recovering it
cross-checks the silx-backed fit against a separate implementation.
"""
model = GaussianModel()
params = model.make_params(amplitude=20.0, center=TRUE_CENTER, sigma=2.0)
return model.eval(params, x=x), TRUE_CENTER, 20.0
[docs]
@pytest.fixture
def lca_x() -> NDArray[np.float64]:
return np.linspace(0, np.pi, 100, dtype=np.float64)
[docs]
@pytest.fixture
def component_list(lca_x) -> NDArray[np.float64]:
x = lca_x
comp_1 = np.sin(x)
comp_2 = np.cos(x)
comp_3 = (comp_1 * comp_2) * 2
comp_4 = (comp_3 * comp_1) * (4 / 3)
comp_5 = (comp_3 * comp_2) * (4 / 3)
# return np.array([comp_1, comp_2, comp_3, comp_4, comp_5])
return np.array(
[
comp_1,
comp_2,
comp_3,
comp_4,
comp_5,
x,
np.pi - x,
np.log(x + 0.01),
np.exp(x),
np.log10(x + 0.01),
1 / (x + 0.01),
np.tan(x),
np.sqrt(x),
np.power(x, 2),
]
)
[docs]
@pytest.fixture
def sample(lca_x, component_list) -> dict[str, NDArray[np.float64]]:
weights = np.array([0.1, 0.2, 0.3, 0.35, 0.05])
components = component_list[:5]
y = np.dot(weights, components)
return {"x": lca_x, "y": y}
[docs]
@pytest.fixture
def components(excess_components) -> dict[str, dict[str, NDArray[np.float64]]]:
return {
key: value
for key, value in excess_components.items()
if int(key.split(" ")[-1]) in (1, 2, 3, 4, 5)
}
[docs]
@pytest.fixture
def excess_components(
lca_x, component_list
) -> dict[str, dict[str, NDArray[np.float64]]]:
return {
f"Component {i + 1}": {"x": lca_x, "y": arr}
for i, arr in enumerate(component_list)
}
[docs]
@pytest.fixture
def rand_components(lca_x, component_list) -> dict[str, dict[str, NDArray[np.float64]]]:
comps = component_list.copy()
rng = np.random.default_rng(seed=1337)
rng.shuffle(comps)
return {f"Component {i + 1}": {"x": lca_x, "y": arr} for i, arr in enumerate(comps)}