From 35ab06ce08b2e6178f5db900c3ffe45afc501c69 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Thu, 17 Jul 2025 15:32:04 -0700 Subject: [PATCH] making convolutional smoother more efficient to address #129 --- pynumdiff/tests/test_diff_methods.py | 8 ++++---- pynumdiff/tests/test_utils.py | 10 ++++++++++ pynumdiff/utils/utility.py | 18 ++++++++---------- 3 files changed, 22 insertions(+), 14 deletions(-) diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index d0bc4b7..731fd00 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -33,7 +33,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) (iterated_first_order, {'num_iterations':2}), (iterated_first_order, [2], {'iterate':True}), (lineardiff, {'order':3, 'gamma':5, 'window_size':11, 'solver':'CLARABEL'}), (lineardiff, [3, 5, 11], {'solver':'CLARABEL'}), (polydiff, {'poly_order':2, 'window_size':3}), (polydiff, [2, 3]), - (savgoldiff, {'poly_order':2, 'window_size':4, 'smoothing_win':4}), (savgoldiff, [2, 4, 4]), + (savgoldiff, {'poly_order':2, 'window_size':5, 'smoothing_win':5}), (savgoldiff, [2, 5, 5]), (spectraldiff, {'high_freq_cutoff':0.1}), (spectraldiff, [0.1]), (mediandiff, {'window_size':3, 'num_iterations':2}), (mediandiff, [3, 2], {'iterate':True}), (meandiff, {'window_size':3, 'num_iterations':2}), (meandiff, [3, 2], {'iterate':True}), @@ -95,13 +95,13 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) [(0, 0), (3, 3), (0, 0), (3, 3)]], savgoldiff: [[(-9, -10), (-13, -14), (0, -1), (0, 0)], [(-9, -10), (-13, -13), (0, -1), (0, 0)], - [(-1, -1), (0, -1), (0, -1), (0, 0)], - [(0, -1), (0, 0), (0, -1), (1, 0)], + [(-2, -2), (-1, -1), (0, -1), (0, 0)], + [(0, -1), (0, 0), (0, 0), (1, 0)], [(1, 1), (2, 2), (1, 1), (2, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]], spectraldiff: [[(-9, -10), (-14, -15), (-1, -1), (0, 0)], [(0, 0), (1, 1), (0, 0), (1, 1)], - [(1, 0), (1, 1), (1, 1), (1, 1)], + [(1, 1), (1, 1), (1, 1), (1, 1)], [(0, 0), (1, 1), (0, 0), (1, 1)], [(1, 1), (2, 2), (1, 1), (2, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]], diff --git a/pynumdiff/tests/test_utils.py b/pynumdiff/tests/test_utils.py index 59ee59b..51f4ae3 100644 --- a/pynumdiff/tests/test_utils.py +++ b/pynumdiff/tests/test_utils.py @@ -33,6 +33,16 @@ def test_estimate_integration_constant(): assert 0.9 < x0 < 1.1 # The result should be close to 1.0, but not exactly due to noise +def test_convolutional_smoother(): + """Ensure the convolutional smoother isn't introducing edge effects""" + x = np.ones(10) + kernel_odd = np.ones(3)/3 + kernel_even = np.ones(4)/4 + + assert np.allclose(utility.convolutional_smoother(x, kernel_odd, num_iterations=3), np.ones(len(x))) + assert np.allclose(utility.convolutional_smoother(x, kernel_even, num_iterations=3), np.ones(len(x))) + + def test_hankel_matrix(): """Ensure Hankel matrix comes back as defined""" assert np.allclose(utility.hankel_matrix([1, 2, 3, 4, 5], 3), [[1, 2, 3],[2, 3, 4],[3, 4, 5]]) diff --git a/pynumdiff/utils/utility.py b/pynumdiff/utils/utility.py index 108f18a..27ecdef 100644 --- a/pynumdiff/utils/utility.py +++ b/pynumdiff/utils/utility.py @@ -140,25 +140,23 @@ def friedrichs_kernel(window_size): ker = np.exp(-1/(1-x**2)) return ker / np.sum(ker) -def convolutional_smoother(x, kernel, iterations=1): +def convolutional_smoother(x, kernel, num_iterations=1): """Perform smoothing by convolving x with a kernel. :param np.array[float] x: 1D data :param np.array[float] kernel: kernel to use in convolution - :param int iterations: number of iterations, >=1 + :param int num_iterations: number of iterations, >=1 :return: **x_hat** (np.array[float]) -- smoothed x """ - x_hat = np.hstack((x[::-1], x, x[::-1])) # pad - w = np.linspace(0, 1, len(x_hat)) # weights + pad_width = len(kernel)//2 + x_hat = x - for _ in range(iterations): - x_hat_f = np.convolve(x_hat, kernel, 'same') - x_hat_b = np.convolve(x_hat[::-1], kernel, 'same')[::-1] - - x_hat = x_hat_f*w + x_hat_b*(1-w) + for i in range(num_iterations): + x_padded = np.pad(x_hat, pad_width, mode='symmetric') # pad with repetition of the edges + x_hat = np.convolve(x_padded, kernel, 'valid')[:len(x)] # 'valid' slices out only full-overlap spots - return x_hat[len(x):len(x)*2] + return x_hat def slide_function(func, x, dt, kernel, *args, stride=1, pass_weights=False, **kwargs):