diff --git a/pynumdiff/basis_fit.py b/pynumdiff/basis_fit.py index e6ae307..32fd679 100644 --- a/pynumdiff/basis_fit.py +++ b/pynumdiff/basis_fit.py @@ -34,7 +34,7 @@ def spectraldiff(x, dt, cutoff_freq, extension='odd', pad_to_flat=False, axis=0) N = x.shape[axis] x = np.moveaxis(x, axis, 0) - x_flat = x.reshape(N, -1) # flat 2D of all the vectors to differentiate + x_flat = x.reshape(N, -1).astype(float, copy=False) # flat 2D of all the vectors to differentiate, float so extension can't overflow pad = 0 if pad_to_flat: # repeat the end values outward, smooth the joins, then restore the original in the middle @@ -163,7 +163,7 @@ def waveletdiff(x, dt, wavelet='db8', mode='symmetric', threshold=2.0, level=Non N = x.shape[axis] x = np.moveaxis(x, axis, 0) - x_flat = x.reshape(N, -1) + x_flat = x.reshape(N, -1).astype(float, copy=False) # float so extension arithmetic can't overflow # Step 0: Lock down adaptive parameters max_level = pywt.dwt_max_level(N, wavelet) # structural ceiling: how many halvings the filter still fits in diff --git a/pynumdiff/finite_difference.py b/pynumdiff/finite_difference.py index 7dc170a..f8ebd68 100644 --- a/pynumdiff/finite_difference.py +++ b/pynumdiff/finite_difference.py @@ -23,7 +23,7 @@ def finitediff(x, dt, num_iterations=1, order=2, axis=0): if order not in [1, 2, 4]: raise ValueError("order must be 1, 2, or 4") x = np.moveaxis(x, axis, 0) # move the axis of interest to the front to simplify differencing indexing - x_hat = np.asarray(x) # allows for array-like. Preserve reference to x, for finding the final constant of integration + x_hat = x.astype(float, copy=False) # Preserve reference to x, for finding the final constant of integration; float so no overflow possible dxdt_hat = np.zeros(x.shape) # preallocate reusable memory # For all but the last iteration, do the differentate->integrate smoothing loop, being careful with endpoints diff --git a/pynumdiff/smooth_finite_difference.py b/pynumdiff/smooth_finite_difference.py index a5b84f3..a5a8908 100644 --- a/pynumdiff/smooth_finite_difference.py +++ b/pynumdiff/smooth_finite_difference.py @@ -58,7 +58,7 @@ def butterdiff(x, dt, filter_order=2, cutoff_freq=0.5, num_iterations=1, axis=0) sos = scipy.signal.butter(filter_order, cutoff_freq, output='sos') # second-order sections rather than the (b, a) transfer # function, whose coefficients lose all precision when many poles bunch up near z = 1 at high order and low cutoff - x_hat = x + x_hat = x.astype(float, copy=False) # float so sosfiltfilt's odd-extension padding can't overflow padlen = min(3*(filter_order + 1), x.shape[axis]-1) # scipy's own default, can overwhelm short data for _ in range(num_iterations): x_hat = scipy.signal.sosfiltfilt(sos, x_hat, axis=axis, padlen=padlen) # applies forward and backward pass so zero phase diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index a550431..1f39044 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -9,7 +9,7 @@ from ..total_variation_regularization import iterative_velocity, smooth_acceleration, tvrdiff from ..kalman_smooth import rtsdiff, robustdiff from ..linear_model import lineardiff -# Function aliases for testing cases where parameters change the behavior in a big way, so error limits can be indexed in dict +# Function aliases for testing cases where parameters change behavior in a big way, so error limits can be indexed in dict def kerneldiff_uniform(*args, **kwargs): return kerneldiff(*args, **kwargs) def kerneldiff_median(*args, **kwargs): return kerneldiff(*args, **kwargs) def kerneldiff_gaussian(*args, **kwargs): return kerneldiff(*args, **kwargs) @@ -256,7 +256,6 @@ def lineardiff_irreg_step(*args, **kwargs): return lineardiff(*args, **kwargs) [(0, -1), (2, 2), (0, 0), (2, 2)]] } - # Essentially run the cartesian product of [diff methods] x [test functions] through this one test @mark.parametrize("diff_method_and_params", diff_methods_and_params) # things like splinediff, with their parameters @mark.parametrize("test_func_and_deriv", test_funcs_and_derivs) # analytic functions, with their true derivatives @@ -328,7 +327,6 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re dt2 = 0.02 # distance between samples in the 2D T grids x = T1**2 * np.sin(3/2 * np.pi * T2) # 2D function - multidim_methods_and_params = [ # single parameterizations, chosen to do well on the noiseless 2D example (kerneldiff, {'kernel': 'gaussian', 'window_size': 5}), (butterdiff, {'filter_order': 3, 'cutoff_freq': 1 - 1e-6}), @@ -417,6 +415,13 @@ def test_multidimensionality(multidim_method_and_params, request): legend = ax3.legend(bbox_to_anchor=(0.7, 0.8)); legend.legend_handles[0].set_facecolor(pyplot.cm.viridis(0.6)) fig.suptitle(f'{diff_method.__name__}', fontsize=16) +@mark.parametrize("multidim_method_and_params", multidim_methods_and_params) +def test_integer_inputs(multidim_method_and_params): + """Ensure integer samples give the same derivative as floats, rather than overflowing or rounding""" + diff_method, params = multidim_method_and_params + n = np.arange(60); x_int = np.stack([200 - 3*n, 20 + 3*n], axis=1).astype(np.uint8) # falling and rising, so integer arithmetic wraps + assert np.allclose(diff_method(x_int, dt, **params)[1], diff_method(x_int.astype(float), dt, **params)[1]) + nan_methods_and_params = [ # List of methods that can handle missing values (splinediff, {'degree': 5, 's': 2}), diff --git a/pynumdiff/utils/utility.py b/pynumdiff/utils/utility.py index e4e579a..0f14908 100644 --- a/pynumdiff/utils/utility.py +++ b/pynumdiff/utils/utility.py @@ -142,8 +142,8 @@ def convolutional_smoother(x, kernel, num_iterations=1, axis=0): :return: **x_hat** (np.array[float]) -- smoothed x """ x_hat = x - for i in range(num_iterations): - x_hat = convolve1d(x_hat, kernel, axis=axis, mode='reflect') # 'reflect' pads the signal with mirrored repeats + for i in range(num_iterations): # 'reflect' pads the signal with mirrored repeats, and float output to + x_hat = convolve1d(x_hat, kernel, axis=axis, mode='reflect', output=float) # avoid truncation return x_hat