From a5a1c8adecd963be4d3d3096310fa75e957fca30 Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Thu, 1 Oct 2026 10:40:17 +0200 Subject: [PATCH 1/9] add rebin_xafs equivalence tests and timing report --- tests/test_xafs_rebin.py | 302 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 302 insertions(+) create mode 100644 tests/test_xafs_rebin.py diff --git a/tests/test_xafs_rebin.py b/tests/test_xafs_rebin.py new file mode 100644 index 000000000..ae48b6ad0 --- /dev/null +++ b/tests/test_xafs_rebin.py @@ -0,0 +1,302 @@ +#!/usr/bin/env python +"""tests for rebin_xafs: equivalence with the original (loop-based) +implementation, and a timing report on the example data + +run with `pytest -s tests/test_xafs_rebin.py` to see the timing report +""" +import timeit +from pathlib import Path + +import numpy as np +import pytest +from scipy.interpolate import CubicSpline + +from larch import Group +from larch.io import read_ascii +from larch.math import index_of, interp1d +from larch.xafs import pre_edge, rebin_xafs +from larch.xafs.xafsutils import ktoe, etok, ETOK + +DATADIR = Path(__file__).parent.parent / 'examples' / 'xafsdata' + +# file name, numerator, denominator (None: column is mu), +# mu = log(num/den) for transmission +EXAMPLE_FILES = [('cu_rt01.xmu', 'mu', None), + ('fe2o3_rt1.xmu', 'mu', None), + ('cu_10k.xmu', 'mu', None), + ('cu_metal_rt.xdi', 'mutrans', None), + ('fe3c_rt.xdi', 'mutrans', None), + ('ni_metal_rt.xdi', 'mutrans', None), + ('pt_metal_rt.xdi', 'i0', 'itrans'), + ('se_na2so4_rt.xdi', 'i0', 'itrans'), + ('v_foil.xdi', 'i0', 'i1'), + ('znse_zn_xafs.001', 'i0', 'i1'), + ('mn_cpmnco3.dat', 'xmu_trans', None)] + +METHODS = ('boxcar', 'centroid', 'spline') + +# famewoks default parameters +FAMEWOKS_KWS = dict(pre1=None, pre2=-30, pre_step=2, xanes_step=None, + exafs1=15, exafs2=None, exafs_kstep=0.05) + + +def _rebin_xafs_reference(energy, mu, group, e0=None, pre1=None, pre2=-30, + pre_step=2, xanes_step=None, exafs1=15, exafs2=None, + exafs_kstep=0.05, method='boxcar'): + """verbatim copy of the original rebin_xafs (xraylarch 8d6cd41c), + without the group argument parsing""" + if e0 is None: + e0 = getattr(group, 'e0', None) + + if e0 is None: + raise ValueError("need e0") + + emin = min(energy) - e0 + emax = max(energy) - e0 + + if pre1 is None: + pre1 = pre_step*int(emin/pre_step) + + if exafs2 is None: + exafs2 = emax + + # determine xanes step size: + # find mean of energy difference within 10 eV of E0 + if xanes_step is None: + xanes_step = 0.05 * max(1, int(e0 / 1250.0)) # E0/25000, round down to 0.05 + + # clip into data range + pre1 = max(pre1, emin) + pre2 = min(max(pre2, pre1 + abs(pre_step)), emax) + exafs1 = min(max(exafs1, pre2 + abs(xanes_step)), emax) + exafs2 = min(max(exafs2, exafs1 + abs(exafs_kstep) * 20), emax) + + # enforce monotonically increasing + if pre2 <= pre1: + pre2 = min(pre1 + abs(pre_step), emax) + if exafs1 <= pre2: + exafs1 = min(pre2 + abs(xanes_step), emax) + if exafs2 <= exafs1: + exafs2 = min(exafs1 + abs(exafs_kstep) * 20, emax) + + # create new energy array from the 3 segments (pre, xanes, exafs) + # pre: (pre1 -> pre2) with pre_step (in E-space) + # xanes: (pre2 -> exafs1) with xanes_step (in E-space) + # exafs: (exafs1 -> exafs2) with exafs_kstep (in k-space) + en = [] + for start, stop, step, isk in ((pre1, pre2, pre_step, False), + (pre2, exafs1, xanes_step, False), + (exafs1, exafs2, exafs_kstep, True)): + if start == stop: + continue + + if isk: + start = etok(start) + stop = etok(stop) + + npts = 1 + int(0.1 + abs(stop - start) / step) + reg = np.linspace(start, stop, npts) + if isk: + reg = ktoe(reg) + en.extend(e0 + reg[:-1]) + + # find the segment boundaries of the old energy array + bounds = [index_of(energy, e) for e in en] + mu_out = [] + err_out = [] + + j0 = 0 + for i in range(len(en)): + if i == len(en) - 1: + j1 = len(energy) - 1 + else: + j1 = int((bounds[i] + bounds[i+1] + 1)/2.0) + if i == 0 and j0 == 0: + j0 = index_of(energy, en[0]-5) + # if not enough points in segment, do interpolation + # print(i, en[i], j0, j1, len(energy[j0:j1])) + if (j1 - j0) < 3: + jx = j1 + 1 + if (jx - j0) < 3: + jx += 1 + + val = interp1d(energy[j0:jx], mu[j0:jx], en[i]) + err = mu[j0:jx].std() + if np.isnan(val): + j0 = max(0, j0-1) + jx = min(len(energy), jx+1) + val = interp1d(energy[j0:jx], mu[j0:jx], en[i]) + err = mu[j0:jx].std() + else: + if method.startswith('box'): + val = mu[j0:j1].mean() + elif method.startswith('spl'): + val = CubicSpline(energy[j0:j1], mu[j0:j1])(en[i]) + else: + val = (mu[j0:j1]*energy[j0:j1]).mean()/energy[j0:j1].mean() + mu_out.append(val) + if j0 == j1: + err_out.append(np.nan) + else: + err_out.append(mu[j0:j1].std()) + j0 = j1 + + newname = group.__name__ + '_rebinned' + group.rebinned = Group(energy=np.array(en), mu=np.array(mu_out), + delta_mu=np.array(err_out), e0=e0, + __name__=newname) + return + + +def read_example(fname, num, den): + """read an example data file, return a group with energy, mu, e0""" + dat = read_ascii(DATADIR / fname) + energy = dat.energy if hasattr(dat, 'energy') else dat.e + if den is None: + mu = getattr(dat, num) + else: + mu = np.log(getattr(dat, num) / getattr(dat, den)) + group = Group(energy=1.0*energy, mu=1.0*mu, __name__=fname) + pre_edge(group) + return group + + +def synthetic_group(e0, npts, nonuniform=False, seed=0): + """XAS-like spectrum: arctan edge + EXAFS-like oscillation + noise""" + rng = np.random.default_rng(seed) + if nonuniform: + # typical step scan: coarse pre-edge, fine XANES, k-spaced EXAFS + pre = np.arange(e0-200, e0-30, 5.0) + xan = np.arange(e0-30, e0+30, 0.3) + kex = np.arange(etok(30), 16, 0.04) + energy = np.concatenate((pre, xan, e0 + ktoe(kex))) + energy = energy + rng.normal(scale=0.01, size=len(energy)) + energy.sort() + else: + energy = np.linspace(e0-200, e0+1000, npts) + xe = energy - e0 + k = np.sqrt(ETOK*np.where(xe > 0, xe, 0)) + mu = (0.2 + 1.0e-4*xe + np.arctan(xe/2.0)/np.pi + 0.5 + + 0.05*np.sin(2*2.5*k)*np.exp(-0.01*k*k)*(xe > 0) + + rng.normal(scale=2e-3, size=len(energy))) + group = Group(energy=energy, mu=mu, __name__=f'syn_{e0}_{len(energy)}') + pre_edge(group) + return group + + +def check_same(group, **kws): + """run the reference and the new rebin_xafs and compare the results. + if the reference raises an exception, the new version must raise the + same exception type""" + ref = Group(__name__=group.__name__, e0=group.e0) + try: + _rebin_xafs_reference(group.energy, group.mu, ref, **kws) + except Exception as exc: + with pytest.raises(type(exc)): + rebin_xafs(group.energy, group.mu, group=group, **kws) + return + rebin_xafs(group.energy, group.mu, group=group, **kws) + assert_same(ref.rebinned, group.rebinned) + + +def assert_same(ref, new): + assert new.__name__ == ref.__name__ + assert new.e0 == ref.e0 + for attr in ('energy', 'mu', 'delta_mu'): + a, b = getattr(ref, attr), getattr(new, attr) + assert a.shape == b.shape, attr + assert np.array_equal(np.isnan(a), np.isnan(b)), attr + np.testing.assert_allclose(b, a, rtol=1e-12, atol=0, + equal_nan=True, err_msg=attr) + + +@pytest.mark.parametrize('method', METHODS) +@pytest.mark.parametrize('fname,num,den', EXAMPLE_FILES) +def test_rebin_examples(fname, num, den, method): + group = read_example(fname, num, den) + check_same(group, method=method, **FAMEWOKS_KWS) + + +SYNTH_KWS = [FAMEWOKS_KWS, + dict(pre1=-150, pre2=-30, pre_step=5, xanes_step=0.3, + exafs1=15, exafs2=None, exafs_kstep=0.05), + dict(pre1=None, pre2=-15, pre_step=2, xanes_step=None, + exafs1=25, exafs2=None, exafs_kstep=0.05), + dict(pre1=-100, pre2=-20, pre_step=1, xanes_step=0.3, + exafs1=30, exafs2=800, exafs_kstep=0.1)] + + +@pytest.mark.parametrize('method', METHODS) +@pytest.mark.parametrize('kws', SYNTH_KWS) +@pytest.mark.parametrize('npts,nonuniform', [(500, False), (3000, False), + (0, True)]) +@pytest.mark.parametrize('e0', [7112.0, 8979.0, 26711.0]) +def test_rebin_synthetic(e0, npts, nonuniform, kws, method): + group = synthetic_group(e0, npts, nonuniform=nonuniform) + check_same(group, method=method, **kws) + + +@pytest.mark.parametrize('method', METHODS) +def test_rebin_edge_cases(method): + # data starting after e0 + pre1 + grp = synthetic_group(8979.0, 800) + check_same(grp, method=method, pre1=-300) + + # data ending before the EXAFS region + grp = synthetic_group(8979.0, 800) + sel = grp.energy < grp.e0 + 10 + grp = Group(energy=grp.energy[sel], mu=grp.mu[sel], e0=grp.e0, + __name__='short') + check_same(grp, method=method) + + # duplicated energy values + grp = synthetic_group(8979.0, 500) + energy = grp.energy.copy() + energy[100:300:7] = energy[99:299:7] + energy[250] = energy[251] = energy[252] + grp.energy = energy + check_same(grp, method=method) + + # non-monotonic energy (fallback path) + grp = synthetic_group(8979.0, 500) + energy = grp.energy.copy() + energy[200], energy[201] = energy[201], energy[200] + energy[400] = energy[350] + grp.energy = energy + check_same(grp, method=method) + + # NaN in mu + for npts in (500, 3000): + grp = synthetic_group(8979.0, npts) + grp.mu[[10, 150, 151, 220, 260, 400]] = np.nan + # (CubicSpline refuses NaN: the spline method may raise) + check_same(grp, method=method, **FAMEWOKS_KWS) + + +def time_call(func, ncalls, repeat=5): + """median time per call in ms""" + timer = timeit.Timer(func) + return 1000*np.median(timer.repeat(repeat=repeat, number=ncalls))/ncalls + + +def test_rebin_timing_report(): + """report time per call of the reference and the new rebin_xafs""" + groups = [read_example(*args) for args in EXAMPLE_FILES] + groups.append(synthetic_group(26711.0, 3000)) + lines = [f"{'data':22s} {'npts':>5s} {'method':9s} " + f"{'ref (ms)':>9s} {'new (ms)':>9s} {'speedup':>8s}"] + for group in groups: + for method in METHODS: + ref = Group(__name__=group.__name__, e0=group.e0) + def run_ref(): + _rebin_xafs_reference(group.energy, group.mu, ref, + method=method, **FAMEWOKS_KWS) + def run_new(): + rebin_xafs(group.energy, group.mu, group=group, + method=method, **FAMEWOKS_KWS) + t_ref = time_call(run_ref, 2, repeat=3) + t_new = time_call(run_new, 20) + lines.append(f"{group.__name__:22s} {len(group.energy):5d} " + f"{method:9s} {t_ref:9.3f} {t_new:9.3f} " + f"{t_ref/t_new:8.1f}") + print('\nrebin_xafs timing (median per call)\n' + '\n'.join(lines)) From b42074aa4ce7796d33e40cfe138e9dde97f1176b Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Thu, 1 Oct 2026 10:40:17 +0200 Subject: [PATCH 2/9] vectorize rebin_xafs binning, keeping results identical --- larch/xafs/rebin_xafs.py | 111 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 104 insertions(+), 7 deletions(-) diff --git a/larch/xafs/rebin_xafs.py b/larch/xafs/rebin_xafs.py index ab6b9b86f..cdf32a4dd 100644 --- a/larch/xafs/rebin_xafs.py +++ b/larch/xafs/rebin_xafs.py @@ -168,6 +168,27 @@ def rebin_xafs(energy, mu=None, group=None, e0=None, pre1=None, pre2=-30, reg = ktoe(reg) en.extend(e0 + reg[:-1]) + energy = np.asarray(energy) + mu = np.asarray(mu) + en = np.asarray(en, dtype=float) + mu_out = None + if len(en) > 0 and len(energy) > 1 and np.all(np.diff(energy) >= 0): + mu_out, err_out = _rebin_sorted(energy, mu, en, method) + if mu_out is None: + mu_out, err_out = _rebin_loop(energy, mu, en, method) + + newname = group.__name__ + '_rebinned' + group.rebinned = Group(energy=np.array(en), mu=np.array(mu_out), + delta_mu=np.array(err_out), e0=e0, + __name__=newname) + return + + +def _rebin_loop(energy, mu, en, method): + """rebin energy, mu onto en, one bin at a time (works for unsorted energy) + + returns lists of mu and delta_mu values + """ # find the segment boundaries of the old energy array bounds = [index_of(energy, e) for e in en] mu_out = [] @@ -189,12 +210,10 @@ def rebin_xafs(energy, mu=None, group=None, e0=None, pre1=None, pre2=-30, jx += 1 val = interp1d(energy[j0:jx], mu[j0:jx], en[i]) - err = mu[j0:jx].std() if np.isnan(val): j0 = max(0, j0-1) jx = min(len(energy), jx+1) val = interp1d(energy[j0:jx], mu[j0:jx], en[i]) - err = mu[j0:jx].std() else: if method.startswith('box'): val = mu[j0:j1].mean() @@ -209,8 +228,86 @@ def rebin_xafs(energy, mu=None, group=None, e0=None, pre1=None, pre2=-30, err_out.append(mu[j0:j1].std()) j0 = j1 - newname = group.__name__ + '_rebinned' - group.rebinned = Group(energy=np.array(en), mu=np.array(mu_out), - delta_mu=np.array(err_out), e0=e0, - __name__=newname) - return + return mu_out, err_out + + +def _interp_bin(energy, mu, en, i, j0, j1): + """linear interpolation for bin i, segment energy[j0:j1] with + fewer than 3 points, as in _rebin_loop + + returns value and (possibly widened) j0 + """ + jx = j1 + 1 + if (jx - j0) < 3: + jx += 1 + val = interp1d(energy[j0:jx], mu[j0:jx], en[i]) + if np.isnan(val): + j0 = max(0, j0-1) + jx = min(len(energy), jx+1) + val = interp1d(energy[j0:jx], mu[j0:jx], en[i]) + return val, j0 + + +def _rebin_sorted(energy, mu, en, method): + """vectorized version of _rebin_loop, for monotonically increasing energy + + gives the same results as _rebin_loop, or (None, None) if + the bins cannot be found this way + """ + npts, nbins = len(energy), len(en) + # bin boundaries, as index_of(energy, en) for sorted energy + bounds = np.clip(np.searchsorted(energy, en, side='right') - 1, 0, None) + j1 = np.empty(nbins, dtype=int) + j1[:-1] = (bounds[:-1] + bounds[1:] + 1) // 2 + j1[-1] = npts - 1 + j0 = np.empty(nbins, dtype=int) + j0[1:] = j1[:-1] + j0[0] = max(0, np.searchsorted(energy, en[0]-5, side='right') - 1) + if not (np.all(np.isfinite(en)) and np.all(j1 >= j0)): + return None, None + + mu_out = np.full(nbins, np.nan) + err_out = np.full(nbins, np.nan) + + # mean and std of mu in all non-empty segments + full = j1 > j0 + if full.any(): + jstart = j0[full] + jend = j1[full][-1] + nseg = (j1 - j0)[full] + offsets = jstart - jstart[0] + mu_seg = mu[jstart[0]:jend] + mean = np.add.reduceat(mu_seg, offsets) / nseg + dev = mu_seg - np.repeat(mean, nseg) + err_out[full] = np.sqrt(np.add.reduceat(dev*dev, offsets) / nseg) + if method.startswith('box'): + mu_out[full] = mean + elif not method.startswith('spl'): + en_seg = energy[jstart[0]:jend] + mu_out[full] = ((np.add.reduceat(mu_seg*en_seg, offsets) / nseg) / + (np.add.reduceat(en_seg, offsets) / nseg)) + + # segments with fewer than 3 points: linear interpolation. + # the interpolation window is energy[j0:jx], use np.interp on the + # full arrays where the bracketing points are inside the window. + small = (j1 - j0) < 3 + jx = j1 + 1 + jx = np.minimum(jx + ((jx - j0) < 3), npts) + jlo = np.searchsorted(energy, en, side='right') - 1 + val = np.interp(en, energy, mu) + easy = small & (jlo >= j0) & (jlo <= jx - 2) & ~np.isnan(val) + mu_out[easy] = val[easy] + + # remaining bins, in order: spline and other interpolations + if method.startswith('spl'): + todo = ~easy + else: + todo = small & ~easy + for i in np.where(todo)[0]: + if small[i]: + mu_out[i], j0w = _interp_bin(energy, mu, en, i, j0[i], j1[i]) + if j0w != j0[i]: + err_out[i] = np.nan if j0w == j1[i] else mu[j0w:j1[i]].std() + else: + mu_out[i] = CubicSpline(energy[j0[i]:j1[i]], mu[j0[i]:j1[i]])(en[i]) + return mu_out, err_out From d3e66db05e978245c659a342f586097c8156e71d Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Fri, 2 Oct 2026 10:02:37 +0200 Subject: [PATCH 3/9] rename to TEST_KWS for better reading --- tests/test_xafs_rebin.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_xafs_rebin.py b/tests/test_xafs_rebin.py index ae48b6ad0..441f941bb 100644 --- a/tests/test_xafs_rebin.py +++ b/tests/test_xafs_rebin.py @@ -36,7 +36,7 @@ METHODS = ('boxcar', 'centroid', 'spline') # famewoks default parameters -FAMEWOKS_KWS = dict(pre1=None, pre2=-30, pre_step=2, xanes_step=None, +TEST_KWS = dict(pre1=None, pre2=-30, pre_step=2, xanes_step=None, exafs1=15, exafs2=None, exafs_kstep=0.05) @@ -214,10 +214,10 @@ def assert_same(ref, new): @pytest.mark.parametrize('fname,num,den', EXAMPLE_FILES) def test_rebin_examples(fname, num, den, method): group = read_example(fname, num, den) - check_same(group, method=method, **FAMEWOKS_KWS) + check_same(group, method=method, **TEST_KWS) -SYNTH_KWS = [FAMEWOKS_KWS, +SYNTH_KWS = [TEST_KWS, dict(pre1=-150, pre2=-30, pre_step=5, xanes_step=0.3, exafs1=15, exafs2=None, exafs_kstep=0.05), dict(pre1=None, pre2=-15, pre_step=2, xanes_step=None, @@ -270,7 +270,7 @@ def test_rebin_edge_cases(method): grp = synthetic_group(8979.0, npts) grp.mu[[10, 150, 151, 220, 260, 400]] = np.nan # (CubicSpline refuses NaN: the spline method may raise) - check_same(grp, method=method, **FAMEWOKS_KWS) + check_same(grp, method=method, **TEST_KWS) def time_call(func, ncalls, repeat=5): @@ -290,10 +290,10 @@ def test_rebin_timing_report(): ref = Group(__name__=group.__name__, e0=group.e0) def run_ref(): _rebin_xafs_reference(group.energy, group.mu, ref, - method=method, **FAMEWOKS_KWS) + method=method, **TEST_KWS) def run_new(): rebin_xafs(group.energy, group.mu, group=group, - method=method, **FAMEWOKS_KWS) + method=method, **TEST_KWS) t_ref = time_call(run_ref, 2, repeat=3) t_new = time_call(run_new, 20) lines.append(f"{group.__name__:22s} {len(group.energy):5d} " From aff0d7e7b3d7eeb70a128acfe0c655ec977c5a41 Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Fri, 2 Oct 2026 11:50:21 +0200 Subject: [PATCH 4/9] update dev build --- .pre-commit-config.yaml | 28 ++++++++++++++++++++++++++++ pyproject.toml | 15 +++++++++++++-- 2 files changed, 41 insertions(+), 2 deletions(-) create mode 100644 .pre-commit-config.yaml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 000000000..7a290a4f1 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,28 @@ +# pre-commit hooks, see https://pre-commit.com +# install once with: pre-commit install +# run on all files: pre-commit run --all-files +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v6.0.0 + hooks: + - id: check-yaml + - id: check-toml + - id: check-merge-conflict + - id: check-added-large-files + + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.15.22 + hooks: + - id: ruff-check + + # mypy runs from the active environment (pip install -e ".[dev]"), + # so it sees the installed numpy and scipy-stubs + - repo: local + hooks: + - id: mypy + name: mypy + entry: mypy + language: system + types: [python] + require_serial: true + exclude: ^(sandbox|doc|examples|history)/ diff --git a/pyproject.toml b/pyproject.toml index d89ec90fe..d2c8d9638 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -93,7 +93,9 @@ Documentation = "https://xraypy.github.io/xraylarch/" Tracker = "https://github.com/xraypy/xraylarch/issues" [project.optional-dependencies] -dev = ["build", "twine", "pre-commit", "pytest", "pytest-xvfb;platform_system=='Linux'", "pytest-cov", "coverage"] +dev = ["build", "twine", "pre-commit", "ruff", "mypy", "scipy-stubs", + "pytest", "pytest-xvfb;platform_system=='Linux'", "pytest-cov", "pytest-xdist", + "coverage", "hypothesis"] doc = ["Sphinx", "numpydoc", "sphinx-copybutton", "sphinx-subfigure", "sphinxcontrib-video", "sphinxcontrib-bibtex", "sphinxcontrib-argdoc"] wxgui = ["wxpython>=4.3.0", "wxmplot>=2026.2.2", "wxutils>=2026.3.0", "sitka-spruce>=0.4.1", "darkdetect", "pyobjc-framework-Cocoa; platform_system == 'Darwin'"] @@ -113,4 +115,13 @@ exclude = ['larch/__init__.py', 'larch/version.py', line-length = 95 [tool.ruff.lint] -ignore = ['BLE001', 'I001'] # 'except Exception, order of imports +ignore = ['BLE001', # except Exception + 'I001', # order of imports + 'E265', # block comment should start with '# ' + 'E501'] # line too long + +[tool.mypy] +python_version = "3.11" +ignore_missing_imports = true +follow_imports = 'silent' +exclude = ['^sandbox/', '^doc/', '^dist/', '^build/', '^examples/', '^history/'] From 83e36819bdd7e08b5f2bc1783cfac0be18a6b18a Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Fri, 2 Oct 2026 11:50:48 +0200 Subject: [PATCH 5/9] moved into pyproject.toml --- .flake8 | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 .flake8 diff --git a/.flake8 b/.flake8 deleted file mode 100644 index 7ad6cd81d..000000000 --- a/.flake8 +++ /dev/null @@ -1,2 +0,0 @@ -[flake8] -ignore = E265, E501 From 48bb12095c370b41bd0ee9e199c2ada8cacd9cc6 Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Fri, 2 Oct 2026 11:53:46 +0200 Subject: [PATCH 6/9] dev-related: avoid mypy errors --- larch/symboltable.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/larch/symboltable.py b/larch/symboltable.py index 903ca2955..27125e3f8 100644 --- a/larch/symboltable.py +++ b/larch/symboltable.py @@ -3,6 +3,7 @@ SymbolTable for Larch interpreter ''' import copy +from typing import TYPE_CHECKING, Any import numpy from pyshortcuts import gformat from . import site_config @@ -41,6 +42,12 @@ def __init__(self, name=None, **kws): for key, val in kws.items(): setattr(self, key, val) + if TYPE_CHECKING: + # Group members are set dynamically: tell type checkers that + # any attribute may exist (no effect at runtime) + def __getattr__(self, name: str) -> Any: ... + def __setattr__(self, name: str, value: Any) -> None: ... + def __len__(self): return len(dir(self)) From d96e490572615127eb3b63d1541c247fb19c032a Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Fri, 2 Oct 2026 11:54:10 +0200 Subject: [PATCH 7/9] add type hinting --- larch/xafs/rebin_xafs.py | 41 +++++++++++++++++++++++++++------------- 1 file changed, 28 insertions(+), 13 deletions(-) diff --git a/larch/xafs/rebin_xafs.py b/larch/xafs/rebin_xafs.py index cdf32a4dd..85eca209e 100644 --- a/larch/xafs/rebin_xafs.py +++ b/larch/xafs/rebin_xafs.py @@ -1,4 +1,5 @@ import numpy as np +from numpy.typing import ArrayLike, NDArray from scipy.interpolate import CubicSpline from larch import Group @@ -8,7 +9,9 @@ from .xafsutils import ktoe, etok, TINY_ENERGY @Make_CallArgs(["energy", "mu"]) -def sort_xafs(energy, mu=None, group=None, fix_repeats=True, remove_nans=True, overwrite=True): +def sort_xafs(energy: ArrayLike | Group, mu: ArrayLike | None = None, + group: Group | None = None, fix_repeats: bool = True, + remove_nans: bool = True, overwrite: bool = True) -> None: """sort energy, mu pair of XAFS data so that energy is monotonically increasing Arguments @@ -33,6 +36,8 @@ def sort_xafs(energy, mu=None, group=None, fix_repeats=True, remove_nans=True, o energy, mu, group = parse_group_args(energy, members=('energy', 'mu'), defaults=(mu,), group=group, fcn_name='sort_xafs') + energy = np.asarray(energy) + mu = np.asarray(mu) indices = np.argsort(energy) new_energy = energy[indices] @@ -54,9 +59,12 @@ def sort_xafs(energy, mu=None, group=None, fix_repeats=True, remove_nans=True, o @Make_CallArgs(["energy", "mu"]) -def rebin_xafs(energy, mu=None, group=None, e0=None, pre1=None, pre2=-30, - pre_step=2, xanes_step=None, exafs1=15, exafs2=None, - exafs_kstep=0.05, method='boxcar'): +def rebin_xafs(energy: ArrayLike | Group, mu: ArrayLike | None = None, + group: Group | None = None, e0: float | None = None, + pre1: float | None = None, pre2: float = -30, + pre_step: float = 2, xanes_step: float | None = None, + exafs1: float = 15, exafs2: float | None = None, + exafs_kstep: float = 0.05, method: str = 'boxcar') -> None: """rebin XAFS energy and mu to a 'standard 3 region XAFS scan' Arguments @@ -151,7 +159,7 @@ def rebin_xafs(energy, mu=None, group=None, e0=None, pre1=None, pre2=-30, # pre: (pre1 -> pre2) with pre_step (in E-space) # xanes: (pre2 -> exafs1) with xanes_step (in E-space) # exafs: (exafs1 -> exafs2) with exafs_kstep (in k-space) - en = [] + elist: list[float] = [] for start, stop, step, isk in ((pre1, pre2, pre_step, False), (pre2, exafs1, xanes_step, False), (exafs1, exafs2, exafs_kstep, True)): @@ -166,12 +174,13 @@ def rebin_xafs(energy, mu=None, group=None, e0=None, pre1=None, pre2=-30, reg = np.linspace(start, stop, npts) if isk: reg = ktoe(reg) - en.extend(e0 + reg[:-1]) + elist.extend(e0 + reg[:-1]) energy = np.asarray(energy) mu = np.asarray(mu) - en = np.asarray(en, dtype=float) - mu_out = None + en = np.asarray(elist, dtype=float) + mu_out: NDArray[np.float64] | list[float] | None = None + err_out: NDArray[np.float64] | list[float] | None = None if len(en) > 0 and len(energy) > 1 and np.all(np.diff(energy) >= 0): mu_out, err_out = _rebin_sorted(energy, mu, en, method) if mu_out is None: @@ -184,15 +193,17 @@ def rebin_xafs(energy, mu=None, group=None, e0=None, pre1=None, pre2=-30, return -def _rebin_loop(energy, mu, en, method): +def _rebin_loop(energy: NDArray[np.float64], mu: NDArray[np.float64], + en: NDArray[np.float64], + method: str) -> tuple[list[float], list[float]]: """rebin energy, mu onto en, one bin at a time (works for unsorted energy) returns lists of mu and delta_mu values """ # find the segment boundaries of the old energy array bounds = [index_of(energy, e) for e in en] - mu_out = [] - err_out = [] + mu_out: list[float] = [] + err_out: list[float] = [] j0 = 0 for i in range(len(en)): @@ -231,7 +242,9 @@ def _rebin_loop(energy, mu, en, method): return mu_out, err_out -def _interp_bin(energy, mu, en, i, j0, j1): +def _interp_bin(energy: NDArray[np.float64], mu: NDArray[np.float64], + en: NDArray[np.float64], i: int, j0: int, + j1: int) -> tuple[float, int]: """linear interpolation for bin i, segment energy[j0:j1] with fewer than 3 points, as in _rebin_loop @@ -248,7 +261,9 @@ def _interp_bin(energy, mu, en, i, j0, j1): return val, j0 -def _rebin_sorted(energy, mu, en, method): +def _rebin_sorted(energy: NDArray[np.float64], mu: NDArray[np.float64], + en: NDArray[np.float64], method: str + ) -> tuple[NDArray[np.float64], NDArray[np.float64]] | tuple[None, None]: """vectorized version of _rebin_loop, for monotonically increasing energy gives the same results as _rebin_loop, or (None, None) if From 3633de0a3eb8771c7467b162ae1fc484bb1a25ff Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Fri, 2 Oct 2026 11:54:37 +0200 Subject: [PATCH 8/9] make the test faster --- tests/test_xafs_rebin.py | 43 +++++++++++++++++++++++++++------------- 1 file changed, 29 insertions(+), 14 deletions(-) diff --git a/tests/test_xafs_rebin.py b/tests/test_xafs_rebin.py index 441f941bb..96adde796 100644 --- a/tests/test_xafs_rebin.py +++ b/tests/test_xafs_rebin.py @@ -2,8 +2,10 @@ """tests for rebin_xafs: equivalence with the original (loop-based) implementation, and a timing report on the example data -run with `pytest -s tests/test_xafs_rebin.py` to see the timing report +the timing report is skipped by default, run it with: + LARCH_REBIN_TIMING=1 pytest -s tests/test_xafs_rebin.py -k timing """ +import os import timeit from pathlib import Path @@ -14,7 +16,7 @@ from larch import Group from larch.io import read_ascii from larch.math import index_of, interp1d -from larch.xafs import pre_edge, rebin_xafs +from larch.xafs import find_e0, rebin_xafs from larch.xafs.xafsutils import ktoe, etok, ETOK DATADIR = Path(__file__).parent.parent / 'examples' / 'xafsdata' @@ -126,7 +128,7 @@ def _rebin_xafs_reference(energy, mu, group, e0=None, pre1=None, pre2=-30, j0 = max(0, j0-1) jx = min(len(energy), jx+1) val = interp1d(energy[j0:jx], mu[j0:jx], en[i]) - err = mu[j0:jx].std() + err = mu[j0:jx].std() # noqa: F841 else: if method.startswith('box'): val = mu[j0:j1].mean() @@ -157,7 +159,7 @@ def read_example(fname, num, den): else: mu = np.log(getattr(dat, num) / getattr(dat, den)) group = Group(energy=1.0*energy, mu=1.0*mu, __name__=fname) - pre_edge(group) + group.e0 = find_e0(group.energy, group.mu) return group @@ -180,7 +182,7 @@ def synthetic_group(e0, npts, nonuniform=False, seed=0): + 0.05*np.sin(2*2.5*k)*np.exp(-0.01*k*k)*(xe > 0) + rng.normal(scale=2e-3, size=len(energy))) group = Group(energy=energy, mu=mu, __name__=f'syn_{e0}_{len(energy)}') - pre_edge(group) + group.e0 = find_e0(group.energy, group.mu) return group @@ -210,8 +212,12 @@ def assert_same(ref, new): equal_nan=True, err_msg=attr) -@pytest.mark.parametrize('method', METHODS) -@pytest.mark.parametrize('fname,num,den', EXAMPLE_FILES) +# each example file with one method, in turn +EXAMPLE_CASES = [(*args, METHODS[i % len(METHODS)]) + for i, args in enumerate(EXAMPLE_FILES)] + + +@pytest.mark.parametrize('fname,num,den,method', EXAMPLE_CASES) def test_rebin_examples(fname, num, den, method): group = read_example(fname, num, den) check_same(group, method=method, **TEST_KWS) @@ -226,14 +232,21 @@ def test_rebin_examples(fname, num, den, method): exafs1=30, exafs2=800, exafs_kstep=0.1)] -@pytest.mark.parametrize('method', METHODS) -@pytest.mark.parametrize('kws', SYNTH_KWS) -@pytest.mark.parametrize('npts,nonuniform', [(500, False), (3000, False), - (0, True)]) -@pytest.mark.parametrize('e0', [7112.0, 8979.0, 26711.0]) -def test_rebin_synthetic(e0, npts, nonuniform, kws, method): +# (e0, npts, nonuniform, kws, method): instead of the full product, each +# e0, data spacing, keywords set and method appears at least once +SYNTH_CASES = [(7112.0, 0, True, 3, 'boxcar'), + (7112.0, 500, False, 1, 'spline'), + (8979.0, 500, False, 0, 'centroid'), + (8979.0, 0, True, 2, 'spline'), + (8979.0, 3000, False, 3, 'boxcar'), + (26711.0, 3000, False, 2, 'spline'), + (26711.0, 0, True, 1, 'centroid')] + + +@pytest.mark.parametrize('e0,npts,nonuniform,ikws,method', SYNTH_CASES) +def test_rebin_synthetic(e0, npts, nonuniform, ikws, method): group = synthetic_group(e0, npts, nonuniform=nonuniform) - check_same(group, method=method, **kws) + check_same(group, method=method, **SYNTH_KWS[ikws]) @pytest.mark.parametrize('method', METHODS) @@ -279,6 +292,8 @@ def time_call(func, ncalls, repeat=5): return 1000*np.median(timer.repeat(repeat=repeat, number=ncalls))/ncalls +@pytest.mark.skipif(not os.environ.get('LARCH_REBIN_TIMING'), + reason='benchmark: set LARCH_REBIN_TIMING=1 to run') def test_rebin_timing_report(): """report time per call of the reference and the new rebin_xafs""" groups = [read_example(*args) for args in EXAMPLE_FILES] From 34be3ef37fd707416ae5152f39fe32d0c533ac8b Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Fri, 2 Oct 2026 11:54:49 +0200 Subject: [PATCH 9/9] add new tests --- tests/run_tests.sh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/run_tests.sh b/tests/run_tests.sh index 6c4931c1d..0e01e5e99 100644 --- a/tests/run_tests.sh +++ b/tests/run_tests.sh @@ -4,5 +4,6 @@ pytest test_athena_addgroup.py test_basic_processing.py \ test_larchexamples_xafs.py test_larchexamples_xray.py \ test_math_deglitch.py test_math_utils.py test_plot_rixsdata_import.py \ test_read_beamlinedata.py \ - test_read_xafsdata.py test_symbol_callbacks.py test_xas_data_source.py \ + test_read_xafsdata.py test_symbol_callbacks.py test_xafs_rebin.py \ + test_xas_data_source.py \ test_xrf_roi.py