feat: relative width smoothing procedures

[jokasimr]

Adds functionality to scippneutron for smoothing 1D curves with kernels that have a width relative to the dimension coordinate.

As an example, in the figure below the step signal is smoothed by a gaussian kernel with width proportional to x:

The computed smoothed signal and the expected exact signal overlap.

[jokasimr]

The idea is that the scaled kernel has constant width on a grid that is the logarithm of the original grid, and when the kernel is constant width it can be computed efficiently using scipy.signal.convolve.

Here's a document describing the implementation in more detail:

relative_kernel_smoothing_math.pdf

[jokasimr]

@SimonHeybrock Do you think this functionality should go here in Scippneutron or in Scipp?

It's mainly useful for analysis work where you have an idealized model and want to apply some smoothing from resolution effects to obtain something that is comparable to the measurement data.

[SimonHeybrock]

Before looking into details, can you compare this to what we already have in https://scipp.github.io/generated/modules/scipy/scipp.scipy.ndimage.gaussian_filter.html (and https://scipp.github.io/generated/modules/scipy/scipp.scipy.ndimage.generic_filter.html)? What is added on top of that? Would it make sense to have a common API?

[jokasimr]

The difference is that this is a more convenient interface for typical smoothing needs, more of a interface for users than for developers.

1. (Grid awareness) The functions implemented in this module are grid aware while gaussian_filter is not grid aware. That is, for gaussian_filter, the effect on the signal will depend on the grid it was defined on. If the intention of the user is to smooth the input with a fixed-width gaussian kernel then they have to make sure their signal is defined on a uniform grid, etc.

2. (Flexible kernel support) This module exposes a number of common smoothing kernels, and any distribution in scipy.stats can be used as a kernel.

3. (Relative width kernels) smooth_relative is a convenience function for smoothing with a kernel having a width that grows proportional to the grid coordinate. That is a common case, think about instruments where uncertainty in a coordinate is proportional to itself (q or wavelength uncertainty grows proportional to q or wavelength). But such operations can be quite hard to implement correctly and efficiently using a tool like generic_filter.

Since the functions implemented in this module use scipy.signal.convolve under the hood they are probably significantly faster than generic_filter in most cases. But I have not made any runtime comparisons.

The functions in this module are explicitly 1D, and not ND like the functions you mentioned. That is a limitation, but it allows a simpler interface and it covers common use cases.

[SimonHeybrock]

Nice work — the numerical core is careful and the test suite is the strongest part of this PR. Analytic reference solutions for the smoothed quadratic, second-order convergence checks against grid spacing, and tail-truncation scaling with a grid-limited plateau is real verification rather than smoke tests. The convolution alignment, the O(n) boundary normalization in _valid_weight_sums, and the direct-vs-FFT switch for non-finite input all check out.

Findings below, roughly in order of importance. One is a genuine crash (OverflowError on the relative path for coordinates that collide in log space), two are smaller error-handling issues, and the rest is structure and documentation.

The largest item is that the relative and translation-invariant paths are one algorithm under a coordinate map — I prototyped the unification to confirm that, details on _translation_invariant_kernel_weights.

Separately, and not a code review point so not raised inline: I still want to settle whether this belongs in scippneutron at all, given it is 672 lines with no neutron-specific content. Let's discuss that on the thread.

[SimonHeybrock]

This and _relative_kernel_weights are the same function under a coordinate map, and the reduction is exact rather than approximate: substitute u(z) = scale*z, z(u) = u/scale, z_domain_min = -inf into the relative version and norm_mass collapses to 1, has_finite_log_support collapses to the plain isfinite check, and log1p/expm1 collapse to linear. One level up, _smooth_relative_kernel and _smooth_kernel_values differ only by log-vs-identity plus the x > 0 check.

