Merge pull request #130 from florisvb/simplify-convolutional-smoother

making convolutional smoother more efficient to address #129

Author: pavelkomarov

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)]],

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]])

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):