From c1032ae58c0c7754b86a8902f6080ec552712daf Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 4 Aug 2026 12:31:20 +0200 Subject: [PATCH 1/6] feat: relative width smoothing procedures --- docs/api-reference/index.md | 1 + src/scippneutron/__init__.py | 1 + src/scippneutron/smoothing.py | 618 ++++++++++++++++++++++++++++++++++ tests/smoothing_test.py | 454 +++++++++++++++++++++++++ 4 files changed, 1074 insertions(+) create mode 100644 src/scippneutron/smoothing.py create mode 100644 tests/smoothing_test.py diff --git a/docs/api-reference/index.md b/docs/api-reference/index.md index 196442178..deb3d33f9 100644 --- a/docs/api-reference/index.md +++ b/docs/api-reference/index.md @@ -65,5 +65,6 @@ and possible confusion of `theta` (from Bragg’s law) with `theta` in spherical logging metadata peaks + smoothing tof ``` diff --git a/src/scippneutron/__init__.py b/src/scippneutron/__init__.py index 4f883682c..717f64b6d 100644 --- a/src/scippneutron/__init__.py +++ b/src/scippneutron/__init__.py @@ -32,6 +32,7 @@ 'data', 'metadata', 'peaks', + 'smoothing', 'tof', ] diff --git a/src/scippneutron/smoothing.py b/src/scippneutron/smoothing.py new file mode 100644 index 000000000..cf1bc51e7 --- /dev/null +++ b/src/scippneutron/smoothing.py @@ -0,0 +1,618 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) + +from __future__ import annotations + +from typing import Any, Protocol, TypeAlias, cast + +import numpy as np +import scipp as sc +from numpy.typing import ArrayLike, NDArray +from scipy.signal import convolve +from scipy.stats import norm, triang, uniform + +__all__ = [ + "smooth_relative_gaussian", + "smooth_relative_kernel", + "smooth_relative_rectangle", + "smooth_relative_triangle", +] + + +class _Distribution(Protocol): + def cdf(self, x: ArrayLike) -> Any: ... + + def ppf(self, probability: ArrayLike) -> Any: ... + + def support(self) -> tuple[Any, Any]: ... + + +_Kernel: TypeAlias = str | _Distribution +_FloatArray: TypeAlias = NDArray[np.float64] +_IntArray: TypeAlias = NDArray[np.int64] +_ScippArray: TypeAlias = sc.DataArray | sc.Variable + +_BUILTIN_KERNELS: dict[str, _Distribution] = { + # Standard Gaussian in relative-coordinate units. + "gaussian": norm(), + "normal": norm(), + # Centered rectangle on [-1, 1]. + # With alpha=a, this has relative support [-a, a]. + "rectangle": uniform(loc=-1.0, scale=2.0), + "rect": uniform(loc=-1.0, scale=2.0), + "box": uniform(loc=-1.0, scale=2.0), + "uniform": uniform(loc=-1.0, scale=2.0), + # Centered triangle on [-1, 1], peak at 0. + # Users can pass triang(c=...) themselves for asymmetric triangles. + "triangle": triang(c=0.5, loc=-1.0, scale=2.0), + "triangular": triang(c=0.5, loc=-1.0, scale=2.0), +} + + +def _as_kernel_distribution(kernel: _Kernel) -> _Distribution: + if isinstance(kernel, str): + try: + return _BUILTIN_KERNELS[kernel.lower()] + except KeyError: + valid = ", ".join(sorted(_BUILTIN_KERNELS)) + raise ValueError( + f"unknown kernel {kernel!r}; expected one of: {valid}" + ) from None + + required = ("cdf", "ppf", "support") + missing = [name for name in required if not hasattr(kernel, name)] + if missing: + raise TypeError( + "kernel must be a scipy.stats distribution-like object with methods " + f"{', '.join(required)}; missing {', '.join(missing)}" + ) + + return kernel + + +def _relative_kernel_weights( + log_spacing: float, + alpha: float, + kernel: _Kernel, + tail: float, + max_offset: int, +) -> tuple[_IntArray, _FloatArray]: + """ + Weights for smoothing on a geometric grid q_i = q0 * exp(i * log_spacing). + + The kernel distribution describes the relative displacement Z: + + q' = q * (1 + alpha * Z) + + Equivalently, + + K(q, q') = 1 / (alpha * q) * f((q' - q) / (alpha * q)) + + where f is the PDF of the supplied distribution. + + """ + if not np.isfinite(alpha) or alpha <= 0: + raise ValueError("alpha must be positive") + if not np.isfinite(log_spacing) or log_spacing <= 0: + raise ValueError("log_spacing must be positive") + if not np.isfinite(tail) or not (0.0 < tail < 1.0): + raise ValueError("tail must be between 0 and 1") + if max_offset < 0: + raise ValueError("max_offset must be non-negative") + + dist = _as_kernel_distribution(kernel) + + # Positive physical domain: + # + # q' > 0 + # q * (1 + alpha * z) > 0 + # z > -1 / alpha + z_domain_min = -1.0 / alpha + + p_domain_min = float(dist.cdf(z_domain_min)) + norm_mass = 1.0 - p_domain_min + + if not np.isfinite(norm_mass) or norm_mass <= 0.0: + raise ValueError("kernel has no positive-domain mass for this alpha") + + try: + support_min, support_max = dist.support() + except TypeError as e: + raise TypeError( + "kernel must be a fully specified distribution. For distributions " + "with shape parameters, pass a frozen distribution such as " + "triang(c=0.5), not triang." + ) from e + + support_min = float(support_min) + support_max = float(support_max) + + # Exact z-support after clipping to q' > 0. + z_left_exact = max(z_domain_min, support_min) + z_right_exact = support_max + + # If the clipped support maps to finite log-space, use it exactly. + # If it touches q' = 0, the log lower bound is -inf, so use a tail cutoff. + has_finite_log_support = ( + np.isfinite(z_left_exact) + and np.isfinite(z_right_exact) + and (1.0 + alpha * z_left_exact > 0.0) + ) + + if has_finite_log_support: + u_left = np.log1p(alpha * z_left_exact) + u_right = np.log1p(alpha * z_right_exact) + else: + probabilities = ( + p_domain_min + np.array([0.5 * tail, 1.0 - 0.5 * tail]) * norm_mass + ) + z_left, z_right = (float(value) for value in dist.ppf(probabilities)) + + u_left = np.log1p(alpha * z_left) + u_right = np.log1p(alpha * z_right) + + if not np.isfinite(u_right): + raise ValueError("right kernel bound is not finite; increase tail") + + # Cells are centered at m*h and span [(m-1/2)h, (m+1/2)h]. + # Include offset zero even for one-sided kernels and clamp the stencil to + # offsets that can contribute to the finite input. + m_min = int(np.clip(np.floor(u_left / log_spacing + 0.5), -max_offset, 0)) + m_max = int(np.clip(np.ceil(u_right / log_spacing - 0.5), 0, max_offset)) + + m = np.arange(m_min, m_max + 1, dtype=np.int64) + + L = (m - 0.5) * log_spacing + U = (m + 0.5) * log_spacing + + # z = (q' - q) / (alpha q) + # = (exp(u) - 1) / alpha + zL = np.expm1(L) / alpha + zU = np.expm1(U) / alpha + + # Exact cell-integrated weights in log-space. + w = (dist.cdf(zU) - dist.cdf(zL)) / norm_mass + w = np.maximum(w, 0.0) + + nonzero = np.flatnonzero(w > 0.0) + if nonzero.size == 0: + # The distribution has no mass within reach of the finite input. + # Returning zero weights lets the caller mark those values as NaN. + return m, w + + # Trim zero-only ends, but preserve offset zero for convolution alignment. + zero = -m_min + first = min(nonzero[0], zero) + last = max(nonzero[-1], zero) + 1 + + # Renormalize after finite tail truncation. + trimmed_weights = w[first:last] + return m[first:last], trimmed_weights / trimmed_weights.sum() + + +def _valid_weight_sums(n: int, m: _IntArray, w: _FloatArray) -> _FloatArray: + """ + Boundary normalization. + + Equivalent to + + np.convolve(np.ones(n), w[::-1], mode="full")[start:start+n] + + but O(n), not another full convolution. + """ + i = np.arange(n, dtype=np.int64) + + m_min = int(m[0]) + m_max = int(m[-1]) + + lower = np.maximum(m_min, -i) + upper = np.minimum(m_max, n - 1 - i) + + cumsum = np.empty(w.size + 1, dtype=float) + cumsum[0] = 0.0 + np.cumsum(w, out=cumsum[1:]) + + return cast( + _FloatArray, + cumsum[upper - m_min + 1] - cumsum[lower - m_min], + ) + + +def _smooth_relative_kernel_on_geomgrid( + y: ArrayLike, + log_spacing: float, + alpha: float, + kernel: _Kernel, + tail: float, +) -> _FloatArray: + y = np.asarray(y, dtype=float) + + if y.ndim != 1: + raise ValueError("y must be one-dimensional") + if not np.isfinite(alpha) or alpha < 0: + raise ValueError("alpha must be non-negative") + if not np.isfinite(log_spacing) or log_spacing <= 0: + raise ValueError("log_spacing must be positive") + if not np.isfinite(tail) or not (0.0 < tail < 1.0): + raise ValueError("tail must be between 0 and 1") + if y.size == 0 or alpha == 0: + return y.copy() + + m, w = _relative_kernel_weights( + log_spacing=log_spacing, + alpha=alpha, + kernel=kernel, + tail=tail, + max_offset=y.size - 1, + ) + + # Desired operation: + # + # out[i] = sum_m w[m] * y[i + m] + # + # scipy convolution reverses the second argument, hence w[::-1]. + full = cast(_FloatArray, convolve(y, w[::-1], mode="full")) + + start = int(m[-1]) + numerator = full[start : start + y.size] + + denom = _valid_weight_sums(y.size, m, w) + + out = np.full_like(numerator, np.nan) + np.divide( + numerator, + denom, + out=out, + where=denom > 0.0, + ) + return out + + +def _smooth_relative_kernel( + x: ArrayLike, + y: ArrayLike, + alpha: float = 1.0, + kernel: _Kernel = "gaussian", + tail: float = 1e-12, + max_grid_points: int = 1_000_000, +) -> _FloatArray: + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + + if x.ndim != 1 or y.ndim != 1: + raise ValueError("x and y must be one-dimensional") + if x.size != y.size: + raise ValueError("x and y must have the same length") + if not np.isfinite(alpha) or alpha < 0: + raise ValueError("alpha must be non-negative") + if not np.isfinite(tail) or not (0.0 < tail < 1.0): + raise ValueError("tail must be between 0 and 1") + if isinstance(max_grid_points, bool | np.bool_) or not isinstance( + max_grid_points, int | np.integer + ): + raise TypeError("max_grid_points must be an integer") + if max_grid_points < 2: + raise ValueError("max_grid_points must be at least 2") + if np.any(~np.isfinite(x)): + raise ValueError("x must contain only finite values") + if np.any(x <= 0): + raise ValueError("x must be positive") + if np.any(np.diff(x) <= 0): + raise ValueError("x must be strictly increasing") + if x.size == 0 or alpha == 0 or x.size == 1: + return y.copy() + + logx = np.log(x) + dlog = np.diff(logx) + log_range = float(logx[-1] - logx[0]) + + # Preserve an existing geometric grid. Otherwise choose a geometric grid + # at least as dense as the smallest input spacing in log-space. + if np.allclose(dlog, dlog[0], rtol=1e-7, atol=0.0): + k = x.size + else: + k = int(np.ceil(log_range / np.min(dlog))) + 1 + + if k > max_grid_points: + raise ValueError( + "geometric resampling would require too many points, exceeding " + f"max_grid_points={max_grid_points:,}. Increase max_grid_points to " + "allow a larger grid." + ) + + xp = np.geomspace(x[0], x[-1], k) + yg = np.interp(xp, x, y) + zg = _smooth_relative_kernel_on_geomgrid( + yg, + alpha=alpha, + log_spacing=log_range / (k - 1), + kernel=kernel, + tail=tail, + ) + return cast(_FloatArray, np.interp(x, xp, zg)) + + +def _smooth_variable( + x: sc.Variable, + y: sc.Variable, + *, + alpha: float, + kernel: _Kernel, + tail: float, + max_grid_points: int, +) -> sc.Variable: + if x.ndim != 1 or y.ndim != 1: + raise sc.DimensionError("x and y must be one-dimensional") + if x.dims != y.dims: + raise sc.DimensionError("x and y must have the same dimension") + if x.is_binned or y.is_binned: + raise sc.DTypeError("x and y must not be binned") + if y.variances is not None: + raise sc.VariancesError( + "Smoothing signals with variances is not supported because it would " + "introduce correlations between data points." + ) + + return sc.array( + dims=y.dims, + values=_smooth_relative_kernel( + x.values, + y.values, + alpha=alpha, + kernel=kernel, + tail=tail, + max_grid_points=max_grid_points, + ), + unit=y.unit, + ) + + +def smooth_relative_kernel( + x: _ScippArray, + y: sc.Variable | None = None, + alpha: float = 1.0, + kernel: _Kernel = "gaussian", + tail: float = 1e-12, + max_grid_points: int = 1_000_000, +) -> _ScippArray: + """Smooth sampled data with a kernel of relative width. + + The kernel describes a distribution of relative displacements ``Z``, with + displaced coordinates given by ``x' = x * (1 + alpha * Z)``. At the + boundaries, the kernel is renormalized over the available finite input + domain. + + Input that is not geometrically spaced is interpolated to a geometric grid, + smoothed, and interpolated back to the original coordinates. + + Parameters + ---------- + x: + One-dimensional data to smooth, or positive, strictly increasing + one-dimensional sample coordinates. A data array must have a + dimension coordinate. + y: + Values to smooth when ``x`` contains the sample coordinates. Must be a + one-dimensional variable with the same dimension as ``x``. Must be + omitted when ``x`` is a data array. + alpha: + Scale factor for the relative-displacement distribution. Set to zero to + return a copy of the input without smoothing. + kernel: + Kernel distribution. Supported names are ``'gaussian'``, ``'rectangle'``, + and ``'triangle'``, including their aliases. Alternatively, provide a + fully specified distribution with ``cdf``, ``ppf``, and ``support`` + methods. + tail: + Total probability omitted when truncating a kernel with unbounded support + or support reaching the nonpositive coordinate domain. + max_grid_points: + Maximum permitted size of the intermediate geometric grid. Raises an + error rather than silently reducing resolution if this limit is exceeded. + + Returns + ------- + : + Smoothed data of the same type as the input. Coordinates and units are + preserved. + + Raises + ------ + ValueError + If the inputs have invalid values, if a data array has masks, if a + string does not identify a supported kernel, or if the required + intermediate grid exceeds ``max_grid_points``. + scipp.DimensionError + If the inputs are not one-dimensional or a pair of variables does not + have matching dimensions. + scipp.CoordError + If a data array has no dimension coordinate or has a bin-edge + coordinate. + scipp.VariancesError + If the signal has variances. + TypeError + If ``kernel`` is not a distribution-like object, or if + ``max_grid_points`` is not an integer. + """ + if isinstance(x, sc.DataArray): + if y is not None: + raise TypeError("y must be omitted when x is a DataArray") + if x.ndim != 1: + raise sc.DimensionError("data must be one-dimensional") + if x.dim not in x.coords: + raise sc.CoordError("data must have a dimension coordinate") + if x.coords.is_edges(x.dim): + raise sc.CoordError("the dimension coordinate must not contain bin edges") + if x.masks: + raise ValueError("smoothing data with masks is not supported") + + out = x.copy(deep=False) + out.data = _smooth_variable( + x.coords[x.dim], + x.data, + alpha=alpha, + kernel=kernel, + tail=tail, + max_grid_points=max_grid_points, + ) + return out + + if not isinstance(y, sc.Variable): + raise TypeError("expected a DataArray or a pair of Variables") + return _smooth_variable( + x, + y, + alpha=alpha, + kernel=kernel, + tail=tail, + max_grid_points=max_grid_points, + ) + + +def smooth_relative_gaussian( + x: _ScippArray, + y: sc.Variable | None = None, + alpha: float = 1.0, + tail: float = 1e-12, + max_grid_points: int = 1_000_000, +) -> _ScippArray: + """Smooth sampled data with a relative Gaussian kernel. + + ``alpha`` is the standard deviation of the Gaussian as a fraction of each + coordinate value. + + Parameters + ---------- + x: + One-dimensional data to smooth, or positive, strictly increasing + one-dimensional sample coordinates. A data array must have a + dimension coordinate. + y: + Values to smooth when ``x`` contains the sample coordinates. Must be a + one-dimensional variable with the same dimension as ``x``. Must be + omitted when ``x`` is a data array. + alpha: + Relative standard deviation of the Gaussian kernel. + tail: + Total Gaussian probability omitted when truncating the kernel. + max_grid_points: + Maximum permitted size of the intermediate geometric grid. + + Returns + ------- + : + Smoothed data of the same type as the input. + + See Also + -------- + smooth_relative_kernel: + Smooth with a named or user-provided relative kernel. + """ + return smooth_relative_kernel( + x, + y, + alpha=alpha, + kernel="gaussian", + tail=tail, + max_grid_points=max_grid_points, + ) + + +def smooth_relative_rectangle( + x: _ScippArray, + y: sc.Variable | None = None, + alpha: float = 1.0, + tail: float = 1e-12, + max_grid_points: int = 1_000_000, +) -> _ScippArray: + """Smooth sampled data with a relative rectangular kernel. + + The relative displacement is uniformly distributed on + ``[-alpha, alpha]``. + + Parameters + ---------- + x: + One-dimensional data to smooth, or positive, strictly increasing + one-dimensional sample coordinates. A data array must have a + dimension coordinate. + y: + Values to smooth when ``x`` contains the sample coordinates. Must be a + one-dimensional variable with the same dimension as ``x``. Must be + omitted when ``x`` is a data array. + alpha: + Half-width of the rectangular kernel relative to each coordinate value. + tail: + Total probability omitted if the kernel reaches the nonpositive + coordinate domain. + max_grid_points: + Maximum permitted size of the intermediate geometric grid. + + Returns + ------- + : + Smoothed data of the same type as the input. + + See Also + -------- + smooth_relative_kernel: + Smooth with a named or user-provided relative kernel. + """ + return smooth_relative_kernel( + x, + y, + alpha=alpha, + kernel="rectangle", + tail=tail, + max_grid_points=max_grid_points, + ) + + +def smooth_relative_triangle( + x: _ScippArray, + y: sc.Variable | None = None, + alpha: float = 1.0, + tail: float = 1e-12, + max_grid_points: int = 1_000_000, +) -> _ScippArray: + """Smooth sampled data with a relative triangular kernel. + + The relative displacement has symmetric triangular support on + ``[-alpha, alpha]`` and its peak at zero. + + Parameters + ---------- + x: + One-dimensional data to smooth, or positive, strictly increasing + one-dimensional sample coordinates. A data array must have a + dimension coordinate. + y: + Values to smooth when ``x`` contains the sample coordinates. Must be a + one-dimensional variable with the same dimension as ``x``. Must be + omitted when ``x`` is a data array. + alpha: + Half-width of the triangular kernel relative to each coordinate value. + tail: + Total probability omitted if the kernel reaches the nonpositive + coordinate domain. + max_grid_points: + Maximum permitted size of the intermediate geometric grid. + + Returns + ------- + : + Smoothed data of the same type as the input. + + See Also + -------- + smooth_relative_kernel: + Smooth with a named or user-provided relative kernel. + """ + return smooth_relative_kernel( + x, + y, + alpha=alpha, + kernel="triangle", + tail=tail, + max_grid_points=max_grid_points, + ) diff --git a/tests/smoothing_test.py b/tests/smoothing_test.py new file mode 100644 index 000000000..727f55807 --- /dev/null +++ b/tests/smoothing_test.py @@ -0,0 +1,454 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) + +import numpy as np +import pytest +import scipp as sc +from scipy.special import ndtr +from scipy.stats import uniform + +import scippneutron as scn +from scippneutron.smoothing import ( + _relative_kernel_weights, + _smooth_relative_kernel_on_geomgrid, + smooth_relative_gaussian, + smooth_relative_kernel, + smooth_relative_rectangle, + smooth_relative_triangle, +) + + +def _normal_pdf(z): + z = np.asarray(z, dtype=float) + return np.exp(-0.5 * z**2) / np.sqrt(2.0 * np.pi) + + +def _quadratic(x): + return 1.0 + 0.3 * x + 0.7 * x**2 + + +def _exact_smoothed_quadratic(x, alpha, lower, upper): + """Exactly smooth ``_quadratic`` with a finite-domain Gaussian.""" + x = np.asarray(x, dtype=float) + sigma = alpha * x + + z_lower = (lower - x) / sigma + z_upper = (upper - x) / sigma + normalization = ndtr(z_upper) - ndtr(z_lower) + + mean_z = (_normal_pdf(z_lower) - _normal_pdf(z_upper)) / normalization + mean_z_squared = ( + 1.0 + + (z_lower * _normal_pdf(z_lower) - z_upper * _normal_pdf(z_upper)) + / normalization + ) + + mean_x = x + sigma * mean_z + mean_x_squared = x**2 + 2.0 * x * sigma * mean_z + sigma**2 * mean_z_squared + return 1.0 + 0.3 * mean_x + 0.7 * mean_x_squared + + +def _geometric_cell_centers(lower, upper, size): + log_step = np.log(upper / lower) / size + return lower * np.exp((np.arange(size) + 0.5) * log_step) + + +def _variables(x, y, *, variances=None): + return ( + sc.array(dims=['x'], values=np.asarray(x), unit='m'), + sc.array( + dims=['x'], + values=np.asarray(y), + variances=None if variances is None else np.asarray(variances), + unit='counts', + ), + ) + + +def _smooth_values(smooth, x, y, **kwargs): + x, y = _variables(x, y) + return smooth(x, y, **kwargs).values + + +def _quadratic_smoothing_error(*, size, alpha, tail=1e-9): + lower = 0.1 + upper = 0.9 + x = _geometric_cell_centers(lower, upper, size) + actual = _smooth_values( + smooth_relative_gaussian, + x, + _quadratic(x), + alpha=alpha, + tail=tail, + ) + expected = _exact_smoothed_quadratic(x, alpha, lower, upper) + return actual - expected + + +def test_gaussian_smoothing_matches_exact_quadratic_on_geometric_grid(): + # The outer cell edges of this grid coincide with the finite integration + # domain used by _quadratic_smoothing_error. + error = _quadratic_smoothing_error(size=4000, alpha=0.1) + + assert np.max(np.abs(error)) < 5e-7 + assert np.sqrt(np.mean(error**2)) < 1e-7 + + +def test_gaussian_smoothing_error_is_second_order_in_log_grid_spacing(): + sizes = np.array([500, 1000, 2000, 4000]) + errors = np.array( + [ + np.max(np.abs(_quadratic_smoothing_error(size=size, alpha=0.1))) + for size in sizes + ] + ) + + # Halving the log-grid spacing should reduce midpoint quadrature error by + # four, corresponding to second-order convergence. + observed_orders = np.log2(errors[:-1] / errors[1:]) + np.testing.assert_allclose(observed_orders, 2.0, atol=0.06) + + +def test_gaussian_smoothing_error_scales_with_inverse_kernel_width(): + size = 4000 + alphas = np.array([0.025, 0.05, 0.1]) + errors = np.array( + [ + np.max(np.abs(_quadratic_smoothing_error(size=size, alpha=alpha))) + for alpha in alphas + ] + ) + + # For a well-resolved narrow Gaussian, the largest error is at a truncated + # boundary and scales as h**2 / alpha. + log_step = np.log(0.9 / 0.1) / size + scaled_errors = errors * alphas / log_step**2 + assert np.all(np.diff(errors) < 0.0) + np.testing.assert_allclose( + scaled_errors, + np.mean(scaled_errors), + rtol=0.08, + ) + + +def test_gaussian_smoothing_error_scales_with_tail_until_grid_error_dominates(): + tails = np.array([1e-2, 1e-3, 1e-4, 1e-5]) + truncation_errors = np.array( + [ + np.max(np.abs(_quadratic_smoothing_error(size=4000, alpha=0.1, tail=tail))) + for tail in tails + ] + ) + + # Gaussian moments in the omitted tails add logarithmic factors, so each + # decade should improve the result by approximately, but not exactly, ten. + reduction_per_decade = truncation_errors[:-1] / truncation_errors[1:] + assert np.all((5.0 < reduction_per_decade) & (reduction_per_decade < 15.0)) + + # Once tail truncation is negligible, reducing it further cannot improve + # the fixed-grid quadrature error. + grid_limited_errors = np.array( + [ + np.max(np.abs(_quadratic_smoothing_error(size=4000, alpha=0.1, tail=tail))) + for tail in (1e-9, 1e-12) + ] + ) + np.testing.assert_allclose( + grid_limited_errors[0], grid_limited_errors[1], rtol=0.05 + ) + + +@pytest.mark.parametrize( + ("smooth", "relative_variance", "max_error"), + [ + (smooth_relative_rectangle, 1.0 / 3.0, 3e-7), + (smooth_relative_triangle, 1.0 / 6.0, 5e-8), + ], +) +def test_compact_symmetric_kernel_matches_exact_interior_quadratic_moments( + smooth, relative_variance, max_error +): + lower = 0.1 + upper = 0.9 + size = 4000 + alpha = 0.1 + x = _geometric_cell_centers(lower, upper, size) + + actual = _smooth_values(smooth, x, _quadratic(x), alpha=alpha) + expected = 1.0 + 0.3 * x + 0.7 * x**2 * (1.0 + alpha**2 * relative_variance) + interior = (x * (1.0 - alpha) >= lower) & (x * (1.0 + alpha) <= upper) + + assert np.max(np.abs(actual[interior] - expected[interior])) < max_error + + +def test_smoothing_module_is_exposed_by_package(): + assert scn.smoothing.smooth_relative_gaussian is smooth_relative_gaussian + + +def test_accepts_data_array_and_preserves_metadata(): + x = sc.geomspace('x', 0.1, 0.9, 100, unit='m') + data = sc.DataArray( + sc.array(dims=['x'], values=_quadratic(x.values), unit='counts'), + coords={ + 'x': x, + 'aux': sc.arange('x', 100, unit='s'), + 'scalar': sc.scalar(1.2, unit='K'), + }, + ) + + actual = smooth_relative_gaussian(data, alpha=0.1) + expected = smooth_relative_gaussian(x, data.data, alpha=0.1) + + assert isinstance(actual, sc.DataArray) + assert sc.identical(actual.data, expected) + assert sc.identical(actual.coords['x'], data.coords['x']) + assert sc.identical(actual.coords['aux'], data.coords['aux']) + assert sc.identical(actual.coords['scalar'], data.coords['scalar']) + + +def test_rejects_variable_with_variances(): + x = sc.geomspace('x', 0.1, 0.9, 100, unit='m') + values = _quadratic(x.values) + y = sc.array( + dims=['x'], + values=values, + variances=2.0 + x.values, + unit='counts', + ) + + with pytest.raises(sc.VariancesError, match="signals with variances"): + smooth_relative_gaussian(x, y, alpha=0.1) + + +def test_rejects_data_array_with_variances(): + x = sc.geomspace('x', 0.1, 0.9, 100, unit='m') + data = sc.DataArray( + sc.array( + dims=['x'], + values=_quadratic(x.values), + variances=2.0 + x.values, + unit='counts', + ), + coords={'x': x}, + ) + + with pytest.raises(sc.VariancesError, match="signals with variances"): + smooth_relative_gaussian(data, alpha=0.1) + + +def test_rejects_numpy_arrays(): + with pytest.raises(TypeError, match="DataArray or a pair of Variables"): + smooth_relative_gaussian(np.arange(1.0, 4.0), np.ones(3)) + + +def test_rejects_variables_with_different_dimensions(): + x = sc.arange('x', 1.0, 4.0) + y = sc.ones(dims=['y'], shape=[3]) + + with pytest.raises(sc.DimensionError, match="same dimension"): + smooth_relative_gaussian(x, y) + + +def test_rejects_data_array_without_dimension_coordinate(): + data = sc.DataArray(sc.ones(dims=['x'], shape=[3])) + + with pytest.raises(sc.CoordError, match="dimension coordinate"): + smooth_relative_gaussian(data) + + +def test_rejects_data_array_with_bin_edge_coordinate(): + data = sc.DataArray( + sc.ones(dims=['x'], shape=[3]), + coords={'x': sc.arange('x', 1.0, 5.0)}, + ) + + with pytest.raises(sc.CoordError, match="bin edges"): + smooth_relative_gaussian(data) + + +def test_rejects_data_array_with_masks(): + data = sc.DataArray( + sc.ones(dims=['x'], shape=[3]), + coords={'x': sc.arange('x', 1.0, 4.0)}, + masks={'bad': sc.array(dims=['x'], values=[False, True, False])}, + ) + + with pytest.raises(ValueError, match="data with masks"): + smooth_relative_gaussian(data) + + +def test_rejects_y_with_data_array(): + x, y = _variables([1.0, 2.0], [3.0, 4.0]) + data = sc.DataArray(y, coords={'x': x}) + + with pytest.raises(TypeError, match="y must be omitted"): + smooth_relative_gaussian(data, y) + + +def test_rejects_geometric_grid_larger_than_limit(): + x = np.array([1.0, 1.01, 2.0]) + y = np.ones_like(x) + + with pytest.raises( + ValueError, + match=r"geometric resampling would require too many points.*max_grid_points=70", + ): + _smooth_values(smooth_relative_kernel, x, y, max_grid_points=70) + + +def test_accepts_geometric_grid_equal_to_limit(): + x = np.array([1.0, 1.01, 2.0]) + y = np.ones_like(x) + + actual = _smooth_values(smooth_relative_kernel, x, y, max_grid_points=71) + + np.testing.assert_allclose(actual, y) + + +def test_geometric_input_does_not_gain_a_point_from_roundoff(): + size = 100 + x = _geometric_cell_centers(0.1, 0.9, size) + y = _quadratic(x) + + actual = _smooth_values(smooth_relative_gaussian, x, y, max_grid_points=size) + + assert actual.shape == y.shape + + +@pytest.mark.parametrize( + "smooth", + [ + smooth_relative_gaussian, + smooth_relative_rectangle, + smooth_relative_triangle, + ], +) +def test_convenience_functions_forward_max_grid_points(smooth): + x = np.array([1.0, 1.01, 2.0]) + y = np.ones_like(x) + + with pytest.raises(ValueError, match="max_grid_points=2"): + _smooth_values(smooth, x, y, max_grid_points=2) + + +def test_pathologically_close_coordinates_fail_before_allocation(): + x = np.array([1.0, np.nextafter(1.0, 2.0), 2.0]) + y = np.ones_like(x) + + with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"): + _smooth_values(smooth_relative_kernel, x, y) + + +def test_wide_coordinate_range_does_not_overflow_grid_construction(): + x = np.array([1e-300, 1.0, 1e300]) + y = np.ones_like(x) + + actual = _smooth_values(smooth_relative_kernel, x, y) + + np.testing.assert_allclose(actual, y) + + +def test_kernel_stencil_is_bounded_before_allocation(): + max_offset = 1000 + + offsets, weights = _relative_kernel_weights( + log_spacing=1e-6, + alpha=1.0, + kernel="gaussian", + tail=1e-12, + max_offset=max_offset, + ) + + assert offsets[0] >= -max_offset + assert offsets[-1] <= max_offset + assert offsets.size <= 2 * max_offset + 1 + assert offsets.size == weights.size + + +def test_asymmetric_kernel_convolution_matches_direct_weighted_sum(): + size = 20 + log_spacing = 0.03 + alpha = 0.2 + kernel = uniform(loc=0.5, scale=1.0) + y = np.arange(size, dtype=float) ** 2 + + offsets, weights = _relative_kernel_weights( + log_spacing=log_spacing, + alpha=alpha, + kernel=kernel, + tail=1e-12, + max_offset=size - 1, + ) + expected = np.empty_like(y) + for i in range(size): + valid = (0 <= i + offsets) & (i + offsets < size) + denominator = np.sum(weights[valid]) + expected[i] = ( + np.dot(weights[valid], y[i + offsets[valid]]) / denominator + if denominator > 0.0 + else np.nan + ) + + actual = _smooth_relative_kernel_on_geomgrid( + y, + log_spacing=log_spacing, + alpha=alpha, + kernel=kernel, + tail=1e-12, + ) + + np.testing.assert_allclose(actual, expected, equal_nan=True) + + +def test_kernel_with_no_reachable_mass_returns_nan(): + actual = _smooth_relative_kernel_on_geomgrid( + np.arange(5.0), + log_spacing=0.1, + alpha=0.1, + kernel=uniform(loc=100.0, scale=1.0), + tail=1e-12, + ) + + assert np.all(np.isnan(actual)) + + +@pytest.mark.parametrize("alpha", [np.nan, np.inf, -np.inf, -1.0]) +def test_rejects_invalid_alpha_before_noop_return(alpha): + x, y = _variables([], []) + with pytest.raises(ValueError, match="alpha must be non-negative"): + smooth_relative_kernel(x, y, alpha=alpha) + + +@pytest.mark.parametrize("tail", [np.nan, np.inf, -np.inf, 0.0, 1.0]) +def test_rejects_invalid_tail_before_noop_return(tail): + x, y = _variables([], []) + with pytest.raises(ValueError, match="tail must be between 0 and 1"): + smooth_relative_kernel(x, y, tail=tail) + + +@pytest.mark.parametrize( + ("x", "message"), + [ + ([-1.0], "x must be positive"), + ([np.nan], "x must contain only finite values"), + ([np.inf], "x must contain only finite values"), + ], +) +def test_rejects_invalid_single_coordinate_before_noop_return(x, message): + x, y = _variables(x, [1.0]) + with pytest.raises(ValueError, match=message): + smooth_relative_kernel(x, y) + + +@pytest.mark.parametrize("max_grid_points", [True, 2.5]) +def test_rejects_non_integer_max_grid_points(max_grid_points): + x, y = _variables([], []) + with pytest.raises(TypeError, match="max_grid_points must be an integer"): + smooth_relative_kernel(x, y, max_grid_points=max_grid_points) + + +@pytest.mark.parametrize("max_grid_points", [-1, 0, 1]) +def test_rejects_too_small_max_grid_points(max_grid_points): + x, y = _variables([], []) + with pytest.raises(ValueError, match="max_grid_points must be at least 2"): + smooth_relative_kernel(x, y, max_grid_points=max_grid_points) From 7f204156793d3d6c9af1bc2f3c9a79ebe6283221 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 4 Aug 2026 16:02:34 +0200 Subject: [PATCH 2/6] feat: add similar functions for fixed-width smoothing kernels --- src/scippneutron/smoothing.py | 462 +++++++++++++++++++++++++++------- tests/smoothing_test.py | 237 ++++++++++++++++- 2 files changed, 598 insertions(+), 101 deletions(-) diff --git a/src/scippneutron/smoothing.py b/src/scippneutron/smoothing.py index cf1bc51e7..e6f7d757e 100644 --- a/src/scippneutron/smoothing.py +++ b/src/scippneutron/smoothing.py @@ -12,10 +12,14 @@ from scipy.stats import norm, triang, uniform __all__ = [ + "smooth_gaussian", + "smooth_kernel", + "smooth_rectangle", "smooth_relative_gaussian", "smooth_relative_kernel", "smooth_relative_rectangle", "smooth_relative_triangle", + "smooth_triangle", ] @@ -33,11 +37,10 @@ def support(self) -> tuple[Any, Any]: ... _ScippArray: TypeAlias = sc.DataArray | sc.Variable _BUILTIN_KERNELS: dict[str, _Distribution] = { - # Standard Gaussian in relative-coordinate units. + # Standard Gaussian. "gaussian": norm(), "normal": norm(), # Centered rectangle on [-1, 1]. - # With alpha=a, this has relative support [-a, a]. "rectangle": uniform(loc=-1.0, scale=2.0), "rect": uniform(loc=-1.0, scale=2.0), "box": uniform(loc=-1.0, scale=2.0), @@ -70,6 +73,33 @@ def _as_kernel_distribution(kernel: _Kernel) -> _Distribution: return kernel +def _kernel_support(dist: _Distribution) -> tuple[float, float]: + try: + support_min, support_max = dist.support() + except TypeError as e: + raise TypeError( + "kernel must be a fully specified distribution. For distributions " + "with shape parameters, pass a frozen distribution such as " + "triang(c=0.5), not triang." + ) from e + return float(support_min), float(support_max) + + +def _trim_kernel_weights( + offsets: _IntArray, weights: _FloatArray +) -> tuple[_IntArray, _FloatArray]: + nonzero = np.flatnonzero(weights > 0.0) + if nonzero.size == 0: + return offsets, weights + + # Trim zero-only ends, but preserve offset zero for convolution alignment. + zero = -offsets[0] + first = min(nonzero[0], zero) + last = max(nonzero[-1], zero) + 1 + weights = weights[first:last] + return offsets[first:last], weights / weights.sum() + + def _relative_kernel_weights( log_spacing: float, alpha: float, @@ -89,7 +119,6 @@ def _relative_kernel_weights( K(q, q') = 1 / (alpha * q) * f((q' - q) / (alpha * q)) where f is the PDF of the supplied distribution. - """ if not np.isfinite(alpha) or alpha <= 0: raise ValueError("alpha must be positive") @@ -115,17 +144,7 @@ def _relative_kernel_weights( if not np.isfinite(norm_mass) or norm_mass <= 0.0: raise ValueError("kernel has no positive-domain mass for this alpha") - try: - support_min, support_max = dist.support() - except TypeError as e: - raise TypeError( - "kernel must be a fully specified distribution. For distributions " - "with shape parameters, pass a frozen distribution such as " - "triang(c=0.5), not triang." - ) from e - - support_min = float(support_min) - support_max = float(support_max) + support_min, support_max = _kernel_support(dist) # Exact z-support after clipping to q' > 0. z_left_exact = max(z_domain_min, support_min) @@ -174,20 +193,44 @@ def _relative_kernel_weights( w = (dist.cdf(zU) - dist.cdf(zL)) / norm_mass w = np.maximum(w, 0.0) - nonzero = np.flatnonzero(w > 0.0) - if nonzero.size == 0: - # The distribution has no mass within reach of the finite input. - # Returning zero weights lets the caller mark those values as NaN. - return m, w + return _trim_kernel_weights(m, w) - # Trim zero-only ends, but preserve offset zero for convolution alignment. - zero = -m_min - first = min(nonzero[0], zero) - last = max(nonzero[-1], zero) + 1 - # Renormalize after finite tail truncation. - trimmed_weights = w[first:last] - return m[first:last], trimmed_weights / trimmed_weights.sum() +def _translation_invariant_kernel_weights( + spacing: float, + width: float, + kernel: _Kernel, + tail: float, + max_offset: int, +) -> tuple[_IntArray, _FloatArray]: + if not np.isfinite(width) or width <= 0: + raise ValueError("width must be positive") + if not np.isfinite(spacing) or spacing <= 0: + raise ValueError("spacing must be positive") + if not np.isfinite(tail) or not (0.0 < tail < 1.0): + raise ValueError("tail must be between 0 and 1") + if max_offset < 0: + raise ValueError("max_offset must be non-negative") + + dist = _as_kernel_distribution(kernel) + z_left, z_right = _kernel_support(dist) + if not np.all(np.isfinite([z_left, z_right])): + z_left, z_right = ( + float(value) for value in dist.ppf([0.5 * tail, 1.0 - 0.5 * tail]) + ) + + bounds = width * np.array([z_left, z_right]) + if not np.all(np.isfinite(bounds)): + raise ValueError("kernel bounds are not finite; increase tail") + + m_min = int(np.clip(np.floor(bounds[0] / spacing + 0.5), -max_offset, 0)) + m_max = int(np.clip(np.ceil(bounds[1] / spacing - 0.5), 0, max_offset)) + m = np.arange(m_min, m_max + 1, dtype=np.int64) + + lower = (m - 0.5) * spacing / width + upper = (m + 0.5) * spacing / width + weights = np.maximum(dist.cdf(upper) - dist.cdf(lower), 0.0) + return _trim_kernel_weights(m, weights) def _valid_weight_sums(n: int, m: _IntArray, w: _FloatArray) -> _FloatArray: @@ -218,6 +261,27 @@ def _valid_weight_sums(n: int, m: _IntArray, w: _FloatArray) -> _FloatArray: ) +def _smooth_with_weights( + y: _FloatArray, offsets: _IntArray, weights: _FloatArray +) -> _FloatArray: + # Desired operation: + # + # out[i] = sum_m weights[m] * y[i + m] + # + # scipy convolution reverses the second argument, hence weights[::-1]. + # FFT convolution would spread a single NaN or infinity over the entire + # output. SciPy recommends the direct method for non-finite inputs. + method = "auto" if np.all(np.isfinite(y)) else "direct" + full = cast(_FloatArray, convolve(y, weights[::-1], mode="full", method=method)) + start = int(offsets[-1]) + numerator = full[start : start + y.size] + denominator = _valid_weight_sums(y.size, offsets, weights) + + out = np.full_like(numerator, np.nan) + np.divide(numerator, denominator, out=out, where=denominator > 0.0) + return out + + def _smooth_relative_kernel_on_geomgrid( y: ArrayLike, log_spacing: float, @@ -238,34 +302,14 @@ def _smooth_relative_kernel_on_geomgrid( if y.size == 0 or alpha == 0: return y.copy() - m, w = _relative_kernel_weights( + offsets, weights = _relative_kernel_weights( log_spacing=log_spacing, alpha=alpha, kernel=kernel, tail=tail, max_offset=y.size - 1, ) - - # Desired operation: - # - # out[i] = sum_m w[m] * y[i + m] - # - # scipy convolution reverses the second argument, hence w[::-1]. - full = cast(_FloatArray, convolve(y, w[::-1], mode="full")) - - start = int(m[-1]) - numerator = full[start : start + y.size] - - denom = _valid_weight_sums(y.size, m, w) - - out = np.full_like(numerator, np.nan) - np.divide( - numerator, - denom, - out=out, - where=denom > 0.0, - ) - return out + return _smooth_with_weights(y, offsets, weights) def _smooth_relative_kernel( @@ -332,15 +376,66 @@ def _smooth_relative_kernel( return cast(_FloatArray, np.interp(x, xp, zg)) -def _smooth_variable( - x: sc.Variable, - y: sc.Variable, - *, - alpha: float, - kernel: _Kernel, - tail: float, - max_grid_points: int, -) -> sc.Variable: +def _smooth_kernel_values( + x: ArrayLike, + y: ArrayLike, + width: float, + kernel: _Kernel = "gaussian", + tail: float = 1e-12, +) -> _FloatArray: + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + + if x.ndim != 1 or y.ndim != 1: + raise ValueError("x and y must be one-dimensional") + if x.size != y.size: + raise ValueError("x and y must have the same length") + if not np.isfinite(width) or width < 0: + raise ValueError("width must be non-negative") + if not np.isfinite(tail) or not (0.0 < tail < 1.0): + raise ValueError("tail must be between 0 and 1") + if np.any(~np.isfinite(x)): + raise ValueError("x must contain only finite values") + if np.any(np.diff(x) <= 0): + raise ValueError("x must be strictly increasing") + if x.size == 0 or width == 0 or x.size == 1: + return y.copy() + + spacing = np.diff(x) + if not np.allclose(spacing, spacing[0], rtol=1e-7, atol=0.0): + raise ValueError("x must be regularly spaced") + + offsets, weights = _translation_invariant_kernel_weights( + spacing=float(spacing[0]), + width=width, + kernel=kernel, + tail=tail, + max_offset=y.size - 1, + ) + return _smooth_with_weights(y, offsets, weights) + + +def _scipp_input( + x: _ScippArray, y: sc.Variable | None +) -> tuple[sc.Variable, sc.Variable, sc.DataArray | None]: + template: sc.DataArray | None = None + if isinstance(x, sc.DataArray): + if y is not None: + raise TypeError("y must be omitted when x is a DataArray") + if x.ndim != 1: + raise sc.DimensionError("data must be one-dimensional") + if x.dim not in x.coords: + raise sc.CoordError("data must have a dimension coordinate") + if x.coords.is_edges(x.dim): + raise sc.CoordError("the dimension coordinate must not contain bin edges") + if x.masks: + raise ValueError("smoothing data with masks is not supported") + template = x + y = x.data + x = x.coords[x.dim] + elif not isinstance(y, sc.Variable): + raise TypeError("expected a DataArray or a pair of Variables") + if x.ndim != 1 or y.ndim != 1: raise sc.DimensionError("x and y must be one-dimensional") if x.dims != y.dims: @@ -352,19 +447,217 @@ def _smooth_variable( "Smoothing signals with variances is not supported because it would " "introduce correlations between data points." ) + return x, y, template + + +def _scipp_output( + template: sc.DataArray | None, y: sc.Variable, values: _FloatArray +) -> _ScippArray: + data = sc.array(dims=y.dims, values=values, unit=y.unit) + if template is None: + return data + + # The new container shares unchanged coordinates, but its data is independent. + out = template.copy(deep=False) + out.data = data + return out + + +def _width_in_coordinate_unit(width: sc.Variable, x: sc.Variable) -> float: + if not isinstance(width, sc.Variable): + raise TypeError("width must be a scipp.Variable") + if width.ndim != 0: + raise sc.DimensionError("width must be a scalar") + if width.variances is not None: + raise sc.VariancesError("kernel widths with variances are not supported") + return float(width.to(unit=x.unit).value) + + +def smooth_kernel( + x: _ScippArray, + y: sc.Variable | None = None, + *, + width: sc.Variable, + kernel: _Kernel = "gaussian", + tail: float = 1e-12, +) -> _ScippArray: + """Smooth regularly sampled data with a translation-invariant kernel. + + The kernel describes a distribution of displacements ``Z``, with displaced + coordinates given by ``x' = x + width * Z``. At the boundaries, the kernel + is renormalized over the available finite input domain. - return sc.array( - dims=y.dims, - values=_smooth_relative_kernel( - x.values, - y.values, - alpha=alpha, - kernel=kernel, - tail=tail, - max_grid_points=max_grid_points, - ), - unit=y.unit, + Parameters + ---------- + x: + One-dimensional data to smooth, or strictly increasing, regularly + spaced sample coordinates. A data array must have a dimension + coordinate. + y: + Values to smooth when ``x`` contains the sample coordinates. Must be a + one-dimensional variable with the same dimension as ``x``. Must be + omitted when ``x`` is a data array. + width: + Scale factor for the displacement distribution. Must be a scalar with a + unit compatible with the coordinate. Set to zero to return a copy of the + input without smoothing. + kernel: + Kernel distribution. Supported names are ``'gaussian'``, ``'rectangle'``, + and ``'triangle'``, including their aliases. Alternatively, provide a + fully specified distribution with ``cdf``, ``ppf``, and ``support`` + methods. + tail: + Total probability omitted when truncating a kernel with unbounded + support. + + Returns + ------- + : + Smoothed data of the same type as the input. Coordinates and units are + preserved. + + Raises + ------ + ValueError + If the coordinate is not regularly spaced or the inputs otherwise have + invalid values, if a data array has masks, or if a string does not + identify a supported kernel. + scipp.DimensionError + If the inputs are not one-dimensional, a pair of variables does not + have matching dimensions, or ``width`` is not scalar. + scipp.CoordError + If a data array has no dimension coordinate or has a bin-edge + coordinate. + scipp.UnitError + If the unit of ``width`` is incompatible with the coordinate unit. + scipp.VariancesError + If the signal or ``width`` has variances. + TypeError + If ``width`` is not a variable or ``kernel`` is not a distribution-like + object. + """ + x, y, template = _scipp_input(x, y) + values = _smooth_kernel_values( + x.values, + y.values, + width=_width_in_coordinate_unit(width, x), + kernel=kernel, + tail=tail, ) + return _scipp_output(template, y, values) + + +def smooth_gaussian( + x: _ScippArray, + y: sc.Variable | None = None, + *, + width: sc.Variable, + tail: float = 1e-12, +) -> _ScippArray: + """Smooth regularly sampled data with a fixed-width Gaussian kernel. + + Parameters + ---------- + x: + One-dimensional data to smooth, or strictly increasing, regularly + spaced sample coordinates. A data array must have a dimension + coordinate. + y: + Values to smooth when ``x`` contains the sample coordinates. Must be a + one-dimensional variable with the same dimension as ``x``. Must be + omitted when ``x`` is a data array. + width: + Standard deviation of the Gaussian. Must be a scalar with a unit + compatible with the coordinate. + tail: + Total Gaussian probability omitted when truncating the kernel. + + Returns + ------- + : + Smoothed data of the same type as the input. + + See Also + -------- + smooth_kernel: + Smooth with a named or user-provided translation-invariant kernel. + """ + return smooth_kernel(x, y, width=width, kernel="gaussian", tail=tail) + + +def smooth_rectangle( + x: _ScippArray, + y: sc.Variable | None = None, + *, + width: sc.Variable, +) -> _ScippArray: + """Smooth regularly sampled data with a fixed-width rectangular kernel. + + The displacement is uniformly distributed on ``[-width, width]``. + + Parameters + ---------- + x: + One-dimensional data to smooth, or strictly increasing, regularly + spaced sample coordinates. A data array must have a dimension + coordinate. + y: + Values to smooth when ``x`` contains the sample coordinates. Must be a + one-dimensional variable with the same dimension as ``x``. Must be + omitted when ``x`` is a data array. + width: + Half-width of the rectangular kernel. Must be a scalar with a unit + compatible with the coordinate. + + Returns + ------- + : + Smoothed data of the same type as the input. + + See Also + -------- + smooth_kernel: + Smooth with a named or user-provided translation-invariant kernel. + """ + return smooth_kernel(x, y, width=width, kernel="rectangle") + + +def smooth_triangle( + x: _ScippArray, + y: sc.Variable | None = None, + *, + width: sc.Variable, +) -> _ScippArray: + """Smooth regularly sampled data with a fixed-width triangular kernel. + + The displacement has symmetric triangular support on ``[-width, width]`` + and its peak at zero. + + Parameters + ---------- + x: + One-dimensional data to smooth, or strictly increasing, regularly + spaced sample coordinates. A data array must have a dimension + coordinate. + y: + Values to smooth when ``x`` contains the sample coordinates. Must be a + one-dimensional variable with the same dimension as ``x``. Must be + omitted when ``x`` is a data array. + width: + Half-width of the triangular kernel. Must be a scalar with a unit + compatible with the coordinate. + + Returns + ------- + : + Smoothed data of the same type as the input. + + See Also + -------- + smooth_kernel: + Smooth with a named or user-provided translation-invariant kernel. + """ + return smooth_kernel(x, y, width=width, kernel="triangle") def smooth_relative_kernel( @@ -434,39 +727,16 @@ def smooth_relative_kernel( If ``kernel`` is not a distribution-like object, or if ``max_grid_points`` is not an integer. """ - if isinstance(x, sc.DataArray): - if y is not None: - raise TypeError("y must be omitted when x is a DataArray") - if x.ndim != 1: - raise sc.DimensionError("data must be one-dimensional") - if x.dim not in x.coords: - raise sc.CoordError("data must have a dimension coordinate") - if x.coords.is_edges(x.dim): - raise sc.CoordError("the dimension coordinate must not contain bin edges") - if x.masks: - raise ValueError("smoothing data with masks is not supported") - - out = x.copy(deep=False) - out.data = _smooth_variable( - x.coords[x.dim], - x.data, - alpha=alpha, - kernel=kernel, - tail=tail, - max_grid_points=max_grid_points, - ) - return out - - if not isinstance(y, sc.Variable): - raise TypeError("expected a DataArray or a pair of Variables") - return _smooth_variable( - x, - y, + x, y, template = _scipp_input(x, y) + values = _smooth_relative_kernel( + x.values, + y.values, alpha=alpha, kernel=kernel, tail=tail, max_grid_points=max_grid_points, ) + return _scipp_output(template, y, values) def smooth_relative_gaussian( diff --git a/tests/smoothing_test.py b/tests/smoothing_test.py index 727f55807..79bccb148 100644 --- a/tests/smoothing_test.py +++ b/tests/smoothing_test.py @@ -11,10 +11,14 @@ from scippneutron.smoothing import ( _relative_kernel_weights, _smooth_relative_kernel_on_geomgrid, + smooth_gaussian, + smooth_kernel, + smooth_rectangle, smooth_relative_gaussian, smooth_relative_kernel, smooth_relative_rectangle, smooth_relative_triangle, + smooth_triangle, ) @@ -27,10 +31,9 @@ def _quadratic(x): return 1.0 + 0.3 * x + 0.7 * x**2 -def _exact_smoothed_quadratic(x, alpha, lower, upper): - """Exactly smooth ``_quadratic`` with a finite-domain Gaussian.""" +def _exact_smoothed_quadratic_with_sigma(x, sigma, lower, upper): x = np.asarray(x, dtype=float) - sigma = alpha * x + sigma = np.asarray(sigma, dtype=float) z_lower = (lower - x) / sigma z_upper = (upper - x) / sigma @@ -48,18 +51,27 @@ def _exact_smoothed_quadratic(x, alpha, lower, upper): return 1.0 + 0.3 * mean_x + 0.7 * mean_x_squared +def _exact_smoothed_quadratic(x, alpha, lower, upper): + """Exactly smooth ``_quadratic`` with a relative finite-domain Gaussian.""" + return _exact_smoothed_quadratic_with_sigma(x, alpha * np.asarray(x), lower, upper) + + def _geometric_cell_centers(lower, upper, size): log_step = np.log(upper / lower) / size return lower * np.exp((np.arange(size) + 0.5) * log_step) -def _variables(x, y, *, variances=None): +def _linear_cell_centers(lower, upper, size): + step = (upper - lower) / size + return lower + (np.arange(size) + 0.5) * step + + +def _variables(x, y): return ( sc.array(dims=['x'], values=np.asarray(x), unit='m'), sc.array( dims=['x'], values=np.asarray(y), - variances=None if variances is None else np.asarray(variances), unit='counts', ), ) @@ -85,6 +97,216 @@ def _quadratic_smoothing_error(*, size, alpha, tail=1e-9): return actual - expected +def _fixed_width_quadratic_smoothing_error(*, size, width, tail=1e-9): + lower = -0.4 + upper = 0.9 + x = _linear_cell_centers(lower, upper, size) + x_var, y_var = _variables(x, _quadratic(x)) + actual = smooth_gaussian( + x_var, + y_var, + width=sc.scalar(width, unit='m'), + tail=tail, + ).values + expected = _exact_smoothed_quadratic_with_sigma(x, width, lower, upper) + return actual - expected + + +def test_fixed_width_gaussian_matches_exact_quadratic_on_regular_grid(): + error = _fixed_width_quadratic_smoothing_error(size=4000, width=0.1) + + assert np.max(np.abs(error)) < 5e-7 + assert np.sqrt(np.mean(error**2)) < 1e-7 + + +def test_fixed_width_gaussian_error_is_second_order_in_grid_spacing(): + sizes = np.array([500, 1000, 2000, 4000]) + errors = np.array( + [ + np.max(np.abs(_fixed_width_quadratic_smoothing_error(size=size, width=0.1))) + for size in sizes + ] + ) + + observed_orders = np.log2(errors[:-1] / errors[1:]) + np.testing.assert_allclose(observed_orders, 2.0, atol=0.06) + + +def test_fixed_width_gaussian_error_scales_with_tail_until_grid_error_dominates(): + tails = np.array([1e-2, 1e-3, 1e-4, 1e-5]) + errors = np.array( + [ + np.max( + np.abs( + _fixed_width_quadratic_smoothing_error( + size=4000, width=0.1, tail=tail + ) + ) + ) + for tail in tails + ] + ) + + reduction_per_decade = errors[:-1] / errors[1:] + assert np.all((5.0 < reduction_per_decade) & (reduction_per_decade < 15.0)) + + grid_limited_errors = np.array( + [ + np.max( + np.abs( + _fixed_width_quadratic_smoothing_error( + size=4000, width=0.1, tail=tail + ) + ) + ) + for tail in (1e-9, 1e-12) + ] + ) + np.testing.assert_allclose( + grid_limited_errors[0], grid_limited_errors[1], rtol=0.05 + ) + + +@pytest.mark.parametrize( + ("smooth", "kernel_variance", "max_error"), + [ + (smooth_rectangle, 1.0 / 3.0, 3e-7), + (smooth_triangle, 1.0 / 6.0, 5e-8), + ], +) +def test_fixed_width_compact_kernel_matches_interior_quadratic_moments( + smooth, kernel_variance, max_error +): + lower = -0.4 + upper = 0.9 + size = 4000 + width = 0.1 + x = _linear_cell_centers(lower, upper, size) + x_var, y_var = _variables(x, _quadratic(x)) + + actual = smooth(x_var, y_var, width=sc.scalar(width, unit='m')).values + expected = 1.0 + 0.3 * x + 0.7 * (x**2 + width**2 * kernel_variance) + interior = (x - width >= lower) & (x + width <= upper) + + assert np.max(np.abs(actual[interior] - expected[interior])) < max_error + + +def test_fixed_width_smoothing_converts_width_to_coordinate_unit(): + x = sc.linspace('x', -0.4, 0.9, 100, unit='m') + y = sc.array(dims=['x'], values=_quadratic(x.values), unit='counts') + + in_meters = smooth_gaussian(x, y, width=sc.scalar(0.1, unit='m')) + in_centimeters = smooth_gaussian(x, y, width=sc.scalar(10.0, unit='cm')) + + assert sc.identical(in_meters, in_centimeters) + + +def test_fixed_width_smoothing_accepts_distribution(): + x = sc.linspace('x', -0.4, 0.9, 100, unit='m') + y = sc.array(dims=['x'], values=_quadratic(x.values), unit='counts') + width = sc.scalar(0.1, unit='m') + + actual = smooth_kernel( + x, + y, + width=width, + kernel=uniform(loc=-1.0, scale=2.0), + ) + expected = smooth_rectangle(x, y, width=width) + + assert sc.identical(actual, expected) + + +def test_fixed_width_asymmetric_kernel_has_correct_direction(): + x = sc.arange('x', 0.0, 2.0, 0.1, unit='m') + y = sc.array(dims=['x'], values=x.values, unit='counts') + + actual = smooth_kernel( + x, + y, + width=sc.scalar(0.2, unit='m'), + kernel=uniform(loc=0.5, scale=1.0), + ) + + np.testing.assert_allclose(actual.values[3:-3], y.values[3:-3] + 0.2) + + +def test_nonfinite_value_only_affects_overlapping_kernel_windows(): + size = 10_000 + x = sc.arange('x', float(size)) + values = np.ones(size) + values[size // 2] = np.nan + y = sc.array(dims=['x'], values=values) + + actual = smooth_gaussian(x, y, width=sc.scalar(150.0)) + + assert np.isfinite(actual.values[0]) + assert np.isnan(actual.values[size // 2]) + assert np.isfinite(actual.values[-1]) + + +def test_fixed_width_smoothing_accepts_data_array(): + x = sc.linspace('x', -0.4, 0.9, 100, unit='m') + data = sc.DataArray( + sc.array(dims=['x'], values=_quadratic(x.values), unit='counts'), + coords={'x': x, 'aux': sc.arange('x', 100)}, + ) + + actual = smooth_gaussian(data, width=sc.scalar(0.1, unit='m')) + expected = smooth_gaussian(x, data.data, width=sc.scalar(0.1, unit='m')) + + assert sc.identical(actual.data, expected) + assert sc.identical(actual.coords['x'], data.coords['x']) + assert sc.identical(actual.coords['aux'], data.coords['aux']) + + +def test_fixed_width_smoothing_rejects_nonuniform_grid(): + x, y = _variables([0.0, 1.0, 2.1], [1.0, 2.0, 3.0]) + + with pytest.raises(ValueError, match="regularly spaced"): + smooth_kernel(x, y, width=sc.scalar(0.1, unit='m')) + + +@pytest.mark.parametrize("width", [np.nan, np.inf, -np.inf, -1.0]) +def test_fixed_width_smoothing_rejects_invalid_width(width): + x, y = _variables([], []) + + with pytest.raises(ValueError, match="width must be non-negative"): + smooth_kernel(x, y, width=sc.scalar(width, unit='m')) + + +def test_fixed_width_smoothing_rejects_non_scalar_width(): + x, y = _variables([0.0, 1.0], [1.0, 2.0]) + + with pytest.raises(sc.DimensionError, match="width must be a scalar"): + smooth_kernel(x, y, width=sc.array(dims=['width'], values=[0.1], unit='m')) + + +def test_fixed_width_smoothing_rejects_non_variable_width(): + x, y = _variables([0.0, 1.0], [1.0, 2.0]) + + with pytest.raises(TypeError, match=r"width must be a scipp\.Variable"): + smooth_kernel(x, y, width=0.1) # type: ignore[arg-type] + + +def test_fixed_width_smoothing_rejects_incompatible_width_unit(): + x, y = _variables([0.0, 1.0], [1.0, 2.0]) + + with pytest.raises(sc.UnitError): + smooth_kernel(x, y, width=sc.scalar(0.1, unit='s')) + + +def test_fixed_width_smoothing_rejects_width_with_variance(): + x, y = _variables([0.0, 1.0], [1.0, 2.0]) + + with pytest.raises(sc.VariancesError, match="widths with variances"): + smooth_kernel( + x, + y, + width=sc.scalar(0.1, variance=0.01, unit='m'), + ) + + def test_gaussian_smoothing_matches_exact_quadratic_on_geometric_grid(): # The outer cell edges of this grid coincide with the finite integration # domain used by _quadratic_smoothing_error. @@ -195,16 +417,21 @@ def test_accepts_data_array_and_preserves_metadata(): 'scalar': sc.scalar(1.2, unit='K'), }, ) + original = data.copy() actual = smooth_relative_gaussian(data, alpha=0.1) expected = smooth_relative_gaussian(x, data.data, alpha=0.1) + assert sc.identical(data, original) assert isinstance(actual, sc.DataArray) assert sc.identical(actual.data, expected) assert sc.identical(actual.coords['x'], data.coords['x']) assert sc.identical(actual.coords['aux'], data.coords['aux']) assert sc.identical(actual.coords['scalar'], data.coords['scalar']) + actual.values[0] = -1.0 + assert sc.identical(data, original) + def test_rejects_variable_with_variances(): x = sc.geomspace('x', 0.1, 0.9, 100, unit='m') From ce1ea58f154f8f68fcdf21c4282ff329400f6a8c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 4 Aug 2026 17:02:17 +0200 Subject: [PATCH 3/6] fix: support non-uniform grids in fixed-width kernel smoothing procedures --- src/scippneutron/smoothing.py | 125 ++++++++++++++++++++++++---------- tests/smoothing_test.py | 106 ++++++++++++++++++++++++++-- 2 files changed, 189 insertions(+), 42 deletions(-) diff --git a/src/scippneutron/smoothing.py b/src/scippneutron/smoothing.py index e6f7d757e..9d52a5955 100644 --- a/src/scippneutron/smoothing.py +++ b/src/scippneutron/smoothing.py @@ -312,6 +312,15 @@ def _smooth_relative_kernel_on_geomgrid( return _smooth_with_weights(y, offsets, weights) +def _validate_max_grid_points(max_grid_points: int) -> None: + if isinstance(max_grid_points, bool | np.bool_) or not isinstance( + max_grid_points, int | np.integer + ): + raise TypeError("max_grid_points must be an integer") + if max_grid_points < 2: + raise ValueError("max_grid_points must be at least 2") + + def _smooth_relative_kernel( x: ArrayLike, y: ArrayLike, @@ -331,12 +340,7 @@ def _smooth_relative_kernel( raise ValueError("alpha must be non-negative") if not np.isfinite(tail) or not (0.0 < tail < 1.0): raise ValueError("tail must be between 0 and 1") - if isinstance(max_grid_points, bool | np.bool_) or not isinstance( - max_grid_points, int | np.integer - ): - raise TypeError("max_grid_points must be an integer") - if max_grid_points < 2: - raise ValueError("max_grid_points must be at least 2") + _validate_max_grid_points(max_grid_points) if np.any(~np.isfinite(x)): raise ValueError("x must contain only finite values") if np.any(x <= 0): @@ -382,6 +386,7 @@ def _smooth_kernel_values( width: float, kernel: _Kernel = "gaussian", tail: float = 1e-12, + max_grid_points: int = 1_000_000, ) -> _FloatArray: x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) @@ -394,6 +399,7 @@ def _smooth_kernel_values( raise ValueError("width must be non-negative") if not np.isfinite(tail) or not (0.0 < tail < 1.0): raise ValueError("tail must be between 0 and 1") + _validate_max_grid_points(max_grid_points) if np.any(~np.isfinite(x)): raise ValueError("x must contain only finite values") if np.any(np.diff(x) <= 0): @@ -402,17 +408,34 @@ def _smooth_kernel_values( return y.copy() spacing = np.diff(x) - if not np.allclose(spacing, spacing[0], rtol=1e-7, atol=0.0): - raise ValueError("x must be regularly spaced") + x_range = float(x[-1] - x[0]) + + # Preserve an existing uniform grid. Otherwise choose a uniform grid at + # least as dense as the smallest input spacing. + if np.allclose(spacing, spacing[0], rtol=1e-7, atol=0.0): + k = x.size + else: + k = int(np.ceil(x_range / np.min(spacing))) + 1 + + if k > max_grid_points: + raise ValueError( + "uniform resampling would require too many points, exceeding " + f"max_grid_points={max_grid_points:,}. Increase max_grid_points to " + "allow a larger grid." + ) offsets, weights = _translation_invariant_kernel_weights( - spacing=float(spacing[0]), + spacing=x_range / (k - 1), width=width, kernel=kernel, tail=tail, - max_offset=y.size - 1, + max_offset=k - 1, ) - return _smooth_with_weights(y, offsets, weights) + xp = np.linspace(x[0], x[-1], k) + yg = np.interp(xp, x, y) + + zg = _smooth_with_weights(yg, offsets, weights) + return cast(_FloatArray, np.interp(x, xp, zg)) def _scipp_input( @@ -480,19 +503,22 @@ def smooth_kernel( width: sc.Variable, kernel: _Kernel = "gaussian", tail: float = 1e-12, + max_grid_points: int = 1_000_000, ) -> _ScippArray: - """Smooth regularly sampled data with a translation-invariant kernel. + """Smooth sampled data with a translation-invariant kernel. The kernel describes a distribution of displacements ``Z``, with displaced coordinates given by ``x' = x + width * Z``. At the boundaries, the kernel is renormalized over the available finite input domain. + Input that is not uniformly spaced is interpolated to a uniform grid, + smoothed, and interpolated back to the original coordinates. + Parameters ---------- x: - One-dimensional data to smooth, or strictly increasing, regularly - spaced sample coordinates. A data array must have a dimension - coordinate. + One-dimensional data to smooth, or strictly increasing sample + coordinates. A data array must have a dimension coordinate. y: Values to smooth when ``x`` contains the sample coordinates. Must be a one-dimensional variable with the same dimension as ``x``. Must be @@ -509,6 +535,9 @@ def smooth_kernel( tail: Total probability omitted when truncating a kernel with unbounded support. + max_grid_points: + Maximum permitted size of the intermediate uniform grid. Raises an + error rather than silently reducing resolution if this limit is exceeded. Returns ------- @@ -519,9 +548,9 @@ def smooth_kernel( Raises ------ ValueError - If the coordinate is not regularly spaced or the inputs otherwise have - invalid values, if a data array has masks, or if a string does not - identify a supported kernel. + If the inputs have invalid values, if a data array has masks, if a string + does not identify a supported kernel, or if the required intermediate + grid exceeds ``max_grid_points``. scipp.DimensionError If the inputs are not one-dimensional, a pair of variables does not have matching dimensions, or ``width`` is not scalar. @@ -533,8 +562,8 @@ def smooth_kernel( scipp.VariancesError If the signal or ``width`` has variances. TypeError - If ``width`` is not a variable or ``kernel`` is not a distribution-like - object. + If ``width`` is not a variable, ``kernel`` is not a distribution-like + object, or ``max_grid_points`` is not an integer. """ x, y, template = _scipp_input(x, y) values = _smooth_kernel_values( @@ -543,6 +572,7 @@ def smooth_kernel( width=_width_in_coordinate_unit(width, x), kernel=kernel, tail=tail, + max_grid_points=max_grid_points, ) return _scipp_output(template, y, values) @@ -553,15 +583,15 @@ def smooth_gaussian( *, width: sc.Variable, tail: float = 1e-12, + max_grid_points: int = 1_000_000, ) -> _ScippArray: - """Smooth regularly sampled data with a fixed-width Gaussian kernel. + """Smooth sampled data with a fixed-width Gaussian kernel. Parameters ---------- x: - One-dimensional data to smooth, or strictly increasing, regularly - spaced sample coordinates. A data array must have a dimension - coordinate. + One-dimensional data to smooth, or strictly increasing sample + coordinates. A data array must have a dimension coordinate. y: Values to smooth when ``x`` contains the sample coordinates. Must be a one-dimensional variable with the same dimension as ``x``. Must be @@ -571,6 +601,8 @@ def smooth_gaussian( compatible with the coordinate. tail: Total Gaussian probability omitted when truncating the kernel. + max_grid_points: + Maximum permitted size of the intermediate uniform grid. Returns ------- @@ -582,7 +614,14 @@ def smooth_gaussian( smooth_kernel: Smooth with a named or user-provided translation-invariant kernel. """ - return smooth_kernel(x, y, width=width, kernel="gaussian", tail=tail) + return smooth_kernel( + x, + y, + width=width, + kernel="gaussian", + tail=tail, + max_grid_points=max_grid_points, + ) def smooth_rectangle( @@ -590,17 +629,17 @@ def smooth_rectangle( y: sc.Variable | None = None, *, width: sc.Variable, + max_grid_points: int = 1_000_000, ) -> _ScippArray: - """Smooth regularly sampled data with a fixed-width rectangular kernel. + """Smooth sampled data with a fixed-width rectangular kernel. The displacement is uniformly distributed on ``[-width, width]``. Parameters ---------- x: - One-dimensional data to smooth, or strictly increasing, regularly - spaced sample coordinates. A data array must have a dimension - coordinate. + One-dimensional data to smooth, or strictly increasing sample + coordinates. A data array must have a dimension coordinate. y: Values to smooth when ``x`` contains the sample coordinates. Must be a one-dimensional variable with the same dimension as ``x``. Must be @@ -608,6 +647,8 @@ def smooth_rectangle( width: Half-width of the rectangular kernel. Must be a scalar with a unit compatible with the coordinate. + max_grid_points: + Maximum permitted size of the intermediate uniform grid. Returns ------- @@ -619,7 +660,13 @@ def smooth_rectangle( smooth_kernel: Smooth with a named or user-provided translation-invariant kernel. """ - return smooth_kernel(x, y, width=width, kernel="rectangle") + return smooth_kernel( + x, + y, + width=width, + kernel="rectangle", + max_grid_points=max_grid_points, + ) def smooth_triangle( @@ -627,8 +674,9 @@ def smooth_triangle( y: sc.Variable | None = None, *, width: sc.Variable, + max_grid_points: int = 1_000_000, ) -> _ScippArray: - """Smooth regularly sampled data with a fixed-width triangular kernel. + """Smooth sampled data with a fixed-width triangular kernel. The displacement has symmetric triangular support on ``[-width, width]`` and its peak at zero. @@ -636,9 +684,8 @@ def smooth_triangle( Parameters ---------- x: - One-dimensional data to smooth, or strictly increasing, regularly - spaced sample coordinates. A data array must have a dimension - coordinate. + One-dimensional data to smooth, or strictly increasing sample + coordinates. A data array must have a dimension coordinate. y: Values to smooth when ``x`` contains the sample coordinates. Must be a one-dimensional variable with the same dimension as ``x``. Must be @@ -646,6 +693,8 @@ def smooth_triangle( width: Half-width of the triangular kernel. Must be a scalar with a unit compatible with the coordinate. + max_grid_points: + Maximum permitted size of the intermediate uniform grid. Returns ------- @@ -657,7 +706,13 @@ def smooth_triangle( smooth_kernel: Smooth with a named or user-provided translation-invariant kernel. """ - return smooth_kernel(x, y, width=width, kernel="triangle") + return smooth_kernel( + x, + y, + width=width, + kernel="triangle", + max_grid_points=max_grid_points, + ) def smooth_relative_kernel( diff --git a/tests/smoothing_test.py b/tests/smoothing_test.py index 79bccb148..e7d454a40 100644 --- a/tests/smoothing_test.py +++ b/tests/smoothing_test.py @@ -260,10 +260,88 @@ def test_fixed_width_smoothing_accepts_data_array(): assert sc.identical(actual.coords['aux'], data.coords['aux']) -def test_fixed_width_smoothing_rejects_nonuniform_grid(): - x, y = _variables([0.0, 1.0, 2.1], [1.0, 2.0, 3.0]) +def test_fixed_width_smoothing_resamples_nonuniform_grid(): + x_values = np.array([0.0, 0.25, 0.75, 1.0]) + y_values = _quadratic(x_values) + width = sc.scalar(0.2, unit='m') + x, y = _variables(x_values, y_values) - with pytest.raises(ValueError, match="regularly spaced"): + actual = smooth_gaussian(x, y, width=width) + + xp_values = np.linspace(0.0, 1.0, 5) + xp, yp = _variables(xp_values, np.interp(xp_values, x_values, y_values)) + smoothed_yp = smooth_gaussian(xp, yp, width=width) + expected = np.interp(x_values, xp_values, smoothed_yp.values) + np.testing.assert_allclose(actual.values, expected) + + +def test_fixed_width_gaussian_matches_exact_quadratic_on_nonuniform_grid(): + lower = -0.4 + upper = 0.9 + size = 2000 + width = 0.1 + spacing = (upper - lower) / size + x = _linear_cell_centers(lower, upper, size) + x += 0.2 * spacing * np.sin(np.linspace(0.0, 8.0 * np.pi, size)) + x_var, y_var = _variables(x, _quadratic(x)) + + actual = smooth_gaussian( + x_var, + y_var, + width=sc.scalar(width, unit='m'), + tail=1e-9, + ).values + + intervals = int(np.ceil((x[-1] - x[0]) / np.min(np.diff(x)))) + resampled_spacing = (x[-1] - x[0]) / intervals + expected = _exact_smoothed_quadratic_with_sigma( + x, + width, + x[0] - 0.5 * resampled_spacing, + x[-1] + 0.5 * resampled_spacing, + ) + + assert np.max(np.abs(actual - expected)) < 5e-7 + + +def test_fixed_width_smoothing_rejects_resampled_grid_larger_than_limit(): + x, y = _variables([0.0, 0.25, 1.0], [1.0, 2.0, 3.0]) + + with pytest.raises( + ValueError, + match=r"uniform resampling would require too many points.*max_grid_points=4", + ): + smooth_kernel(x, y, width=sc.scalar(0.1, unit='m'), max_grid_points=4) + + +def test_fixed_width_smoothing_accepts_resampled_grid_equal_to_limit(): + x, y = _variables([0.0, 0.25, 1.0], [1.0, 1.0, 1.0]) + + actual = smooth_kernel( + x, + y, + width=sc.scalar(0.1, unit='m'), + max_grid_points=5, + ) + + np.testing.assert_allclose(actual.values, y.values) + + +@pytest.mark.parametrize("smooth", [smooth_gaussian, smooth_rectangle, smooth_triangle]) +def test_fixed_width_convenience_functions_forward_max_grid_points(smooth): + x, y = _variables([0.0, 0.25, 1.0], [1.0, 2.0, 3.0]) + + with pytest.raises(ValueError, match="max_grid_points=2"): + smooth(x, y, width=sc.scalar(0.1, unit='m'), max_grid_points=2) + + +def test_fixed_width_pathologically_close_coordinates_fail_before_allocation(): + x, y = _variables( + [1.0, np.nextafter(1.0, 2.0), 2.0], + [1.0, 1.0, 1.0], + ) + + with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"): smooth_kernel(x, y, width=sc.scalar(0.1, unit='m')) @@ -667,15 +745,29 @@ def test_rejects_invalid_single_coordinate_before_noop_return(x, message): smooth_relative_kernel(x, y) +@pytest.mark.parametrize( + ("smooth", "kwargs"), + [ + (smooth_relative_kernel, {}), + (smooth_kernel, {"width": sc.scalar(0.1, unit='m')}), + ], +) @pytest.mark.parametrize("max_grid_points", [True, 2.5]) -def test_rejects_non_integer_max_grid_points(max_grid_points): +def test_rejects_non_integer_max_grid_points(smooth, kwargs, max_grid_points): x, y = _variables([], []) with pytest.raises(TypeError, match="max_grid_points must be an integer"): - smooth_relative_kernel(x, y, max_grid_points=max_grid_points) + smooth(x, y, max_grid_points=max_grid_points, **kwargs) +@pytest.mark.parametrize( + ("smooth", "kwargs"), + [ + (smooth_relative_kernel, {}), + (smooth_kernel, {"width": sc.scalar(0.1, unit='m')}), + ], +) @pytest.mark.parametrize("max_grid_points", [-1, 0, 1]) -def test_rejects_too_small_max_grid_points(max_grid_points): +def test_rejects_too_small_max_grid_points(smooth, kwargs, max_grid_points): x, y = _variables([], []) with pytest.raises(ValueError, match="max_grid_points must be at least 2"): - smooth_relative_kernel(x, y, max_grid_points=max_grid_points) + smooth(x, y, max_grid_points=max_grid_points, **kwargs) From bd515a03827a322e300424933a51f69af3f347b7 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 4 Aug 2026 17:29:10 +0200 Subject: [PATCH 4/6] refactor: only expose two functions, smooth and smooth_relative, remove the kernel specific helpers --- src/scippneutron/smoothing.py | 465 +++++++--------------------------- tests/smoothing_test.py | 255 ++++++++++--------- 2 files changed, 233 insertions(+), 487 deletions(-) diff --git a/src/scippneutron/smoothing.py b/src/scippneutron/smoothing.py index 9d52a5955..adb8cec8d 100644 --- a/src/scippneutron/smoothing.py +++ b/src/scippneutron/smoothing.py @@ -3,6 +3,7 @@ from __future__ import annotations +from numbers import Real from typing import Any, Protocol, TypeAlias, cast import numpy as np @@ -12,14 +13,8 @@ from scipy.stats import norm, triang, uniform __all__ = [ - "smooth_gaussian", - "smooth_kernel", - "smooth_rectangle", - "smooth_relative_gaussian", - "smooth_relative_kernel", - "smooth_relative_rectangle", - "smooth_relative_triangle", - "smooth_triangle", + "smooth", + "smooth_relative", ] @@ -44,6 +39,7 @@ def support(self) -> tuple[Any, Any]: ... "rectangle": uniform(loc=-1.0, scale=2.0), "rect": uniform(loc=-1.0, scale=2.0), "box": uniform(loc=-1.0, scale=2.0), + "boxcar": uniform(loc=-1.0, scale=2.0), "uniform": uniform(loc=-1.0, scale=2.0), # Centered triangle on [-1, 1], peak at 0. # Users can pass triang(c=...) themselves for asymmetric triangles. @@ -102,7 +98,7 @@ def _trim_kernel_weights( def _relative_kernel_weights( log_spacing: float, - alpha: float, + scale: float, kernel: _Kernel, tail: float, max_offset: int, @@ -112,16 +108,16 @@ def _relative_kernel_weights( The kernel distribution describes the relative displacement Z: - q' = q * (1 + alpha * Z) + q' = q * (1 + scale * Z) Equivalently, - K(q, q') = 1 / (alpha * q) * f((q' - q) / (alpha * q)) + K(q, q') = 1 / (scale * q) * f((q' - q) / (scale * q)) where f is the PDF of the supplied distribution. """ - if not np.isfinite(alpha) or alpha <= 0: - raise ValueError("alpha must be positive") + if not np.isfinite(scale) or scale <= 0: + raise ValueError("scale must be positive") if not np.isfinite(log_spacing) or log_spacing <= 0: raise ValueError("log_spacing must be positive") if not np.isfinite(tail) or not (0.0 < tail < 1.0): @@ -134,15 +130,15 @@ def _relative_kernel_weights( # Positive physical domain: # # q' > 0 - # q * (1 + alpha * z) > 0 - # z > -1 / alpha - z_domain_min = -1.0 / alpha + # q * (1 + scale * z) > 0 + # z > -1 / scale + z_domain_min = -1.0 / scale p_domain_min = float(dist.cdf(z_domain_min)) norm_mass = 1.0 - p_domain_min if not np.isfinite(norm_mass) or norm_mass <= 0.0: - raise ValueError("kernel has no positive-domain mass for this alpha") + raise ValueError("kernel has no positive-domain mass for this scale") support_min, support_max = _kernel_support(dist) @@ -155,20 +151,20 @@ def _relative_kernel_weights( has_finite_log_support = ( np.isfinite(z_left_exact) and np.isfinite(z_right_exact) - and (1.0 + alpha * z_left_exact > 0.0) + and (1.0 + scale * z_left_exact > 0.0) ) if has_finite_log_support: - u_left = np.log1p(alpha * z_left_exact) - u_right = np.log1p(alpha * z_right_exact) + u_left = np.log1p(scale * z_left_exact) + u_right = np.log1p(scale * z_right_exact) else: probabilities = ( p_domain_min + np.array([0.5 * tail, 1.0 - 0.5 * tail]) * norm_mass ) z_left, z_right = (float(value) for value in dist.ppf(probabilities)) - u_left = np.log1p(alpha * z_left) - u_right = np.log1p(alpha * z_right) + u_left = np.log1p(scale * z_left) + u_right = np.log1p(scale * z_right) if not np.isfinite(u_right): raise ValueError("right kernel bound is not finite; increase tail") @@ -184,10 +180,10 @@ def _relative_kernel_weights( L = (m - 0.5) * log_spacing U = (m + 0.5) * log_spacing - # z = (q' - q) / (alpha q) - # = (exp(u) - 1) / alpha - zL = np.expm1(L) / alpha - zU = np.expm1(U) / alpha + # z = (q' - q) / (scale q) + # = (exp(u) - 1) / scale + zL = np.expm1(L) / scale + zU = np.expm1(U) / scale # Exact cell-integrated weights in log-space. w = (dist.cdf(zU) - dist.cdf(zL)) / norm_mass @@ -198,13 +194,13 @@ def _relative_kernel_weights( def _translation_invariant_kernel_weights( spacing: float, - width: float, + scale: float, kernel: _Kernel, tail: float, max_offset: int, ) -> tuple[_IntArray, _FloatArray]: - if not np.isfinite(width) or width <= 0: - raise ValueError("width must be positive") + if not np.isfinite(scale) or scale <= 0: + raise ValueError("scale must be positive") if not np.isfinite(spacing) or spacing <= 0: raise ValueError("spacing must be positive") if not np.isfinite(tail) or not (0.0 < tail < 1.0): @@ -219,7 +215,7 @@ def _translation_invariant_kernel_weights( float(value) for value in dist.ppf([0.5 * tail, 1.0 - 0.5 * tail]) ) - bounds = width * np.array([z_left, z_right]) + bounds = scale * np.array([z_left, z_right]) if not np.all(np.isfinite(bounds)): raise ValueError("kernel bounds are not finite; increase tail") @@ -227,8 +223,8 @@ def _translation_invariant_kernel_weights( m_max = int(np.clip(np.ceil(bounds[1] / spacing - 0.5), 0, max_offset)) m = np.arange(m_min, m_max + 1, dtype=np.int64) - lower = (m - 0.5) * spacing / width - upper = (m + 0.5) * spacing / width + lower = (m - 0.5) * spacing / scale + upper = (m + 0.5) * spacing / scale weights = np.maximum(dist.cdf(upper) - dist.cdf(lower), 0.0) return _trim_kernel_weights(m, weights) @@ -285,7 +281,7 @@ def _smooth_with_weights( def _smooth_relative_kernel_on_geomgrid( y: ArrayLike, log_spacing: float, - alpha: float, + scale: float, kernel: _Kernel, tail: float, ) -> _FloatArray: @@ -293,18 +289,18 @@ def _smooth_relative_kernel_on_geomgrid( if y.ndim != 1: raise ValueError("y must be one-dimensional") - if not np.isfinite(alpha) or alpha < 0: - raise ValueError("alpha must be non-negative") + if not np.isfinite(scale) or scale < 0: + raise ValueError("scale must be non-negative") if not np.isfinite(log_spacing) or log_spacing <= 0: raise ValueError("log_spacing must be positive") if not np.isfinite(tail) or not (0.0 < tail < 1.0): raise ValueError("tail must be between 0 and 1") - if y.size == 0 or alpha == 0: + if y.size == 0 or scale == 0: return y.copy() offsets, weights = _relative_kernel_weights( log_spacing=log_spacing, - alpha=alpha, + scale=scale, kernel=kernel, tail=tail, max_offset=y.size - 1, @@ -324,7 +320,7 @@ def _validate_max_grid_points(max_grid_points: int) -> None: def _smooth_relative_kernel( x: ArrayLike, y: ArrayLike, - alpha: float = 1.0, + scale: float, kernel: _Kernel = "gaussian", tail: float = 1e-12, max_grid_points: int = 1_000_000, @@ -336,8 +332,8 @@ def _smooth_relative_kernel( raise ValueError("x and y must be one-dimensional") if x.size != y.size: raise ValueError("x and y must have the same length") - if not np.isfinite(alpha) or alpha < 0: - raise ValueError("alpha must be non-negative") + if not np.isfinite(scale) or scale < 0: + raise ValueError("scale must be non-negative") if not np.isfinite(tail) or not (0.0 < tail < 1.0): raise ValueError("tail must be between 0 and 1") _validate_max_grid_points(max_grid_points) @@ -347,7 +343,7 @@ def _smooth_relative_kernel( raise ValueError("x must be positive") if np.any(np.diff(x) <= 0): raise ValueError("x must be strictly increasing") - if x.size == 0 or alpha == 0 or x.size == 1: + if x.size == 0 or scale == 0 or x.size == 1: return y.copy() logx = np.log(x) @@ -372,7 +368,7 @@ def _smooth_relative_kernel( yg = np.interp(xp, x, y) zg = _smooth_relative_kernel_on_geomgrid( yg, - alpha=alpha, + scale=scale, log_spacing=log_range / (k - 1), kernel=kernel, tail=tail, @@ -383,7 +379,7 @@ def _smooth_relative_kernel( def _smooth_kernel_values( x: ArrayLike, y: ArrayLike, - width: float, + scale: float, kernel: _Kernel = "gaussian", tail: float = 1e-12, max_grid_points: int = 1_000_000, @@ -395,8 +391,8 @@ def _smooth_kernel_values( raise ValueError("x and y must be one-dimensional") if x.size != y.size: raise ValueError("x and y must have the same length") - if not np.isfinite(width) or width < 0: - raise ValueError("width must be non-negative") + if not np.isfinite(scale) or scale < 0: + raise ValueError("scale must be non-negative") if not np.isfinite(tail) or not (0.0 < tail < 1.0): raise ValueError("tail must be between 0 and 1") _validate_max_grid_points(max_grid_points) @@ -404,7 +400,7 @@ def _smooth_kernel_values( raise ValueError("x must contain only finite values") if np.any(np.diff(x) <= 0): raise ValueError("x must be strictly increasing") - if x.size == 0 or width == 0 or x.size == 1: + if x.size == 0 or scale == 0 or x.size == 1: return y.copy() spacing = np.diff(x) @@ -426,7 +422,7 @@ def _smooth_kernel_values( offsets, weights = _translation_invariant_kernel_weights( spacing=x_range / (k - 1), - width=width, + scale=scale, kernel=kernel, tail=tail, max_offset=k - 1, @@ -486,21 +482,33 @@ def _scipp_output( return out -def _width_in_coordinate_unit(width: sc.Variable, x: sc.Variable) -> float: - if not isinstance(width, sc.Variable): - raise TypeError("width must be a scipp.Variable") - if width.ndim != 0: - raise sc.DimensionError("width must be a scalar") - if width.variances is not None: - raise sc.VariancesError("kernel widths with variances are not supported") - return float(width.to(unit=x.unit).value) +def _scale_in_coordinate_unit(scale: sc.Variable, x: sc.Variable) -> float: + if not isinstance(scale, sc.Variable): + raise TypeError("scale must be a scipp.Variable") + if scale.ndim != 0: + raise sc.DimensionError("scale must be a scalar") + if scale.variances is not None: + raise sc.VariancesError("kernel scales with variances are not supported") + return float(scale.to(unit=x.unit).value) -def smooth_kernel( +def _dimensionless_scale(scale: object) -> float: + if not isinstance(scale, sc.Variable): + if not isinstance(scale, Real): + raise TypeError("scale must be a real number or a scipp.Variable") + return float(scale) + if scale.ndim != 0: + raise sc.DimensionError("scale must be a scalar") + if scale.variances is not None: + raise sc.VariancesError("kernel scales with variances are not supported") + return float(scale.to(unit=sc.units.dimensionless).value) + + +def smooth( x: _ScippArray, y: sc.Variable | None = None, *, - width: sc.Variable, + scale: sc.Variable, kernel: _Kernel = "gaussian", tail: float = 1e-12, max_grid_points: int = 1_000_000, @@ -508,7 +516,7 @@ def smooth_kernel( """Smooth sampled data with a translation-invariant kernel. The kernel describes a distribution of displacements ``Z``, with displaced - coordinates given by ``x' = x + width * Z``. At the boundaries, the kernel + coordinates given by ``x' = x + scale * Z``. At the boundaries, the kernel is renormalized over the available finite input domain. Input that is not uniformly spaced is interpolated to a uniform grid, @@ -523,15 +531,17 @@ def smooth_kernel( Values to smooth when ``x`` contains the sample coordinates. Must be a one-dimensional variable with the same dimension as ``x``. Must be omitted when ``x`` is a data array. - width: + scale: Scale factor for the displacement distribution. Must be a scalar with a unit compatible with the coordinate. Set to zero to return a copy of the input without smoothing. kernel: - Kernel distribution. Supported names are ``'gaussian'``, ``'rectangle'``, - and ``'triangle'``, including their aliases. Alternatively, provide a - fully specified distribution with ``cdf``, ``ppf``, and ``support`` - methods. + Kernel distribution. The canonical names are ``'gaussian'``, + ``'boxcar'``, and ``'triangular'``. They represent a standard normal + distribution, a uniform distribution on [-1, 1], and a symmetric + triangular distribution on [-1, 1], respectively. Other aliases are + accepted. Alternatively, provide a fully specified distribution with + ``cdf``, ``ppf``, and ``support`` methods. tail: Total probability omitted when truncating a kernel with unbounded support. @@ -553,23 +563,23 @@ def smooth_kernel( grid exceeds ``max_grid_points``. scipp.DimensionError If the inputs are not one-dimensional, a pair of variables does not - have matching dimensions, or ``width`` is not scalar. + have matching dimensions, or ``scale`` is not scalar. scipp.CoordError If a data array has no dimension coordinate or has a bin-edge coordinate. scipp.UnitError - If the unit of ``width`` is incompatible with the coordinate unit. + If the unit of ``scale`` is incompatible with the coordinate unit. scipp.VariancesError - If the signal or ``width`` has variances. + If the signal or ``scale`` has variances. TypeError - If ``width`` is not a variable, ``kernel`` is not a distribution-like + If ``scale`` is not a variable, ``kernel`` is not a distribution-like object, or ``max_grid_points`` is not an integer. """ x, y, template = _scipp_input(x, y) values = _smooth_kernel_values( x.values, y.values, - width=_width_in_coordinate_unit(width, x), + scale=_scale_in_coordinate_unit(scale, x), kernel=kernel, tail=tail, max_grid_points=max_grid_points, @@ -577,148 +587,11 @@ def smooth_kernel( return _scipp_output(template, y, values) -def smooth_gaussian( - x: _ScippArray, - y: sc.Variable | None = None, - *, - width: sc.Variable, - tail: float = 1e-12, - max_grid_points: int = 1_000_000, -) -> _ScippArray: - """Smooth sampled data with a fixed-width Gaussian kernel. - - Parameters - ---------- - x: - One-dimensional data to smooth, or strictly increasing sample - coordinates. A data array must have a dimension coordinate. - y: - Values to smooth when ``x`` contains the sample coordinates. Must be a - one-dimensional variable with the same dimension as ``x``. Must be - omitted when ``x`` is a data array. - width: - Standard deviation of the Gaussian. Must be a scalar with a unit - compatible with the coordinate. - tail: - Total Gaussian probability omitted when truncating the kernel. - max_grid_points: - Maximum permitted size of the intermediate uniform grid. - - Returns - ------- - : - Smoothed data of the same type as the input. - - See Also - -------- - smooth_kernel: - Smooth with a named or user-provided translation-invariant kernel. - """ - return smooth_kernel( - x, - y, - width=width, - kernel="gaussian", - tail=tail, - max_grid_points=max_grid_points, - ) - - -def smooth_rectangle( - x: _ScippArray, - y: sc.Variable | None = None, - *, - width: sc.Variable, - max_grid_points: int = 1_000_000, -) -> _ScippArray: - """Smooth sampled data with a fixed-width rectangular kernel. - - The displacement is uniformly distributed on ``[-width, width]``. - - Parameters - ---------- - x: - One-dimensional data to smooth, or strictly increasing sample - coordinates. A data array must have a dimension coordinate. - y: - Values to smooth when ``x`` contains the sample coordinates. Must be a - one-dimensional variable with the same dimension as ``x``. Must be - omitted when ``x`` is a data array. - width: - Half-width of the rectangular kernel. Must be a scalar with a unit - compatible with the coordinate. - max_grid_points: - Maximum permitted size of the intermediate uniform grid. - - Returns - ------- - : - Smoothed data of the same type as the input. - - See Also - -------- - smooth_kernel: - Smooth with a named or user-provided translation-invariant kernel. - """ - return smooth_kernel( - x, - y, - width=width, - kernel="rectangle", - max_grid_points=max_grid_points, - ) - - -def smooth_triangle( +def smooth_relative( x: _ScippArray, y: sc.Variable | None = None, *, - width: sc.Variable, - max_grid_points: int = 1_000_000, -) -> _ScippArray: - """Smooth sampled data with a fixed-width triangular kernel. - - The displacement has symmetric triangular support on ``[-width, width]`` - and its peak at zero. - - Parameters - ---------- - x: - One-dimensional data to smooth, or strictly increasing sample - coordinates. A data array must have a dimension coordinate. - y: - Values to smooth when ``x`` contains the sample coordinates. Must be a - one-dimensional variable with the same dimension as ``x``. Must be - omitted when ``x`` is a data array. - width: - Half-width of the triangular kernel. Must be a scalar with a unit - compatible with the coordinate. - max_grid_points: - Maximum permitted size of the intermediate uniform grid. - - Returns - ------- - : - Smoothed data of the same type as the input. - - See Also - -------- - smooth_kernel: - Smooth with a named or user-provided translation-invariant kernel. - """ - return smooth_kernel( - x, - y, - width=width, - kernel="triangle", - max_grid_points=max_grid_points, - ) - - -def smooth_relative_kernel( - x: _ScippArray, - y: sc.Variable | None = None, - alpha: float = 1.0, + scale: float | sc.Variable, kernel: _Kernel = "gaussian", tail: float = 1e-12, max_grid_points: int = 1_000_000, @@ -726,7 +599,7 @@ def smooth_relative_kernel( """Smooth sampled data with a kernel of relative width. The kernel describes a distribution of relative displacements ``Z``, with - displaced coordinates given by ``x' = x * (1 + alpha * Z)``. At the + displaced coordinates given by ``x' = x * (1 + scale * Z)``. At the boundaries, the kernel is renormalized over the available finite input domain. @@ -743,14 +616,17 @@ def smooth_relative_kernel( Values to smooth when ``x`` contains the sample coordinates. Must be a one-dimensional variable with the same dimension as ``x``. Must be omitted when ``x`` is a data array. - alpha: - Scale factor for the relative-displacement distribution. Set to zero to + scale: + Dimensionless scale factor for the relative-displacement distribution. + May be a real number or a scalar, dimensionless variable. Set to zero to return a copy of the input without smoothing. kernel: - Kernel distribution. Supported names are ``'gaussian'``, ``'rectangle'``, - and ``'triangle'``, including their aliases. Alternatively, provide a - fully specified distribution with ``cdf``, ``ppf``, and ``support`` - methods. + Kernel distribution. The canonical names are ``'gaussian'``, + ``'boxcar'``, and ``'triangular'``. They represent a standard normal + distribution, a uniform distribution on [-1, 1], and a symmetric + triangular distribution on [-1, 1], respectively. Other aliases are + accepted. Alternatively, provide a fully specified distribution with + ``cdf``, ``ppf``, and ``support`` methods. tail: Total probability omitted when truncating a kernel with unbounded support or support reaching the nonpositive coordinate domain. @@ -772,172 +648,25 @@ def smooth_relative_kernel( intermediate grid exceeds ``max_grid_points``. scipp.DimensionError If the inputs are not one-dimensional or a pair of variables does not - have matching dimensions. + have matching dimensions, or ``scale`` is not scalar. scipp.CoordError If a data array has no dimension coordinate or has a bin-edge coordinate. + scipp.UnitError + If ``scale`` is a variable with a non-dimensionless unit. scipp.VariancesError - If the signal has variances. + If the signal or ``scale`` has variances. TypeError - If ``kernel`` is not a distribution-like object, or if - ``max_grid_points`` is not an integer. + If ``scale`` is neither a real number nor a variable, ``kernel`` is not + a distribution-like object, or ``max_grid_points`` is not an integer. """ x, y, template = _scipp_input(x, y) values = _smooth_relative_kernel( x.values, y.values, - alpha=alpha, + scale=_dimensionless_scale(scale), kernel=kernel, tail=tail, max_grid_points=max_grid_points, ) return _scipp_output(template, y, values) - - -def smooth_relative_gaussian( - x: _ScippArray, - y: sc.Variable | None = None, - alpha: float = 1.0, - tail: float = 1e-12, - max_grid_points: int = 1_000_000, -) -> _ScippArray: - """Smooth sampled data with a relative Gaussian kernel. - - ``alpha`` is the standard deviation of the Gaussian as a fraction of each - coordinate value. - - Parameters - ---------- - x: - One-dimensional data to smooth, or positive, strictly increasing - one-dimensional sample coordinates. A data array must have a - dimension coordinate. - y: - Values to smooth when ``x`` contains the sample coordinates. Must be a - one-dimensional variable with the same dimension as ``x``. Must be - omitted when ``x`` is a data array. - alpha: - Relative standard deviation of the Gaussian kernel. - tail: - Total Gaussian probability omitted when truncating the kernel. - max_grid_points: - Maximum permitted size of the intermediate geometric grid. - - Returns - ------- - : - Smoothed data of the same type as the input. - - See Also - -------- - smooth_relative_kernel: - Smooth with a named or user-provided relative kernel. - """ - return smooth_relative_kernel( - x, - y, - alpha=alpha, - kernel="gaussian", - tail=tail, - max_grid_points=max_grid_points, - ) - - -def smooth_relative_rectangle( - x: _ScippArray, - y: sc.Variable | None = None, - alpha: float = 1.0, - tail: float = 1e-12, - max_grid_points: int = 1_000_000, -) -> _ScippArray: - """Smooth sampled data with a relative rectangular kernel. - - The relative displacement is uniformly distributed on - ``[-alpha, alpha]``. - - Parameters - ---------- - x: - One-dimensional data to smooth, or positive, strictly increasing - one-dimensional sample coordinates. A data array must have a - dimension coordinate. - y: - Values to smooth when ``x`` contains the sample coordinates. Must be a - one-dimensional variable with the same dimension as ``x``. Must be - omitted when ``x`` is a data array. - alpha: - Half-width of the rectangular kernel relative to each coordinate value. - tail: - Total probability omitted if the kernel reaches the nonpositive - coordinate domain. - max_grid_points: - Maximum permitted size of the intermediate geometric grid. - - Returns - ------- - : - Smoothed data of the same type as the input. - - See Also - -------- - smooth_relative_kernel: - Smooth with a named or user-provided relative kernel. - """ - return smooth_relative_kernel( - x, - y, - alpha=alpha, - kernel="rectangle", - tail=tail, - max_grid_points=max_grid_points, - ) - - -def smooth_relative_triangle( - x: _ScippArray, - y: sc.Variable | None = None, - alpha: float = 1.0, - tail: float = 1e-12, - max_grid_points: int = 1_000_000, -) -> _ScippArray: - """Smooth sampled data with a relative triangular kernel. - - The relative displacement has symmetric triangular support on - ``[-alpha, alpha]`` and its peak at zero. - - Parameters - ---------- - x: - One-dimensional data to smooth, or positive, strictly increasing - one-dimensional sample coordinates. A data array must have a - dimension coordinate. - y: - Values to smooth when ``x`` contains the sample coordinates. Must be a - one-dimensional variable with the same dimension as ``x``. Must be - omitted when ``x`` is a data array. - alpha: - Half-width of the triangular kernel relative to each coordinate value. - tail: - Total probability omitted if the kernel reaches the nonpositive - coordinate domain. - max_grid_points: - Maximum permitted size of the intermediate geometric grid. - - Returns - ------- - : - Smoothed data of the same type as the input. - - See Also - -------- - smooth_relative_kernel: - Smooth with a named or user-provided relative kernel. - """ - return smooth_relative_kernel( - x, - y, - alpha=alpha, - kernel="triangle", - tail=tail, - max_grid_points=max_grid_points, - ) diff --git a/tests/smoothing_test.py b/tests/smoothing_test.py index e7d454a40..e37004bae 100644 --- a/tests/smoothing_test.py +++ b/tests/smoothing_test.py @@ -11,14 +11,8 @@ from scippneutron.smoothing import ( _relative_kernel_weights, _smooth_relative_kernel_on_geomgrid, - smooth_gaussian, - smooth_kernel, - smooth_rectangle, - smooth_relative_gaussian, - smooth_relative_kernel, - smooth_relative_rectangle, - smooth_relative_triangle, - smooth_triangle, + smooth, + smooth_relative, ) @@ -87,10 +81,10 @@ def _quadratic_smoothing_error(*, size, alpha, tail=1e-9): upper = 0.9 x = _geometric_cell_centers(lower, upper, size) actual = _smooth_values( - smooth_relative_gaussian, + smooth_relative, x, _quadratic(x), - alpha=alpha, + scale=alpha, tail=tail, ) expected = _exact_smoothed_quadratic(x, alpha, lower, upper) @@ -102,10 +96,10 @@ def _fixed_width_quadratic_smoothing_error(*, size, width, tail=1e-9): upper = 0.9 x = _linear_cell_centers(lower, upper, size) x_var, y_var = _variables(x, _quadratic(x)) - actual = smooth_gaussian( + actual = smooth( x_var, y_var, - width=sc.scalar(width, unit='m'), + scale=sc.scalar(width, unit='m'), tail=tail, ).values expected = _exact_smoothed_quadratic_with_sigma(x, width, lower, upper) @@ -168,14 +162,14 @@ def test_fixed_width_gaussian_error_scales_with_tail_until_grid_error_dominates( @pytest.mark.parametrize( - ("smooth", "kernel_variance", "max_error"), + ("kernel", "kernel_variance", "max_error"), [ - (smooth_rectangle, 1.0 / 3.0, 3e-7), - (smooth_triangle, 1.0 / 6.0, 5e-8), + ("boxcar", 1.0 / 3.0, 3e-7), + ("triangular", 1.0 / 6.0, 5e-8), ], ) def test_fixed_width_compact_kernel_matches_interior_quadratic_moments( - smooth, kernel_variance, max_error + kernel, kernel_variance, max_error ): lower = -0.4 upper = 0.9 @@ -184,19 +178,24 @@ def test_fixed_width_compact_kernel_matches_interior_quadratic_moments( x = _linear_cell_centers(lower, upper, size) x_var, y_var = _variables(x, _quadratic(x)) - actual = smooth(x_var, y_var, width=sc.scalar(width, unit='m')).values + actual = smooth( + x_var, + y_var, + scale=sc.scalar(width, unit='m'), + kernel=kernel, + ).values expected = 1.0 + 0.3 * x + 0.7 * (x**2 + width**2 * kernel_variance) interior = (x - width >= lower) & (x + width <= upper) assert np.max(np.abs(actual[interior] - expected[interior])) < max_error -def test_fixed_width_smoothing_converts_width_to_coordinate_unit(): +def test_fixed_width_smoothing_converts_scale_to_coordinate_unit(): x = sc.linspace('x', -0.4, 0.9, 100, unit='m') y = sc.array(dims=['x'], values=_quadratic(x.values), unit='counts') - in_meters = smooth_gaussian(x, y, width=sc.scalar(0.1, unit='m')) - in_centimeters = smooth_gaussian(x, y, width=sc.scalar(10.0, unit='cm')) + in_meters = smooth(x, y, scale=sc.scalar(0.1, unit='m')) + in_centimeters = smooth(x, y, scale=sc.scalar(10.0, unit='cm')) assert sc.identical(in_meters, in_centimeters) @@ -204,15 +203,15 @@ def test_fixed_width_smoothing_converts_width_to_coordinate_unit(): def test_fixed_width_smoothing_accepts_distribution(): x = sc.linspace('x', -0.4, 0.9, 100, unit='m') y = sc.array(dims=['x'], values=_quadratic(x.values), unit='counts') - width = sc.scalar(0.1, unit='m') + scale = sc.scalar(0.1, unit='m') - actual = smooth_kernel( + actual = smooth( x, y, - width=width, + scale=scale, kernel=uniform(loc=-1.0, scale=2.0), ) - expected = smooth_rectangle(x, y, width=width) + expected = smooth(x, y, scale=scale, kernel="boxcar") assert sc.identical(actual, expected) @@ -221,10 +220,10 @@ def test_fixed_width_asymmetric_kernel_has_correct_direction(): x = sc.arange('x', 0.0, 2.0, 0.1, unit='m') y = sc.array(dims=['x'], values=x.values, unit='counts') - actual = smooth_kernel( + actual = smooth( x, y, - width=sc.scalar(0.2, unit='m'), + scale=sc.scalar(0.2, unit='m'), kernel=uniform(loc=0.5, scale=1.0), ) @@ -238,7 +237,7 @@ def test_nonfinite_value_only_affects_overlapping_kernel_windows(): values[size // 2] = np.nan y = sc.array(dims=['x'], values=values) - actual = smooth_gaussian(x, y, width=sc.scalar(150.0)) + actual = smooth(x, y, scale=sc.scalar(150.0)) assert np.isfinite(actual.values[0]) assert np.isnan(actual.values[size // 2]) @@ -252,8 +251,8 @@ def test_fixed_width_smoothing_accepts_data_array(): coords={'x': x, 'aux': sc.arange('x', 100)}, ) - actual = smooth_gaussian(data, width=sc.scalar(0.1, unit='m')) - expected = smooth_gaussian(x, data.data, width=sc.scalar(0.1, unit='m')) + actual = smooth(data, scale=sc.scalar(0.1, unit='m')) + expected = smooth(x, data.data, scale=sc.scalar(0.1, unit='m')) assert sc.identical(actual.data, expected) assert sc.identical(actual.coords['x'], data.coords['x']) @@ -263,14 +262,14 @@ def test_fixed_width_smoothing_accepts_data_array(): def test_fixed_width_smoothing_resamples_nonuniform_grid(): x_values = np.array([0.0, 0.25, 0.75, 1.0]) y_values = _quadratic(x_values) - width = sc.scalar(0.2, unit='m') + scale = sc.scalar(0.2, unit='m') x, y = _variables(x_values, y_values) - actual = smooth_gaussian(x, y, width=width) + actual = smooth(x, y, scale=scale) xp_values = np.linspace(0.0, 1.0, 5) xp, yp = _variables(xp_values, np.interp(xp_values, x_values, y_values)) - smoothed_yp = smooth_gaussian(xp, yp, width=width) + smoothed_yp = smooth(xp, yp, scale=scale) expected = np.interp(x_values, xp_values, smoothed_yp.values) np.testing.assert_allclose(actual.values, expected) @@ -285,10 +284,10 @@ def test_fixed_width_gaussian_matches_exact_quadratic_on_nonuniform_grid(): x += 0.2 * spacing * np.sin(np.linspace(0.0, 8.0 * np.pi, size)) x_var, y_var = _variables(x, _quadratic(x)) - actual = smooth_gaussian( + actual = smooth( x_var, y_var, - width=sc.scalar(width, unit='m'), + scale=sc.scalar(width, unit='m'), tail=1e-9, ).values @@ -311,30 +310,22 @@ def test_fixed_width_smoothing_rejects_resampled_grid_larger_than_limit(): ValueError, match=r"uniform resampling would require too many points.*max_grid_points=4", ): - smooth_kernel(x, y, width=sc.scalar(0.1, unit='m'), max_grid_points=4) + smooth(x, y, scale=sc.scalar(0.1, unit='m'), max_grid_points=4) def test_fixed_width_smoothing_accepts_resampled_grid_equal_to_limit(): x, y = _variables([0.0, 0.25, 1.0], [1.0, 1.0, 1.0]) - actual = smooth_kernel( + actual = smooth( x, y, - width=sc.scalar(0.1, unit='m'), + scale=sc.scalar(0.1, unit='m'), max_grid_points=5, ) np.testing.assert_allclose(actual.values, y.values) -@pytest.mark.parametrize("smooth", [smooth_gaussian, smooth_rectangle, smooth_triangle]) -def test_fixed_width_convenience_functions_forward_max_grid_points(smooth): - x, y = _variables([0.0, 0.25, 1.0], [1.0, 2.0, 3.0]) - - with pytest.raises(ValueError, match="max_grid_points=2"): - smooth(x, y, width=sc.scalar(0.1, unit='m'), max_grid_points=2) - - def test_fixed_width_pathologically_close_coordinates_fail_before_allocation(): x, y = _variables( [1.0, np.nextafter(1.0, 2.0), 2.0], @@ -342,46 +333,46 @@ def test_fixed_width_pathologically_close_coordinates_fail_before_allocation(): ) with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"): - smooth_kernel(x, y, width=sc.scalar(0.1, unit='m')) + smooth(x, y, scale=sc.scalar(0.1, unit='m')) -@pytest.mark.parametrize("width", [np.nan, np.inf, -np.inf, -1.0]) -def test_fixed_width_smoothing_rejects_invalid_width(width): +@pytest.mark.parametrize("scale", [np.nan, np.inf, -np.inf, -1.0]) +def test_fixed_width_smoothing_rejects_invalid_scale(scale): x, y = _variables([], []) - with pytest.raises(ValueError, match="width must be non-negative"): - smooth_kernel(x, y, width=sc.scalar(width, unit='m')) + with pytest.raises(ValueError, match="scale must be non-negative"): + smooth(x, y, scale=sc.scalar(scale, unit='m')) -def test_fixed_width_smoothing_rejects_non_scalar_width(): +def test_fixed_width_smoothing_rejects_non_scalar_scale(): x, y = _variables([0.0, 1.0], [1.0, 2.0]) - with pytest.raises(sc.DimensionError, match="width must be a scalar"): - smooth_kernel(x, y, width=sc.array(dims=['width'], values=[0.1], unit='m')) + with pytest.raises(sc.DimensionError, match="scale must be a scalar"): + smooth(x, y, scale=sc.array(dims=['scale'], values=[0.1], unit='m')) -def test_fixed_width_smoothing_rejects_non_variable_width(): +def test_fixed_width_smoothing_rejects_non_variable_scale(): x, y = _variables([0.0, 1.0], [1.0, 2.0]) - with pytest.raises(TypeError, match=r"width must be a scipp\.Variable"): - smooth_kernel(x, y, width=0.1) # type: ignore[arg-type] + with pytest.raises(TypeError, match=r"scale must be a scipp\.Variable"): + smooth(x, y, scale=0.1) # type: ignore[arg-type] -def test_fixed_width_smoothing_rejects_incompatible_width_unit(): +def test_fixed_width_smoothing_rejects_incompatible_scale_unit(): x, y = _variables([0.0, 1.0], [1.0, 2.0]) with pytest.raises(sc.UnitError): - smooth_kernel(x, y, width=sc.scalar(0.1, unit='s')) + smooth(x, y, scale=sc.scalar(0.1, unit='s')) -def test_fixed_width_smoothing_rejects_width_with_variance(): +def test_fixed_width_smoothing_rejects_scale_with_variance(): x, y = _variables([0.0, 1.0], [1.0, 2.0]) - with pytest.raises(sc.VariancesError, match="widths with variances"): - smooth_kernel( + with pytest.raises(sc.VariancesError, match="scales with variances"): + smooth( x, y, - width=sc.scalar(0.1, variance=0.01, unit='m'), + scale=sc.scalar(0.1, variance=0.01, unit='m'), ) @@ -459,14 +450,14 @@ def test_gaussian_smoothing_error_scales_with_tail_until_grid_error_dominates(): @pytest.mark.parametrize( - ("smooth", "relative_variance", "max_error"), + ("kernel", "relative_variance", "max_error"), [ - (smooth_relative_rectangle, 1.0 / 3.0, 3e-7), - (smooth_relative_triangle, 1.0 / 6.0, 5e-8), + ("boxcar", 1.0 / 3.0, 3e-7), + ("triangular", 1.0 / 6.0, 5e-8), ], ) def test_compact_symmetric_kernel_matches_exact_interior_quadratic_moments( - smooth, relative_variance, max_error + kernel, relative_variance, max_error ): lower = 0.1 upper = 0.9 @@ -474,7 +465,13 @@ def test_compact_symmetric_kernel_matches_exact_interior_quadratic_moments( alpha = 0.1 x = _geometric_cell_centers(lower, upper, size) - actual = _smooth_values(smooth, x, _quadratic(x), alpha=alpha) + actual = _smooth_values( + smooth_relative, + x, + _quadratic(x), + scale=alpha, + kernel=kernel, + ) expected = 1.0 + 0.3 * x + 0.7 * x**2 * (1.0 + alpha**2 * relative_variance) interior = (x * (1.0 - alpha) >= lower) & (x * (1.0 + alpha) <= upper) @@ -482,7 +479,7 @@ def test_compact_symmetric_kernel_matches_exact_interior_quadratic_moments( def test_smoothing_module_is_exposed_by_package(): - assert scn.smoothing.smooth_relative_gaussian is smooth_relative_gaussian + assert scn.smoothing.smooth_relative is smooth_relative def test_accepts_data_array_and_preserves_metadata(): @@ -497,8 +494,8 @@ def test_accepts_data_array_and_preserves_metadata(): ) original = data.copy() - actual = smooth_relative_gaussian(data, alpha=0.1) - expected = smooth_relative_gaussian(x, data.data, alpha=0.1) + actual = smooth_relative(data, scale=0.1) + expected = smooth_relative(x, data.data, scale=0.1) assert sc.identical(data, original) assert isinstance(actual, sc.DataArray) @@ -511,6 +508,36 @@ def test_accepts_data_array_and_preserves_metadata(): assert sc.identical(data, original) +def test_relative_smoothing_accepts_dimensionless_variable_scale(): + x, y = _variables([1.0, 2.0, 3.0], [1.0, 4.0, 9.0]) + + actual = smooth_relative(x, y, scale=sc.scalar(0.1)) + expected = smooth_relative(x, y, scale=0.1) + + assert sc.identical(actual, expected) + + +def test_relative_smoothing_rejects_non_scalar_scale(): + x, y = _variables([1.0, 2.0], [1.0, 2.0]) + + with pytest.raises(sc.DimensionError, match="scale must be a scalar"): + smooth_relative(x, y, scale=sc.array(dims=['scale'], values=[0.1])) + + +def test_relative_smoothing_rejects_scale_with_unit(): + x, y = _variables([1.0, 2.0], [1.0, 2.0]) + + with pytest.raises(sc.UnitError): + smooth_relative(x, y, scale=sc.scalar(0.1, unit='m')) + + +def test_relative_smoothing_rejects_scale_with_variance(): + x, y = _variables([1.0, 2.0], [1.0, 2.0]) + + with pytest.raises(sc.VariancesError, match="scales with variances"): + smooth_relative(x, y, scale=sc.scalar(0.1, variance=0.01)) + + def test_rejects_variable_with_variances(): x = sc.geomspace('x', 0.1, 0.9, 100, unit='m') values = _quadratic(x.values) @@ -522,7 +549,7 @@ def test_rejects_variable_with_variances(): ) with pytest.raises(sc.VariancesError, match="signals with variances"): - smooth_relative_gaussian(x, y, alpha=0.1) + smooth_relative(x, y, scale=0.1) def test_rejects_data_array_with_variances(): @@ -538,12 +565,12 @@ def test_rejects_data_array_with_variances(): ) with pytest.raises(sc.VariancesError, match="signals with variances"): - smooth_relative_gaussian(data, alpha=0.1) + smooth_relative(data, scale=0.1) def test_rejects_numpy_arrays(): with pytest.raises(TypeError, match="DataArray or a pair of Variables"): - smooth_relative_gaussian(np.arange(1.0, 4.0), np.ones(3)) + smooth_relative(np.arange(1.0, 4.0), np.ones(3), scale=0.1) def test_rejects_variables_with_different_dimensions(): @@ -551,14 +578,14 @@ def test_rejects_variables_with_different_dimensions(): y = sc.ones(dims=['y'], shape=[3]) with pytest.raises(sc.DimensionError, match="same dimension"): - smooth_relative_gaussian(x, y) + smooth_relative(x, y, scale=0.1) def test_rejects_data_array_without_dimension_coordinate(): data = sc.DataArray(sc.ones(dims=['x'], shape=[3])) with pytest.raises(sc.CoordError, match="dimension coordinate"): - smooth_relative_gaussian(data) + smooth_relative(data, scale=0.1) def test_rejects_data_array_with_bin_edge_coordinate(): @@ -568,7 +595,7 @@ def test_rejects_data_array_with_bin_edge_coordinate(): ) with pytest.raises(sc.CoordError, match="bin edges"): - smooth_relative_gaussian(data) + smooth_relative(data, scale=0.1) def test_rejects_data_array_with_masks(): @@ -579,7 +606,7 @@ def test_rejects_data_array_with_masks(): ) with pytest.raises(ValueError, match="data with masks"): - smooth_relative_gaussian(data) + smooth_relative(data, scale=0.1) def test_rejects_y_with_data_array(): @@ -587,7 +614,7 @@ def test_rejects_y_with_data_array(): data = sc.DataArray(y, coords={'x': x}) with pytest.raises(TypeError, match="y must be omitted"): - smooth_relative_gaussian(data, y) + smooth_relative(data, y, scale=0.1) def test_rejects_geometric_grid_larger_than_limit(): @@ -598,14 +625,14 @@ def test_rejects_geometric_grid_larger_than_limit(): ValueError, match=r"geometric resampling would require too many points.*max_grid_points=70", ): - _smooth_values(smooth_relative_kernel, x, y, max_grid_points=70) + _smooth_values(smooth_relative, x, y, scale=0.1, max_grid_points=70) def test_accepts_geometric_grid_equal_to_limit(): x = np.array([1.0, 1.01, 2.0]) y = np.ones_like(x) - actual = _smooth_values(smooth_relative_kernel, x, y, max_grid_points=71) + actual = _smooth_values(smooth_relative, x, y, scale=0.1, max_grid_points=71) np.testing.assert_allclose(actual, y) @@ -615,40 +642,30 @@ def test_geometric_input_does_not_gain_a_point_from_roundoff(): x = _geometric_cell_centers(0.1, 0.9, size) y = _quadratic(x) - actual = _smooth_values(smooth_relative_gaussian, x, y, max_grid_points=size) + actual = _smooth_values( + smooth_relative, + x, + y, + scale=0.1, + max_grid_points=size, + ) assert actual.shape == y.shape -@pytest.mark.parametrize( - "smooth", - [ - smooth_relative_gaussian, - smooth_relative_rectangle, - smooth_relative_triangle, - ], -) -def test_convenience_functions_forward_max_grid_points(smooth): - x = np.array([1.0, 1.01, 2.0]) - y = np.ones_like(x) - - with pytest.raises(ValueError, match="max_grid_points=2"): - _smooth_values(smooth, x, y, max_grid_points=2) - - def test_pathologically_close_coordinates_fail_before_allocation(): x = np.array([1.0, np.nextafter(1.0, 2.0), 2.0]) y = np.ones_like(x) with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"): - _smooth_values(smooth_relative_kernel, x, y) + _smooth_values(smooth_relative, x, y, scale=0.1) def test_wide_coordinate_range_does_not_overflow_grid_construction(): x = np.array([1e-300, 1.0, 1e300]) y = np.ones_like(x) - actual = _smooth_values(smooth_relative_kernel, x, y) + actual = _smooth_values(smooth_relative, x, y, scale=0.1) np.testing.assert_allclose(actual, y) @@ -658,7 +675,7 @@ def test_kernel_stencil_is_bounded_before_allocation(): offsets, weights = _relative_kernel_weights( log_spacing=1e-6, - alpha=1.0, + scale=1.0, kernel="gaussian", tail=1e-12, max_offset=max_offset, @@ -679,7 +696,7 @@ def test_asymmetric_kernel_convolution_matches_direct_weighted_sum(): offsets, weights = _relative_kernel_weights( log_spacing=log_spacing, - alpha=alpha, + scale=alpha, kernel=kernel, tail=1e-12, max_offset=size - 1, @@ -697,7 +714,7 @@ def test_asymmetric_kernel_convolution_matches_direct_weighted_sum(): actual = _smooth_relative_kernel_on_geomgrid( y, log_spacing=log_spacing, - alpha=alpha, + scale=alpha, kernel=kernel, tail=1e-12, ) @@ -709,7 +726,7 @@ def test_kernel_with_no_reachable_mass_returns_nan(): actual = _smooth_relative_kernel_on_geomgrid( np.arange(5.0), log_spacing=0.1, - alpha=0.1, + scale=0.1, kernel=uniform(loc=100.0, scale=1.0), tail=1e-12, ) @@ -717,18 +734,18 @@ def test_kernel_with_no_reachable_mass_returns_nan(): assert np.all(np.isnan(actual)) -@pytest.mark.parametrize("alpha", [np.nan, np.inf, -np.inf, -1.0]) -def test_rejects_invalid_alpha_before_noop_return(alpha): +@pytest.mark.parametrize("scale", [np.nan, np.inf, -np.inf, -1.0]) +def test_rejects_invalid_relative_scale_before_noop_return(scale): x, y = _variables([], []) - with pytest.raises(ValueError, match="alpha must be non-negative"): - smooth_relative_kernel(x, y, alpha=alpha) + with pytest.raises(ValueError, match="scale must be non-negative"): + smooth_relative(x, y, scale=scale) @pytest.mark.parametrize("tail", [np.nan, np.inf, -np.inf, 0.0, 1.0]) def test_rejects_invalid_tail_before_noop_return(tail): x, y = _variables([], []) with pytest.raises(ValueError, match="tail must be between 0 and 1"): - smooth_relative_kernel(x, y, tail=tail) + smooth_relative(x, y, scale=0.1, tail=tail) @pytest.mark.parametrize( @@ -742,32 +759,32 @@ def test_rejects_invalid_tail_before_noop_return(tail): def test_rejects_invalid_single_coordinate_before_noop_return(x, message): x, y = _variables(x, [1.0]) with pytest.raises(ValueError, match=message): - smooth_relative_kernel(x, y) + smooth_relative(x, y, scale=0.1) @pytest.mark.parametrize( - ("smooth", "kwargs"), + ("function", "kwargs"), [ - (smooth_relative_kernel, {}), - (smooth_kernel, {"width": sc.scalar(0.1, unit='m')}), + (smooth_relative, {"scale": 0.1}), + (smooth, {"scale": sc.scalar(0.1, unit='m')}), ], ) @pytest.mark.parametrize("max_grid_points", [True, 2.5]) -def test_rejects_non_integer_max_grid_points(smooth, kwargs, max_grid_points): +def test_rejects_non_integer_max_grid_points(function, kwargs, max_grid_points): x, y = _variables([], []) with pytest.raises(TypeError, match="max_grid_points must be an integer"): - smooth(x, y, max_grid_points=max_grid_points, **kwargs) + function(x, y, max_grid_points=max_grid_points, **kwargs) @pytest.mark.parametrize( - ("smooth", "kwargs"), + ("function", "kwargs"), [ - (smooth_relative_kernel, {}), - (smooth_kernel, {"width": sc.scalar(0.1, unit='m')}), + (smooth_relative, {"scale": 0.1}), + (smooth, {"scale": sc.scalar(0.1, unit='m')}), ], ) @pytest.mark.parametrize("max_grid_points", [-1, 0, 1]) -def test_rejects_too_small_max_grid_points(smooth, kwargs, max_grid_points): +def test_rejects_too_small_max_grid_points(function, kwargs, max_grid_points): x, y = _variables([], []) with pytest.raises(ValueError, match="max_grid_points must be at least 2"): - smooth(x, y, max_grid_points=max_grid_points, **kwargs) + function(x, y, max_grid_points=max_grid_points, **kwargs) From 7a820602afc5c26b8767c4a7aab00fee7a7c7d76 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 11 Aug 2026 11:14:34 +0200 Subject: [PATCH 5/6] fix: allow large inputs + convert grid size to int after validation + remove unnecessary validation --- src/scippneutron/smoothing.py | 99 +++++++++++++---------------------- tests/smoothing_test.py | 86 ++++++++++++++++++++++-------- 2 files changed, 98 insertions(+), 87 deletions(-) diff --git a/src/scippneutron/smoothing.py b/src/scippneutron/smoothing.py index adb8cec8d..55f5887c2 100644 --- a/src/scippneutron/smoothing.py +++ b/src/scippneutron/smoothing.py @@ -1,5 +1,16 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +"""Kernel smoothing of one-dimensional sampled signals. + +``smooth`` applies a translation-invariant displacement distribution, for a kernel +whose width is constant in coordinate units. ``smooth_relative`` instead applies a +distribution of fractional displacements, so its width scales with the coordinate. + +Both procedures accept irregularly spaced samples. They interpolate onto a uniform +or geometric working grid that is no coarser than the tightest input spacing, +integrate the kernel probability over grid cells, renormalize the convolution at +finite boundaries, and interpolate the result back to the original coordinates. +""" from __future__ import annotations @@ -116,15 +127,6 @@ def _relative_kernel_weights( where f is the PDF of the supplied distribution. """ - if not np.isfinite(scale) or scale <= 0: - raise ValueError("scale must be positive") - if not np.isfinite(log_spacing) or log_spacing <= 0: - raise ValueError("log_spacing must be positive") - if not np.isfinite(tail) or not (0.0 < tail < 1.0): - raise ValueError("tail must be between 0 and 1") - if max_offset < 0: - raise ValueError("max_offset must be non-negative") - dist = _as_kernel_distribution(kernel) # Positive physical domain: @@ -199,15 +201,6 @@ def _translation_invariant_kernel_weights( tail: float, max_offset: int, ) -> tuple[_IntArray, _FloatArray]: - if not np.isfinite(scale) or scale <= 0: - raise ValueError("scale must be positive") - if not np.isfinite(spacing) or spacing <= 0: - raise ValueError("spacing must be positive") - if not np.isfinite(tail) or not (0.0 < tail < 1.0): - raise ValueError("tail must be between 0 and 1") - if max_offset < 0: - raise ValueError("max_offset must be non-negative") - dist = _as_kernel_distribution(kernel) z_left, z_right = _kernel_support(dist) if not np.all(np.isfinite([z_left, z_right])): @@ -278,36 +271,6 @@ def _smooth_with_weights( return out -def _smooth_relative_kernel_on_geomgrid( - y: ArrayLike, - log_spacing: float, - scale: float, - kernel: _Kernel, - tail: float, -) -> _FloatArray: - y = np.asarray(y, dtype=float) - - if y.ndim != 1: - raise ValueError("y must be one-dimensional") - if not np.isfinite(scale) or scale < 0: - raise ValueError("scale must be non-negative") - if not np.isfinite(log_spacing) or log_spacing <= 0: - raise ValueError("log_spacing must be positive") - if not np.isfinite(tail) or not (0.0 < tail < 1.0): - raise ValueError("tail must be between 0 and 1") - if y.size == 0 or scale == 0: - return y.copy() - - offsets, weights = _relative_kernel_weights( - log_spacing=log_spacing, - scale=scale, - kernel=kernel, - tail=tail, - max_offset=y.size - 1, - ) - return _smooth_with_weights(y, offsets, weights) - - def _validate_max_grid_points(max_grid_points: int) -> None: if isinstance(max_grid_points, bool | np.bool_) or not isinstance( max_grid_points, int | np.integer @@ -317,7 +280,7 @@ def _validate_max_grid_points(max_grid_points: int) -> None: raise ValueError("max_grid_points must be at least 2") -def _smooth_relative_kernel( +def _smooth_relative_values( x: ArrayLike, y: ArrayLike, scale: float, @@ -355,28 +318,33 @@ def _smooth_relative_kernel( if np.allclose(dlog, dlog[0], rtol=1e-7, atol=0.0): k = x.size else: - k = int(np.ceil(log_range / np.min(dlog))) + 1 + k = np.ceil(log_range / np.min(dlog)) + 1.0 - if k > max_grid_points: + if not (k <= max_grid_points) and k > 2 * x.size: raise ValueError( "geometric resampling would require too many points, exceeding " f"max_grid_points={max_grid_points:,}. Increase max_grid_points to " "allow a larger grid." ) + k = int(k) + log_spacing = log_range / (k - 1) + if not np.isfinite(log_spacing) or log_spacing <= 0: + raise ValueError("log_spacing must be positive") xp = np.geomspace(x[0], x[-1], k) yg = np.interp(xp, x, y) - zg = _smooth_relative_kernel_on_geomgrid( - yg, + offsets, weights = _relative_kernel_weights( scale=scale, - log_spacing=log_range / (k - 1), + log_spacing=log_spacing, kernel=kernel, tail=tail, + max_offset=k - 1, ) + zg = _smooth_with_weights(yg, offsets, weights) return cast(_FloatArray, np.interp(x, xp, zg)) -def _smooth_kernel_values( +def _smooth_values( x: ArrayLike, y: ArrayLike, scale: float, @@ -413,15 +381,18 @@ def _smooth_kernel_values( else: k = int(np.ceil(x_range / np.min(spacing))) + 1 - if k > max_grid_points: + if k > max_grid_points and k > 2 * x.size: raise ValueError( "uniform resampling would require too many points, exceeding " f"max_grid_points={max_grid_points:,}. Increase max_grid_points to " "allow a larger grid." ) + spacing = x_range / (k - 1) + if not np.isfinite(spacing) or spacing <= 0: + raise ValueError("spacing must be positive") offsets, weights = _translation_invariant_kernel_weights( - spacing=x_range / (k - 1), + spacing=spacing, scale=scale, kernel=kernel, tail=tail, @@ -452,7 +423,7 @@ def _scipp_input( template = x y = x.data x = x.coords[x.dim] - elif not isinstance(y, sc.Variable): + elif not isinstance(x, sc.Variable) or not isinstance(y, sc.Variable): raise TypeError("expected a DataArray or a pair of Variables") if x.ndim != 1 or y.ndim != 1: @@ -546,8 +517,8 @@ def smooth( Total probability omitted when truncating a kernel with unbounded support. max_grid_points: - Maximum permitted size of the intermediate uniform grid. Raises an - error rather than silently reducing resolution if this limit is exceeded. + Intermediate uniform grids no larger than this are always allowed. Larger + grids may be rejected to guard against excessive resampling. Returns ------- @@ -576,7 +547,7 @@ def smooth( object, or ``max_grid_points`` is not an integer. """ x, y, template = _scipp_input(x, y) - values = _smooth_kernel_values( + values = _smooth_values( x.values, y.values, scale=_scale_in_coordinate_unit(scale, x), @@ -631,8 +602,8 @@ def smooth_relative( Total probability omitted when truncating a kernel with unbounded support or support reaching the nonpositive coordinate domain. max_grid_points: - Maximum permitted size of the intermediate geometric grid. Raises an - error rather than silently reducing resolution if this limit is exceeded. + Intermediate geometric grids no larger than this are always allowed. Larger + grids may be rejected to guard against excessive resampling. Returns ------- @@ -661,7 +632,7 @@ def smooth_relative( a distribution-like object, or ``max_grid_points`` is not an integer. """ x, y, template = _scipp_input(x, y) - values = _smooth_relative_kernel( + values = _smooth_relative_values( x.values, y.values, scale=_dimensionless_scale(scale), diff --git a/tests/smoothing_test.py b/tests/smoothing_test.py index e37004bae..2a73a4ee3 100644 --- a/tests/smoothing_test.py +++ b/tests/smoothing_test.py @@ -8,12 +8,7 @@ from scipy.stats import uniform import scippneutron as scn -from scippneutron.smoothing import ( - _relative_kernel_weights, - _smooth_relative_kernel_on_geomgrid, - smooth, - smooth_relative, -) +from scippneutron.smoothing import _relative_kernel_weights, smooth, smooth_relative def _normal_pdf(z): @@ -304,7 +299,7 @@ def test_fixed_width_gaussian_matches_exact_quadratic_on_nonuniform_grid(): def test_fixed_width_smoothing_rejects_resampled_grid_larger_than_limit(): - x, y = _variables([0.0, 0.25, 1.0], [1.0, 2.0, 3.0]) + x, y = _variables([0.0, 0.125, 1.0], [1.0, 2.0, 3.0]) with pytest.raises( ValueError, @@ -326,6 +321,40 @@ def test_fixed_width_smoothing_accepts_resampled_grid_equal_to_limit(): np.testing.assert_allclose(actual.values, y.values) +@pytest.mark.parametrize( + ("function", "x", "kwargs"), + [ + ( + smooth, + np.linspace(0.0, 1.0, 11), + {"scale": sc.scalar(0.1, unit='m')}, + ), + (smooth_relative, np.geomspace(1.0, 2.0, 11), {"scale": 0.1}), + ], +) +def test_input_grid_larger_than_limit_is_accepted(function, x, kwargs): + y = np.ones_like(x) + + actual = _smooth_values(function, x, y, max_grid_points=10, **kwargs) + + np.testing.assert_allclose(actual, y) + + +@pytest.mark.parametrize( + ("function", "x", "kwargs"), + [ + (smooth, np.array([0.0, 0.26, 1.0]), {"scale": sc.scalar(0.1, unit='m')}), + (smooth_relative, 2.0 ** np.array([0.0, 0.26, 1.0]), {"scale": 0.1}), + ], +) +def test_modest_resampling_larger_than_limit_is_accepted(function, x, kwargs): + y = np.ones_like(x) + + actual = _smooth_values(function, x, y, max_grid_points=4, **kwargs) + + np.testing.assert_allclose(actual, y) + + def test_fixed_width_pathologically_close_coordinates_fail_before_allocation(): x, y = _variables( [1.0, np.nextafter(1.0, 2.0), 2.0], @@ -568,9 +597,16 @@ def test_rejects_data_array_with_variances(): smooth_relative(data, scale=0.1) -def test_rejects_numpy_arrays(): +@pytest.mark.parametrize( + ("x", "y"), + [ + (np.arange(1.0, 4.0), sc.ones(dims=['x'], shape=[3])), + (sc.arange('x', 1.0, 4.0), np.ones(3)), + ], +) +def test_rejects_non_scipp_pair_member(x, y): with pytest.raises(TypeError, match="DataArray or a pair of Variables"): - smooth_relative(np.arange(1.0, 4.0), np.ones(3), scale=0.1) + smooth_relative(x, y, scale=0.1) def test_rejects_variables_with_different_dimensions(): @@ -653,8 +689,10 @@ def test_geometric_input_does_not_gain_a_point_from_roundoff(): assert actual.shape == y.shape -def test_pathologically_close_coordinates_fail_before_allocation(): - x = np.array([1.0, np.nextafter(1.0, 2.0), 2.0]) +@pytest.mark.filterwarnings("ignore:divide by zero encountered in scalar divide") +@pytest.mark.parametrize("x0", [1.0, 3.0, 10.0, 0.1]) +def test_pathologically_close_coordinates_fail_before_allocation(x0): + x = np.array([x0, np.nextafter(x0, 2.0 * x0), 2.0 * x0]) y = np.ones_like(x) with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"): @@ -693,14 +731,12 @@ def test_asymmetric_kernel_convolution_matches_direct_weighted_sum(): alpha = 0.2 kernel = uniform(loc=0.5, scale=1.0) y = np.arange(size, dtype=float) ** 2 + x = np.geomspace(1.0, np.exp((size - 1) * log_spacing), size) - offsets, weights = _relative_kernel_weights( - log_spacing=log_spacing, - scale=alpha, - kernel=kernel, - tail=1e-12, - max_offset=size - 1, - ) + offsets = np.arange(-size + 1, size) + lower = np.expm1((offsets - 0.5) * log_spacing) / alpha + upper = np.expm1((offsets + 0.5) * log_spacing) / alpha + weights = kernel.cdf(upper) - kernel.cdf(lower) expected = np.empty_like(y) for i in range(size): valid = (0 <= i + offsets) & (i + offsets < size) @@ -711,9 +747,10 @@ def test_asymmetric_kernel_convolution_matches_direct_weighted_sum(): else np.nan ) - actual = _smooth_relative_kernel_on_geomgrid( + actual = _smooth_values( + smooth_relative, + x, y, - log_spacing=log_spacing, scale=alpha, kernel=kernel, tail=1e-12, @@ -723,9 +760,12 @@ def test_asymmetric_kernel_convolution_matches_direct_weighted_sum(): def test_kernel_with_no_reachable_mass_returns_nan(): - actual = _smooth_relative_kernel_on_geomgrid( - np.arange(5.0), - log_spacing=0.1, + y = np.arange(5.0) + x = np.geomspace(1.0, np.exp((y.size - 1) * 0.1), y.size) + actual = _smooth_values( + smooth_relative, + x, + y, scale=0.1, kernel=uniform(loc=100.0, scale=1.0), tail=1e-12, From d9110a7015d46e6c09b9d53389f221079515e054 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Fri, 14 Aug 2026 06:28:46 +0000 Subject: [PATCH 6/6] refactor: unify the two smoothing paths behind a grid geometry The relative and translation-invariant procedures are the same algorithm under a coordinate map. Substituting u(z) = scale*z, z(u) = u/scale and z_min = -inf into the relative version reduces it exactly to the translation-invariant one: the reachable-mass normalization collapses to one, log1p/expm1 collapse to linear, and the branch selecting between exact support and a tail cutoff collapses to a plain isfinite check on the bounds in u. Both weight computations and both resample-smooth-interpolate drivers therefore become one of each, parameterized on a geometry supplying the working coordinate and the map between a displacement and its offset in that coordinate. This removes a class of bug rather than a number of lines. The grid size overflowing to infinity for coordinates that collide in the working coordinate was fixed on the relative path only; the uniform path still raised OverflowError for inputs whose spacing ratio overflows, such as [0.0, 5e-324, 1.0]. With a single path the two cannot diverge again. Move the stencil-clamping test onto the public API so it keeps testing the behavior rather than the helper that happens to implement it. --- src/scippneutron/smoothing.py | 302 ++++++++++++++++------------------ tests/smoothing_test.py | 44 ++--- 2 files changed, 165 insertions(+), 181 deletions(-) diff --git a/src/scippneutron/smoothing.py b/src/scippneutron/smoothing.py index 55f5887c2..ad1440ed6 100644 --- a/src/scippneutron/smoothing.py +++ b/src/scippneutron/smoothing.py @@ -107,117 +107,148 @@ def _trim_kernel_weights( return offsets[first:last], weights / weights.sum() -def _relative_kernel_weights( - log_spacing: float, - scale: float, - kernel: _Kernel, - tail: float, - max_offset: int, -) -> tuple[_IntArray, _FloatArray]: +class _Geometry(Protocol): + """ + A grid on which the kernel is translation invariant. + + Both smoothing procedures displace a coordinate x by a distribution Z. The + displacement law differs, but in both cases there is a working coordinate + ``u`` in which the displacement is an additive offset independent of x, so + that a single set of weights applies at every grid point. ``offset`` and + ``displacement`` convert between a displacement z and its offset in u; they + are inverses of each other. """ - Weights for smoothing on a geometric grid q_i = q0 * exp(i * log_spacing). - The kernel distribution describes the relative displacement Z: + #: Names the grid in user-facing messages. + name: str - q' = q * (1 + scale * Z) + def check(self, x: _FloatArray) -> None: + """Reject coordinates outside the domain of the displacement law.""" - Equivalently, + def coordinate(self, x: _FloatArray) -> _FloatArray: + """Working coordinate u(x).""" - K(q, q') = 1 / (scale * q) * f((q' - q) / (scale * q)) + def points(self, start: float, stop: float, count: int) -> _FloatArray: + """``count`` samples from ``start`` to ``stop``, uniform in u.""" - where f is the PDF of the supplied distribution. - """ - dist = _as_kernel_distribution(kernel) + def displacement_min(self, scale: float) -> float: + """Smallest displacement z that keeps x within the valid domain.""" - # Positive physical domain: - # - # q' > 0 - # q * (1 + scale * z) > 0 - # z > -1 / scale - z_domain_min = -1.0 / scale + def offset(self, scale: float, z: ArrayLike) -> Any: + """Offset in u produced by displacement z.""" - p_domain_min = float(dist.cdf(z_domain_min)) - norm_mass = 1.0 - p_domain_min + def displacement(self, scale: float, u: ArrayLike) -> Any: + """Displacement z producing an offset u.""" - if not np.isfinite(norm_mass) or norm_mass <= 0.0: - raise ValueError("kernel has no positive-domain mass for this scale") - support_min, support_max = _kernel_support(dist) +class _Uniform: + """Constant kernel width: ``x' = x + scale * Z``, additive in x itself.""" - # Exact z-support after clipping to q' > 0. - z_left_exact = max(z_domain_min, support_min) - z_right_exact = support_max + name = "uniform" - # If the clipped support maps to finite log-space, use it exactly. - # If it touches q' = 0, the log lower bound is -inf, so use a tail cutoff. - has_finite_log_support = ( - np.isfinite(z_left_exact) - and np.isfinite(z_right_exact) - and (1.0 + scale * z_left_exact > 0.0) - ) + def check(self, x: _FloatArray) -> None: + pass - if has_finite_log_support: - u_left = np.log1p(scale * z_left_exact) - u_right = np.log1p(scale * z_right_exact) - else: - probabilities = ( - p_domain_min + np.array([0.5 * tail, 1.0 - 0.5 * tail]) * norm_mass - ) - z_left, z_right = (float(value) for value in dist.ppf(probabilities)) + def coordinate(self, x: _FloatArray) -> _FloatArray: + return x - u_left = np.log1p(scale * z_left) - u_right = np.log1p(scale * z_right) + def points(self, start: float, stop: float, count: int) -> _FloatArray: + return np.linspace(start, stop, count) - if not np.isfinite(u_right): - raise ValueError("right kernel bound is not finite; increase tail") + def displacement_min(self, scale: float) -> float: + return -np.inf - # Cells are centered at m*h and span [(m-1/2)h, (m+1/2)h]. - # Include offset zero even for one-sided kernels and clamp the stencil to - # offsets that can contribute to the finite input. - m_min = int(np.clip(np.floor(u_left / log_spacing + 0.5), -max_offset, 0)) - m_max = int(np.clip(np.ceil(u_right / log_spacing - 0.5), 0, max_offset)) + def offset(self, scale: float, z: ArrayLike) -> Any: + return scale * np.asarray(z, dtype=float) - m = np.arange(m_min, m_max + 1, dtype=np.int64) + def displacement(self, scale: float, u: ArrayLike) -> Any: + return np.asarray(u, dtype=float) / scale + + +class _Geometric: + """ + Relative kernel width: ``x' = x * (1 + scale * Z)``, additive in ``log(x)``. + + Equivalently, the kernel is + + K(x, x') = 1 / (scale * x) * f((x' - x) / (scale * x)) - L = (m - 0.5) * log_spacing - U = (m + 0.5) * log_spacing + for a distribution with PDF f. Displacements are restricted to ``x' > 0``, + that is ``1 + scale * z > 0``. + """ + + name = "geometric" + + def check(self, x: _FloatArray) -> None: + if np.any(x <= 0): + raise ValueError("x must be positive") - # z = (q' - q) / (scale q) - # = (exp(u) - 1) / scale - zL = np.expm1(L) / scale - zU = np.expm1(U) / scale + def coordinate(self, x: _FloatArray) -> _FloatArray: + return cast(_FloatArray, np.log(x)) - # Exact cell-integrated weights in log-space. - w = (dist.cdf(zU) - dist.cdf(zL)) / norm_mass - w = np.maximum(w, 0.0) + def points(self, start: float, stop: float, count: int) -> _FloatArray: + return np.geomspace(start, stop, count) - return _trim_kernel_weights(m, w) + def displacement_min(self, scale: float) -> float: + return -1.0 / scale + def offset(self, scale: float, z: ArrayLike) -> Any: + # log1p(-1) is -inf, which _kernel_weights treats as unbounded support. + with np.errstate(divide="ignore", invalid="ignore"): + return np.log1p(scale * np.asarray(z, dtype=float)) -def _translation_invariant_kernel_weights( + def displacement(self, scale: float, u: ArrayLike) -> Any: + return np.expm1(np.asarray(u, dtype=float)) / scale + + +def _kernel_weights( + geometry: _Geometry, spacing: float, scale: float, kernel: _Kernel, tail: float, max_offset: int, ) -> tuple[_IntArray, _FloatArray]: + """ + Cell-integrated kernel weights on a grid of the given spacing in u. + + Weights are exact integrals of the kernel probability over grid cells, so + the result is a proper quadrature of the smoothing integral rather than a + point sampling of the kernel. + """ dist = _as_kernel_distribution(kernel) - z_left, z_right = _kernel_support(dist) - if not np.all(np.isfinite([z_left, z_right])): - z_left, z_right = ( - float(value) for value in dist.ppf([0.5 * tail, 1.0 - 0.5 * tail]) - ) - bounds = scale * np.array([z_left, z_right]) - if not np.all(np.isfinite(bounds)): - raise ValueError("kernel bounds are not finite; increase tail") + z_min = geometry.displacement_min(scale) + p_min = float(dist.cdf(z_min)) + reachable_mass = 1.0 - p_min + if not np.isfinite(reachable_mass) or reachable_mass <= 0.0: + raise ValueError("kernel has no mass in the valid domain for this scale") - m_min = int(np.clip(np.floor(bounds[0] / spacing + 0.5), -max_offset, 0)) - m_max = int(np.clip(np.ceil(bounds[1] / spacing - 0.5), 0, max_offset)) + support_min, support_max = _kernel_support(dist) + u_left = float(geometry.offset(scale, max(z_min, support_min))) + u_right = float(geometry.offset(scale, support_max)) + + # Unbounded support, or support reaching the edge of the valid domain, + # gives an infinite bound in u. Truncate at the requested tail instead. + if not (np.isfinite(u_left) and np.isfinite(u_right)): + probabilities = ( + p_min + np.array([0.5 * tail, 1.0 - 0.5 * tail]) * reachable_mass + ) + u_left, u_right = ( + float(value) for value in geometry.offset(scale, dist.ppf(probabilities)) + ) + if not np.isfinite(u_right): + raise ValueError("right kernel bound is not finite; increase tail") + + # Cells are centered at m*h and span [(m-1/2)h, (m+1/2)h]. + # Include offset zero even for one-sided kernels and clamp the stencil to + # offsets that can contribute to the finite input. + m_min = int(np.clip(np.floor(u_left / spacing + 0.5), -max_offset, 0)) + m_max = int(np.clip(np.ceil(u_right / spacing - 0.5), 0, max_offset)) m = np.arange(m_min, m_max + 1, dtype=np.int64) - lower = (m - 0.5) * spacing / scale - upper = (m + 0.5) * spacing / scale + lower = geometry.displacement(scale, (m - 0.5) * spacing) + upper = geometry.displacement(scale, (m + 0.5) * spacing) weights = np.maximum(dist.cdf(upper) - dist.cdf(lower), 0.0) return _trim_kernel_weights(m, weights) @@ -280,13 +311,14 @@ def _validate_max_grid_points(max_grid_points: int) -> None: raise ValueError("max_grid_points must be at least 2") -def _smooth_relative_values( +def _smooth_values( + geometry: _Geometry, x: ArrayLike, y: ArrayLike, scale: float, - kernel: _Kernel = "gaussian", - tail: float = 1e-12, - max_grid_points: int = 1_000_000, + kernel: _Kernel, + tail: float, + max_grid_points: int, ) -> _FloatArray: x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) @@ -302,112 +334,56 @@ def _smooth_relative_values( _validate_max_grid_points(max_grid_points) if np.any(~np.isfinite(x)): raise ValueError("x must contain only finite values") - if np.any(x <= 0): - raise ValueError("x must be positive") + geometry.check(x) if np.any(np.diff(x) <= 0): raise ValueError("x must be strictly increasing") - if x.size == 0 or scale == 0 or x.size == 1: + if x.size < 2 or scale == 0: return y.copy() - logx = np.log(x) - dlog = np.diff(logx) - log_range = float(logx[-1] - logx[0]) + u = geometry.coordinate(x) + du = np.diff(u) + u_range = float(u[-1] - u[0]) - # Preserve an existing geometric grid. Otherwise choose a geometric grid - # at least as dense as the smallest input spacing in log-space. - if np.allclose(dlog, dlog[0], rtol=1e-7, atol=0.0): - k = x.size + # Preserve an existing regular grid. Otherwise choose a grid at least as + # dense as the smallest input spacing, measured in the working coordinate. + if np.allclose(du, du[0], rtol=1e-7, atol=0.0): + k = float(x.size) else: - k = np.ceil(log_range / np.min(dlog)) + 1.0 + # Coordinates that are distinct but collide, or nearly so, in the + # working coordinate give an infinite point count. Let it propagate to + # the guard below rather than raising here. + with np.errstate(divide="ignore", over="ignore"): + k = np.ceil(u_range / np.min(du)) + 1.0 if not (k <= max_grid_points) and k > 2 * x.size: raise ValueError( - "geometric resampling would require too many points, exceeding " + f"{geometry.name} resampling would require too many points, exceeding " f"max_grid_points={max_grid_points:,}. Increase max_grid_points to " "allow a larger grid." ) k = int(k) - log_spacing = log_range / (k - 1) - if not np.isfinite(log_spacing) or log_spacing <= 0: - raise ValueError("log_spacing must be positive") - - xp = np.geomspace(x[0], x[-1], k) - yg = np.interp(xp, x, y) - offsets, weights = _relative_kernel_weights( - scale=scale, - log_spacing=log_spacing, - kernel=kernel, - tail=tail, - max_offset=k - 1, - ) - zg = _smooth_with_weights(yg, offsets, weights) - return cast(_FloatArray, np.interp(x, xp, zg)) - - -def _smooth_values( - x: ArrayLike, - y: ArrayLike, - scale: float, - kernel: _Kernel = "gaussian", - tail: float = 1e-12, - max_grid_points: int = 1_000_000, -) -> _FloatArray: - x = np.asarray(x, dtype=float) - y = np.asarray(y, dtype=float) - - if x.ndim != 1 or y.ndim != 1: - raise ValueError("x and y must be one-dimensional") - if x.size != y.size: - raise ValueError("x and y must have the same length") - if not np.isfinite(scale) or scale < 0: - raise ValueError("scale must be non-negative") - if not np.isfinite(tail) or not (0.0 < tail < 1.0): - raise ValueError("tail must be between 0 and 1") - _validate_max_grid_points(max_grid_points) - if np.any(~np.isfinite(x)): - raise ValueError("x must contain only finite values") - if np.any(np.diff(x) <= 0): - raise ValueError("x must be strictly increasing") - if x.size == 0 or scale == 0 or x.size == 1: - return y.copy() - spacing = np.diff(x) - x_range = float(x[-1] - x[0]) - - # Preserve an existing uniform grid. Otherwise choose a uniform grid at - # least as dense as the smallest input spacing. - if np.allclose(spacing, spacing[0], rtol=1e-7, atol=0.0): - k = x.size - else: - k = int(np.ceil(x_range / np.min(spacing))) + 1 - - if k > max_grid_points and k > 2 * x.size: - raise ValueError( - "uniform resampling would require too many points, exceeding " - f"max_grid_points={max_grid_points:,}. Increase max_grid_points to " - "allow a larger grid." - ) - - spacing = x_range / (k - 1) + spacing = u_range / (k - 1) if not np.isfinite(spacing) or spacing <= 0: - raise ValueError("spacing must be positive") - offsets, weights = _translation_invariant_kernel_weights( + raise ValueError("grid spacing must be positive") + + offsets, weights = _kernel_weights( + geometry, spacing=spacing, scale=scale, kernel=kernel, tail=tail, max_offset=k - 1, ) - xp = np.linspace(x[0], x[-1], k) - yg = np.interp(xp, x, y) - - zg = _smooth_with_weights(yg, offsets, weights) + xp = geometry.points(float(x[0]), float(x[-1]), k) + zg = _smooth_with_weights(np.interp(xp, x, y), offsets, weights) return cast(_FloatArray, np.interp(x, xp, zg)) def _scipp_input( - x: _ScippArray, y: sc.Variable | None + x: object, y: object ) -> tuple[sc.Variable, sc.Variable, sc.DataArray | None]: + """Validate untyped user input and reduce it to a coordinate/data pair.""" template: sc.DataArray | None = None if isinstance(x, sc.DataArray): if y is not None: @@ -548,6 +524,7 @@ def smooth( """ x, y, template = _scipp_input(x, y) values = _smooth_values( + _Uniform(), x.values, y.values, scale=_scale_in_coordinate_unit(scale, x), @@ -632,7 +609,8 @@ def smooth_relative( a distribution-like object, or ``max_grid_points`` is not an integer. """ x, y, template = _scipp_input(x, y) - values = _smooth_relative_values( + values = _smooth_values( + _Geometric(), x.values, y.values, scale=_dimensionless_scale(scale), diff --git a/tests/smoothing_test.py b/tests/smoothing_test.py index 2a73a4ee3..0c5f8f5c0 100644 --- a/tests/smoothing_test.py +++ b/tests/smoothing_test.py @@ -8,7 +8,7 @@ from scipy.stats import uniform import scippneutron as scn -from scippneutron.smoothing import _relative_kernel_weights, smooth, smooth_relative +from scippneutron.smoothing import smooth, smooth_relative def _normal_pdf(z): @@ -355,11 +355,19 @@ def test_modest_resampling_larger_than_limit_is_accepted(function, x, kwargs): np.testing.assert_allclose(actual, y) -def test_fixed_width_pathologically_close_coordinates_fail_before_allocation(): - x, y = _variables( +@pytest.mark.parametrize( + "coordinates", + [ [1.0, np.nextafter(1.0, 2.0), 2.0], - [1.0, 1.0, 1.0], - ) + # The spacing ratio overflows to infinity rather than merely being large. + [0.0, 5e-324, 1.0], + [1.0, np.nextafter(1.0, 2.0), 1e300], + ], +) +def test_fixed_width_pathologically_close_coordinates_fail_before_allocation( + coordinates, +): + x, y = _variables(coordinates, [1.0, 1.0, 1.0]) with pytest.raises(ValueError, match="exceeding max_grid_points=1,000,000"): smooth(x, y, scale=sc.scalar(0.1, unit='m')) @@ -709,20 +717,18 @@ def test_wide_coordinate_range_does_not_overflow_grid_construction(): def test_kernel_stencil_is_bounded_before_allocation(): - max_offset = 1000 - - offsets, weights = _relative_kernel_weights( - log_spacing=1e-6, - scale=1.0, - kernel="gaussian", - tail=1e-12, - max_offset=max_offset, - ) - - assert offsets[0] >= -max_offset - assert offsets[-1] <= max_offset - assert offsets.size <= 2 * max_offset + 1 - assert offsets.size == weights.size + # A gaussian truncated at tail=1e-12 spans roughly 14 sigma, which is about + # 1e7 offsets on this grid. The stencil must be clamped to offsets that can + # reach the input rather than allocated at its nominal width. + size = 21 + x = np.geomspace(1.0, np.exp((size - 1) * 1e-6), size) + y = np.arange(size, dtype=float) + + actual = _smooth_values(smooth_relative, x, y, scale=1.0) + + # The kernel is flat to within 1e-5 over the whole input, so every point + # averages the entire array. + np.testing.assert_allclose(actual, np.full(size, y.mean()), rtol=1e-4) def test_asymmetric_kernel_convolution_matches_direct_weighted_sum():