I prototyped the unified version to check this isn't hand-waving: the four functions become two, parameterized on a small geometry object supplying coordinate, points, displacement_min, offset and displacement. The core drops from 436 to ~290 lines, all 73 tests in this PR pass unchanged, and across 180 combinations of grid (uniform, geometric, jittered, 4-point, 1e-6..1e6), kernel (gaussian, boxcar, triangular, asymmetric uniform, asymmetric triangular, expon) and scale, the outputs agree with this branch to 1e-10.

Two things fall out for free. The truncation condition becomes isfinite(u_left) and isfinite(u_right), which is provably equivalent to has_finite_log_support but easier to check by eye. And computing k in float for the shared max_grid_points guard fixes the overflow reported on line 358 by construction.

Happy to push the prototype if useful.

[jokasimr]

I don't think this is worth DRYing up by unifying the two cases. It's better to keep them separate because I think that makes it easier to understand, it's already complex enough.

[SimonHeybrock]

But less duplication is easier to review.

[jokasimr]

I don't agree that's obviously the case.
I've already said what I think about it. Do you think I should do the refactor suggested by the model?

[SimonHeybrock]

yes

[jokasimr]

I tried a refactor to unify _relative_kernel_weights and _translation_invariant_kernel_weights, but it did not result in any real code reduction, while being more complex.

It's possible there's a better shape that actually reduces code size and complexity. You are welcome to push such a version to this branch if you prefer it.

[SimonHeybrock]

Two properties are worth stating explicitly here, because both public functions are normalized weighted averages rather than convolutions:

• A constant is preserved exactly, but a sum is not. Boundary renormalization pulls mass back into the domain, so a Gaussian peak's total changes by ~2% at scale=1.0 on a 200-point grid. That is the right behavior for the use case in the PR description (smearing an idealized model to compare against data), but it makes the function unsuitable for smearing raw counts — and the tests use unit='counts' throughout, which invites exactly that reading.
• The result may contain NaN where no kernel mass is reachable: at the boundary for one-sided kernels, and for the whole array when the kernel has no reachable mass at all (the nonzero.size == 0 path at line 88). This is deliberate and tested, but neither Returns section mentions it.

[jokasimr]

I think the first point has to be clarified a bit. Exactly what was the test case and what was the result of that test?

It's not correct to say that the normalization behavior "makes the function unsuitable for smearing raw counts". That is a far to general statement. It is just context dependent.

The smoothing operation needs to have a strategy for handling boundaries. Here the strategy used is to re-normalize by the mass of the kernel that falls inside the boundary. This can be understood as an assumption that the (kernel-)weighted mean of the signal outside the boundary is the same as the weighted mean of the signal inside the boundary.
In practice that is a conservative assumption, we don't expect anything drastic happening to the signal exactly at the boundary.
But of course that assumption will be wrong sometimes, any assumption is.

How to deal with that as a user

As a user we might know something about the behavior of the signal outside the domain, for example, we might know it is zero, or we know it decays following a certain patterns, or something else.
In almost all such cases the smoothed signal near the boundary will not be what it would be if we had taken our extra knowledge into account properly.

If the user wants perfect boundary behavior the best option for them is to extend the signal that they pass to the smoother with the "tails" that they assume it has outside of the bounds of the signal. Then they can crop out the portion of the smoothed signal that overlaps with the real signal, and that entire portion will have been unaffected by "boundary effects".

Another (simpler) option is to smooth the signal and then cut out a center section of the smoothed signal that is unaffected by boundary effects.

[SimonHeybrock]

No module docstring. The geometric-grid trick — that a relative-width kernel is translation invariant in log(x), which is what makes scipy.signal.convolve applicable — is the central idea here, and it currently lives only in a PDF attached to a PR comment. That should be in the module docstring so it survives.

[jokasimr]

Makes sense.

[SimonHeybrock]

Pushed d9110a70. Fair warning that the line count barely moves (643 → 619, 190 → 177 statements) — you were right that unification is not primarily a size win. The case for it turned out to be elsewhere.

