"""Tests for the PeakFitting task.
The task wraps ewoksxas.fit.peak_fit.PeakModel: it fits every spectrum in the
input table and returns the fitted model curves with the peak descriptors as
metadata. Spectra are exact (noiseless) peaks on a known baseline, so the
recovered center/area and fit quality (r2 ~ 1) are known in advance. Gaussians
are generated with lmfit (an independent oracle, identical maths to silx);
shapes that differ between the two libraries (pseudo-Voigt) use silx's own
profile functions.
"""
import numpy as np
import pytest
from ewoksorange.tests.utils import execute_task
from lmfit.models import GaussianModel
from silx.math.fit import functions
from ewoksxas.converters.orange import Converter, VarType
from ewoksxas.fit.data import XYData
from ewoksxas.fit.peak_fit import PeakModel
from ewoksxas.tasks.peak_fit import PeakFitting, descriptor_name
# The default (Gaussian) shape's metas: a single fwhm, no per-side widths/eta.
DESCRIPTOR_METAS = (
"center",
"center_stderr",
"height",
"area",
"area_stderr",
"fwhm",
"r2",
"rmse",
)
[docs]
@pytest.fixture
def x():
return np.linspace(-15.0, 15.0, 801, dtype=np.float64)
def _gaussian(x, center, amplitude=20.0, sigma=2.0):
model = GaussianModel()
params = model.make_params(amplitude=amplitude, center=center, sigma=sigma)
return model.eval(params, x=x)
def _sample_table(x, spectra, names=None):
converter = Converter().add_features(x, np.atleast_2d(spectra))
if names is not None:
converter.add_meta("Scan Name", np.asarray(names), var_type=VarType.TEXT)
return converter.to_table()
def _run(table, **parameters):
inputs = {"Data": table, "parameters": parameters}
return execute_task(PeakFitting, inputs=inputs)["Data"]
def _meta(table, name):
return np.asarray(Converter.from_table(table).get_meta_values(name), dtype=float)
[docs]
def test_recovers_peaks_and_reconstructs_spectra(x):
centers = [1.0, 3.0, -2.0]
spectra = np.stack([_gaussian(x, c) + 5.0 for c in centers])
table = _sample_table(x, spectra, [f"S{i}" for i in range(3)])
output = _run(table, shape="gaussian", baseline="offset")
np.testing.assert_allclose(_meta(output, "center"), centers, atol=1e-3)
np.testing.assert_allclose(_meta(output, "area"), 20.0, rtol=1e-3)
np.testing.assert_allclose(_meta(output, "r2"), 1.0, atol=1e-6)
# The model curves reconstruct the input spectra on the same x grid.
model_x, model_y = Converter.from_table(output).features
np.testing.assert_allclose(model_x, x, atol=1e-9)
np.testing.assert_allclose(model_y, spectra, atol=1e-2)
[docs]
def test_outputs_match_direct_fit(x):
spectrum = _gaussian(x, 1.5) + (0.3 * x + 4.0)
table = _sample_table(x, spectrum)
output = _run(table, shape="gaussian", baseline="linear")
fit = PeakModel(XYData(x=x, y=spectrum), shape="gaussian", baseline="linear").fit()
np.testing.assert_allclose(_meta(output, "center")[0], fit.center, atol=1e-9)
np.testing.assert_allclose(_meta(output, "area")[0], fit.area, atol=1e-9)
_, model_y = Converter.from_table(output).features
np.testing.assert_allclose(model_y[0], fit.model.y, atol=1e-9)
[docs]
def test_voigt_reports_mixing(x):
# silx fits a pseudo-Voigt and recovers its mixing fraction directly.
spectrum = functions.sum_pvoigt(x, 8.0, 1.0, 4.0, 0.4)
table = _sample_table(x, spectrum)
output = _run(table, shape="voigt")
assert _meta(output, "eta")[0] == pytest.approx(0.4, abs=2e-2)
assert _meta(output, "fwhm")[0] > 0.0
[docs]
def test_roi_restricts_grid(x):
table = _sample_table(x, _gaussian(x, 1.0))
output = _run(table, xmin=-5.0, xmax=9.0)
model_x, _ = Converter.from_table(output).features
assert model_x.min() >= -5.0
assert model_x.max() <= 9.0
assert model_x.size < x.size
np.testing.assert_allclose(_meta(output, "center"), 1.0, atol=1e-3)
[docs]
def test_asymmetric_shape(x):
# Split Lorentzian: low-side FWHM 4, high-side FWHM 6.
spectrum = functions.sum_splitlorentz(x, 8.0, 1.0, 4.0, 6.0)
table = _sample_table(x, spectrum)
output = _run(table, shape="split lorentzian")
np.testing.assert_allclose(_meta(output, "center"), 1.0, atol=1e-2)
np.testing.assert_allclose(_meta(output, "fwhm_low"), 4.0, atol=1e-2)
np.testing.assert_allclose(_meta(output, "fwhm_high"), 6.0, atol=1e-2)
def _two_gaussians(x, c0, c1):
return _gaussian(x, c0, amplitude=20.0) + _gaussian(x, c1, amplitude=15.0)
[docs]
def test_emits_per_peak_columns(x):
# A two-peak spectrum: auto peak search detects both peaks.
table = _sample_table(x, _two_gaussians(x, -6.0, 6.0))
output = _run(table, shape="gaussian")
# Per-peak columns exist and recover the (center-ordered) peaks.
center_0 = _meta(output, descriptor_name("center", 0))
center_1 = _meta(output, descriptor_name("center", 1))
assert center_0[0] == pytest.approx(-6.0, abs=1e-2)
assert center_1[0] == pytest.approx(6.0, abs=1e-2)
# The bare alias mirrors the first (leftmost) peak.
np.testing.assert_allclose(_meta(output, "center"), center_0)
np.testing.assert_allclose(_meta(output, "n_peaks"), 2.0)
[docs]
def test_auto_n_emits_correct_count(x):
table = _sample_table(x, _two_gaussians(x, -6.0, 6.0))
output = _run(table, shape="gaussian") # auto by default
np.testing.assert_allclose(_meta(output, "n_peaks"), 2.0)
centers = sorted([_meta(output, "center_0")[0], _meta(output, "center_1")[0]])
np.testing.assert_allclose(centers, [-6.0, 6.0], atol=1e-2)
[docs]
def test_ragged_peak_counts_are_nan_padded(x):
# One single-peak and one two-peak spectrum -> the table widens to two peaks
# and the single-peak row pads its second peak with NaN.
spectra = np.stack([_gaussian(x, 0.0), _two_gaussians(x, -6.0, 6.0)])
table = _sample_table(x, spectra)
output = _run(table, shape="gaussian") # auto detects 1 then 2 peaks
n_peaks = _meta(output, "n_peaks")
assert n_peaks[0] == 1.0
assert n_peaks[1] == 2.0
center_1 = _meta(output, descriptor_name("center", 1))
assert np.isnan(center_1[0]) # single-peak row has no second peak
assert center_1[1] == pytest.approx(6.0, abs=1e-2)