diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index 4299555..9cd7f0e 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -137,7 +137,7 @@ def spline_irreg_step(*args, **kwargs): return splinediff(*args, **kwargs) [(0, 0), (1, 1), (0, 0), (1, 1)], [(1, 0), (2, 2), (1, 0), (2, 2)], [(1, 0), (3, 3), (1, 0), (3, 3)]], - polydiff: [[(-14, -15), (-14, -14), (0, -1), (1, 1)], + polydiff: [[(-14, -15), (-13, -14), (0, -1), (1, 1)], [(-14, -14), (-13, -13), (0, -1), (1, 1)], [(-14, -14), (-13, -13), (0, -1), (1, 1)], [(-2, -2), (0, 0), (0, -1), (1, 1)], @@ -179,11 +179,11 @@ def spline_irreg_step(*args, **kwargs): return splinediff(*args, **kwargs) [(0, 0), (1, 0), (0, -1), (1, 0)], [(1, 1), (2, 2), (1, 1), (2, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]], - jerk_sliding: [[(-15, -15), (-16, -16), (0, -1), (1, 0)], + jerk_sliding: [[(-25, -25), (-16, -16), (0, -1), (1, 0)], [(-14, -14), (-14, -14), (0, -1), (0, 0)], [(-14, -14), (-14, -14), (0, -1), (0, 0)], - [(-1, -1), (0, 0), (0, -1), (1, 0)], - [(0, 0), (2, 2), (0, 0), (2, 2)], + [(-1, -1), (0, 0), (0, -1), (0, 0)], + [(1, 0), (2, 2), (1, 0), (2, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]], constant_velocity: [[(-25, -25), (-25, -25), (0, -1), (1, 1)], [(-4, -5), (-3, -3), (0, -1), (1, 1)], diff --git a/pynumdiff/utils/utility.py b/pynumdiff/utils/utility.py index 77daf64..f7c3034 100644 --- a/pynumdiff/utils/utility.py +++ b/pynumdiff/utils/utility.py @@ -181,26 +181,27 @@ def slide_function(func, x, dt, kernel, *args, stride=1, pass_weights=False, **k if len(kernel) % 2 == 0: raise ValueError("Kernel window size should be odd.") half_window_size = (len(kernel) - 1)//2 # int because len(kernel) is always odd - weights = np.zeros((int(np.ceil(len(x)/stride)), len(x))) # Could be more space efficient - x_hats = np.zeros(weights.shape) - dxdt_hats = np.zeros(weights.shape) + x_hat = np.zeros(x.shape) + dxdt_hat = np.zeros(x.shape) + weight_sum = np.zeros(x.shape) for i,midpoint in enumerate(range(0, len(x), stride)): # iterate window midpoints # find where to index data and kernel, taking care at edges - window = slice(max(0, midpoint - half_window_size), - min(len(x), midpoint + half_window_size + 1)) # +1 because slicing is exclusive of end - kslice = slice(max(0, half_window_size - midpoint), - min(len(kernel), len(kernel) - (midpoint + half_window_size + 1 - len(x)))) + start = max(0, midpoint - half_window_size) + end = min(len(x), midpoint + half_window_size + 1) # +1 because slicing is exclusive of end + window = slice(start, end) - # weights need to be renormalized if running off an edge - weights[i, window] = kernel if kslice.stop - kslice.stop == len(kernel) else kernel[kslice]/np.sum(kernel[kslice]) - if pass_weights: kwargs['weights'] = weights[i, window] + kstart = max(0, half_window_size - midpoint) + kend = kstart + (end - start) + kslice = slice(kstart, kend) - # run the function on the window and save results - x_hats[i,window], dxdt_hats[i,window] = func(x[window], dt, *args, **kwargs) + w = kernel if (end-start) == len(kernel) else kernel[kslice]/np.sum(kernel[kslice]) + if pass_weights: kwargs['weights'] = w - weights /= weights.sum(axis=0, keepdims=True) # normalize the weights - x_hat = np.sum(weights*x_hats, axis=0) - dxdt_hat = np.sum(weights*dxdt_hats, axis=0) + # run the function on the window and add weighted results to cumulative answers + x_window_hat, dxdt_window_hat = func(x[window], dt, *args, **kwargs) + x_hat[window] += w * x_window_hat + dxdt_hat[window] += w * dxdt_window_hat + weight_sum[window] += w # save sum of weights for normalization at the end - return x_hat, dxdt_hat + return x_hat/weight_sum, dxdt_hat/weight_sum