The reduction is exact rather than approximate, which is what makes the shape work: substituting u(z) = scale*z, z(u) = u/scale and z_min = -inf into _relative_kernel_weights gives _translation_invariant_kernel_weights identically. norm_mass collapses to 1, log1p/expm1 collapse to linear, and has_finite_log_support collapses to a plain isfinite check on the bounds in u. So the geometry only has to supply the working coordinate and the map between a displacement and its offset in that coordinate — five small methods, no configuration and no branching inside the shared code.

What convinced me it was worth pushing rather than dropping: the paths had already diverged again. 7a820602 fixed the infinite grid size on the relative path, but the uniform path still raised OverflowError for inputs whose spacing ratio overflows rather than merely being large:

smooth(sc.array(dims=['x'], values=np.array([0.0, 5e-324, 1.0]), unit='m'), y, scale=sc.scalar(0.1, unit='m'))
# OverflowError: cannot convert float infinity to integer

Same for [1.0, nextafter(1.0), 1e300]. Both now raise the intended ValueError, and with one path they cannot drift apart again. Added those as cases on the existing test.

Verification, since this is a numerics change and the diff is not small:

• All existing tests pass unchanged apart from the two noted below.
• 896 combinations of grid (uniform, geometric, jittered, 4-point, 2-point, 1e-6..1e6, 1..1e12, clustered) × kernel (gaussian, boxcar, triangular, asymmetric uniform, asymmetric triangular, expon, shifted normal) × scale × tail × both functions, compared against 7a820602: identical values to 1e-11 and identical exception types, with the OverflowError → ValueError cases above the only intended differences.
• mypy on this file goes from three errors to none. Two were the spacing name being reused for both an array and a float; the third was the isinstance(x, sc.Variable) guard being statically unreachable, which I fixed by typing _scipp_input as taking object, since it is the boundary where untyped user input is validated.

Two test changes, both flagged in the review as testing implementation rather than behavior. test_kernel_stencil_is_bounded_before_allocation now goes through smooth_relative on a grid where the gaussian is flat across the whole input, so it asserts every point averages the entire array — that fails just as loudly if the stencil clamp is removed, without naming a private helper. That removes the last private import from the test module.

Untouched, since they are yours to decide rather than mechanical: the float32 → float64 widening, and documenting that these are normalized weighted averages, so a constant is preserved exactly but a sum is not (~2% on a peak at scale=1.0), and that the result can legitimately contain NaN.

[SimonHeybrock]

@jokasimr What do you think about the pushed change? I like the _Geometry protocol, since it aids understanding.

[jokasimr]

I think testing the internal helper was simpler and clearer here. It needed less explanatory comments. I'm also not sure if this tests the same things.

Doesn't seem worth it to me just to get the private helper out of the tests.

[jokasimr]

@jokasimr What do you think about the pushed change? I like the _Geometry protocol, since it aids understanding.

I think it looks good 👍

I'm just a little bit worried about the changed test, it would be good to either re-add that, or test it locally manually once to make sure there's no behavioral change.

[SimonHeybrock]

This review was created with the help of AI. If anything below strikes you as unpleasantly formulated (tone, verbosity, ...), as too nit-picky, or as otherwise improvable, please tell me -- I am actively trying to improve this.

The algorithm looks correct to me. I checked the two convolution helpers independently against a brute-force double loop over 300 random asymmetric stencils, and both the alignment (start = offsets[-1]) and the O(n) boundary normalisation come out exact.

Nice work on the tests! Comparing against an analytically smoothed quadratic via truncated-normal moments, rather than golden values, means they would catch a wrong answer and not merely a changed one; test_asymmetric_kernel_convolution_matches_direct_weighted_sum brute-forces the convolution, which is exactly the test I would have asked for; and the nextafter / 5e-324 / 1e300 cases are more thorough than most numerical code gets.

Comments below are mostly small. The one I would like settled before merge is the input API (see the comment on _scipp_input), since it is public and awkward to walk back later.

Aside, not for this PR: docs/user-guide/algorithms-background/ would be a natural home for the figure in the PR description at some point -- it currently only lives here, and it illustrates the behaviour well.

[SimonHeybrock]

Both entry points accept either a DataArray or an (x, y) pair of Variables. I would suggest dropping the pair form and taking a DataArray only, matching scipp.scipy: ndimage.gaussian_filter(x, /, *, sigma=...) and signal.sosfiltfilt(obj, dim, *, sos=...) both take a single data argument rather than a coord/data pair.

The pair form saves the caller a sc.DataArray(y, coords={'x': x}). In exchange, this function re-validates dimensionality, matching dims and binned-ness that a DataArray enforces structurally; y becomes a positional-optional with a "must be omitted" TypeError; and the first parameter is named x even when it holds the data.

[SimonHeybrock]

_kernel_support has a helpful message for unfrozen distributions:

kernel must be a fully specified distribution. For distributions with shape parameters, pass a frozen distribution such as triang(c=0.5), not triang.

It is unreachable. dist.cdf(z_min) above runs first, so kernel=triang gives the user TypeError: _parse_args() missing 1 required positional argument: 'c' instead. Distributions without shape parameters (norm, uniform) pass the CDF call and never reach the branch either, so it never fires at all.

Taking the support first fixes it:

Suggested change

	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")
	support_min, support_max = _kernel_support(dist)
	z_min = geometry.displacement_min(scale)
	support_min, support_max = _kernel_support(dist)
	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")

[SimonHeybrock]

A binned DataArray hits x.coords.is_edges(x.dim) first and gets CoordError: the dimension coordinate must not contain bin edges, which points at the wrong thing.

This check is only reachable through the Variable-pair path, so if that path goes away it becomes dead code unless it moves ahead of the is_edges check.

[SimonHeybrock]

This part was implemented by myself, looking at it again now I see that scale is threaded through displacement_min, offset and displacement, and is the same value for the whole call. Binding it at construction (_Geometric(scale)) would drop the parameter from three signatures, take _kernel_weights from six parameters to five, and make a mismatched scale between the offset/displacement inverse pair impossible to express.

[SimonHeybrock]

Worth a comment here: there is no matching finiteness check on u_left, which instead clamps to -max_offset via np.clip below. That is right on both sides -- the left bound is bounded by the valid domain, the right one is not, and the geometric case with scale=1 relies on the clamp -- but the asymmetry reads as an oversight without a line saying why.

Minor, same block: when only one side of the support is infinite, both bounds get recomputed from ppf, so tail/2 is discarded off a side that was already finite. Within budget, just slightly wasteful.

[SimonHeybrock]

Suggested change

	def test_gaussian_smoothing_matches_exact_quadratic_on_geometric_grid():
	def test_relative_width_gaussian_matches_exact_quadratic_on_geometric_grid():

[SimonHeybrock]

Suggested change

	def test_gaussian_smoothing_error_is_second_order_in_log_grid_spacing():
	def test_relative_width_gaussian_error_is_second_order_in_log_grid_spacing():

[SimonHeybrock]

Suggested change

	def test_gaussian_smoothing_error_scales_with_inverse_kernel_width():
	def test_relative_width_gaussian_error_scales_with_inverse_kernel_width():

[SimonHeybrock]

Suggested change

	def test_gaussian_smoothing_error_scales_with_tail_until_grid_error_dominates():
	def test_relative_width_gaussian_error_scales_with_tail_until_grid_error_dominates():

[SimonHeybrock]

Suggested change

	def test_compact_symmetric_kernel_matches_exact_interior_quadratic_moments(
	kernel, relative_variance, max_error
	):
	def test_relative_width_compact_kernel_matches_interior_quadratic_moments(
	kernel, relative_variance, max_error
	):