diff --git a/pynumdiff/__init__.py b/pynumdiff/__init__.py index 729f8cf..70975a9 100644 --- a/pynumdiff/__init__.py +++ b/pynumdiff/__init__.py @@ -7,4 +7,4 @@ from .polynomial_fit import splinediff, polydiff, savgoldiff from .total_variation_regularization import tvrdiff, velocity, acceleration, jerk, iterative_velocity, smooth_acceleration, jerk_sliding from .kalman_smooth import rtsdiff, constant_velocity, constant_acceleration, constant_jerk, known_dynamics -from .linear_model import spectraldiff, lineardiff \ No newline at end of file +from .linear_model import spectraldiff, lineardiff, rbfdiff \ No newline at end of file diff --git a/pynumdiff/linear_model/__init__.py b/pynumdiff/linear_model/__init__.py index b8097a5..387c686 100644 --- a/pynumdiff/linear_model/__init__.py +++ b/pynumdiff/linear_model/__init__.py @@ -8,8 +8,8 @@ warn("Limited Linear Model Support Detected! CVXPY is not installed. " + "Install CVXPY to use lineardiff derivatives. You can still use other methods.") -from ._linear_model import savgoldiff, polydiff, spectraldiff +from ._linear_model import savgoldiff, polydiff, spectraldiff, rbfdiff -__all__ = ['lineardiff', 'spectraldiff'] # So these get treated as direct members of the module by sphinx -# polydiff and savgoldiff are still imported for now for backwards compatibility but are not -# documented as part of this module, since they've moved +__all__ = ['lineardiff', 'spectraldiff', 'rbfdiff'] # So these get treated as direct members of the +# module by sphinx polydiff and savgoldiff are still imported for now for backwards compatibility +# but are not documented as part of this module, since they've moved diff --git a/pynumdiff/linear_model/_linear_model.py b/pynumdiff/linear_model/_linear_model.py index 2122fe5..dfb065b 100644 --- a/pynumdiff/linear_model/_linear_model.py +++ b/pynumdiff/linear_model/_linear_model.py @@ -1,5 +1,6 @@ import math, scipy import numpy as np +from scipy import sparse from warnings import warn from pynumdiff.finite_difference import second_order as finite_difference @@ -22,9 +23,6 @@ def polydiff(*args, **kwargs): return _polydiff(*args, **kwargs) -############### -# Linear diff # -############### def _solve_for_A_and_C_given_X_and_Xdot(X, Xdot, num_integrations, dt, gammaC=1e-1, gammaA=1e-6, solver=None, A_known=None, epsilon=1e-6): """Given state and the derivative, find the system evolution and measurement matrices.""" @@ -161,9 +159,6 @@ def _lineardiff(x, dt, order, gamma, solver=None): return x_hat, dxdt_hat -############################### -# Fourier Spectral derivative # -############################### def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_extension=True, pad_to_zero_dxdt=True): """Take a derivative in the fourier domain, with high frequency attentuation. @@ -233,3 +228,53 @@ def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_e x_hat = x_hat + x0 return x_hat, dxdt_hat + + +def rbfdiff(x, _t, sigma=1, lmbd=0.01): + """Find smoothed function and derivative estimates by fitting noisy data against a radial-basis-function-filled + kernel (naively NxN but made sparse by truncating tiny values). This function uses a Gaussian kernel, thereby + assuming the signal is generated from a Gaussian process, effectively the same assumption the Kalman filter makes. + It can handle variable step size just as easily as uniform step size. + + :param np.array[float] x: data to differentiate + :param float or array[float] _t: This function supports variable step size. This parameter is either the constant + step size if given as a single float, or data locations if given as an array of same length as :code:`x`. + :param float sigma: controls width of radial basis function + :param float lmbd: controls strength of bias toward data + + :return: tuple[np.array, np.array] of\n + - **x_hat** -- estimated (smoothed) x + - **dxdt_hat** -- estimated derivative of x + """ + if isinstance(_t, (np.ndarray, list)): # support variable step size for this function + if len(x) != len(_t): raise ValueError("If `_t` is given as array-like, must have same length as `x`.") + t = _t + else: + t = np.arange(len(x))*_t + + # The below does the approximate equivalent of this code, but sparsely in O(N sigma), since the rbf falls off rapidly + # t_i, t_j = np.meshgrid(t,t) + # r = t_j - t_i # radius + # rbf = np.exp(-(r**2) / (2 * sigma**2)) # radial basis function kernel + # drbfdt = -(r / sigma**2) * rbf # derivative of kernel + # rbf_regularized = rbf + lmbd*np.eye(len(t)) + # alpha = np.linalg.solve(rbf_regularized, x) + + cutoff = np.sqrt(-2 * sigma**2 * np.log(1e-4)) + rows, cols, vals, dvals = [], [], [], [] + for n in range(len(t)): + # Only consider points within a cutoff. Gaussian drops below eps at distance ~ sqrt(-2*sigma^2 log eps) + l = np.searchsorted(t, t[n] - cutoff) # O(log N) to find indices of points within cutoff + r = np.searchsorted(t, t[n] + cutoff) # finds index where new value should be inserted + for j in range(l, r): # width of this is dependent on sigma. [l, r) is correct inclusion/exclusion + radius = t[n] - t[j] + v = np.exp(-radius**2 / (2 * sigma**2)) + dv = -radius / sigma**2 * v + rows.append(n); cols.append(j); vals.append(v); dvals.append(dv) + + rbf = sparse.csr_matrix((vals, (rows, cols)), shape=(len(t), len(t))) # Build sparse kernels + drbfdt = sparse.csr_matrix((dvals, (rows, cols)), shape=(len(t), len(t))) + rbf_regularized = rbf + lmbd*sparse.eye(len(t), format="csr") + alpha = sparse.linalg.spsolve(rbf_regularized, x) # solve sparse system + + return rbf @ alpha, drbfdt @ alpha diff --git a/pynumdiff/polynomial_fit/_polynomial_fit.py b/pynumdiff/polynomial_fit/_polynomial_fit.py index 9e319a8..4bbcdf0 100644 --- a/pynumdiff/polynomial_fit/_polynomial_fit.py +++ b/pynumdiff/polynomial_fit/_polynomial_fit.py @@ -5,13 +5,13 @@ from pynumdiff.utils import utility -def splinediff(x, dt, params=None, options={}, degree=3, s=None, num_iterations=1): +def splinediff(x, _t, params=None, options={}, degree=3, s=None, num_iterations=1): """Find smoothed data and derivative estimates by fitting a smoothing spline to the data with - scipy.interpolate.UnivariateSpline. + scipy.interpolate.UnivariateSpline. Variable step size is supported with equal ease as uniform step size. :param np.array[float] x: data to differentiate - :param float or array[float] dt: This function supports variable step size. This parameter is either the constant - step size if given as a single float, or data locations if given as an array of same length as :code:`x`. + :param float or array[float] _t: This function supports variable step size. This parameter is either the constant + step size if given as a single float, or data locations if given as an array of same length as :code:`x`. :param list params: (**deprecated**, prefer :code:`degree`, :code:`cutoff_freq`, and :code:`num_iterations`) :param dict options: (**deprecated**, prefer :code:`num_iterations`) an empty dictionary or {'iterate': (bool)} :param int degree: polynomial degree of the spline. A kth degree spline can be differentiated k times. @@ -30,11 +30,11 @@ def splinediff(x, dt, params=None, options={}, degree=3, s=None, num_iterations= if 'iterate' in options and options['iterate']: num_iterations = params[2] - if isinstance(dt, (np.ndarray, list)): # support variable step size for this function - if len(x) != len(dt): raise ValueError("If `dt` is given as array-like, must have same length as `x`.") - t = dt + if isinstance(_t, (np.ndarray, list)): # support variable step size for this function + if len(x) != len(_t): raise ValueError("If `_t` is given as array-like, must have same length as `x`.") + t = _t else: - t = np.arange(len(x))*dt + t = np.arange(len(x))*_t x_hat = x for _ in range(num_iterations): diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index f6a001a..bcd006f 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -3,7 +3,7 @@ from warnings import warn from ..finite_difference import first_order, second_order, fourth_order -from ..linear_model import lineardiff, spectraldiff +from ..linear_model import lineardiff, spectraldiff, rbfdiff from ..polynomial_fit import polydiff, savgoldiff, splinediff from ..total_variation_regularization import velocity, acceleration, jerk, iterative_velocity, smooth_acceleration, jerk_sliding from ..kalman_smooth import constant_velocity, constant_acceleration, constant_jerk @@ -54,7 +54,8 @@ def spline_irreg_step(*args, **kwargs): return splinediff(*args, **kwargs) (smooth_acceleration, {'gamma':2, 'window_size':5}), (smooth_acceleration, [2, 5]), (jerk_sliding, {'gamma':1, 'window_size':15}), (jerk_sliding, [1], {'window_size':15}), (spectraldiff, {'high_freq_cutoff':0.2}), (spectraldiff, [0.2]), - (lineardiff, {'order':3, 'gamma':5, 'window_size':11, 'solver':'CLARABEL'}), (lineardiff, [3, 5, 11], {'solver':'CLARABEL'}) + (lineardiff, {'order':3, 'gamma':5, 'window_size':11, 'solver':'CLARABEL'}), (lineardiff, [3, 5, 11], {'solver':'CLARABEL'}), + (rbfdiff, {'sigma':0.5, 'lmbd':0.001}) ] # All the testing methodology follows the exact same pattern; the only thing that changes is the @@ -212,7 +213,13 @@ def spline_irreg_step(*args, **kwargs): return splinediff(*args, **kwargs) [(1, 0), (2, 2), (1, 0), (2, 2)], [(1, 0), (2, 1), (1, 0), (2, 1)], [(1, 1), (2, 2), (1, 1), (2, 2)], - [(1, 1), (3, 3), (1, 1), (3, 3)]] + [(1, 1), (3, 3), (1, 1), (3, 3)]], + rbfdiff: [[(-2, -2), (0, 0), (0, -1), (0, 0)], + [(-1, -1), (1, 0), (0, -1), (1, 1)], + [(-1, -1), (1, 1), (0, -1), (1, 1)], + [(-2, -2), (0, 0), (0, -1), (0, 0)], + [(0, 0), (2, 2), (0, 0), (2, 2)], + [(1, 1), (3, 3), (1, 1), (3, 3)]] } # Essentially run the cartesian product of [diff methods] x [test functions] through this one test @@ -231,7 +238,7 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re except: warn(f"Cannot import cvxpy, skipping {diff_method} test."); return # sample the true function and true derivative, and make noisy samples - if diff_method in [spline_irreg_step]: # list that can handle variable dt + if diff_method in [spline_irreg_step, rbfdiff]: # list that can handle variable dt x = f(t_irreg) dxdt = df(t_irreg) _t = t_irreg @@ -252,7 +259,7 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re # plotting code if request.config.getoption("--plot") and not isinstance(params, list): # Get the plot flag from pytest configuration fig, axes = request.config.plots[diff_method] # get the appropriate plot, set up by the store_plots fixture in conftest.py - t_ = t_irreg if diff_method in [spline_irreg_step] else t + t_ = t_irreg if diff_method in [spline_irreg_step, rbfdiff] else t axes[i, 0].plot(t_, f(t_)) axes[i, 0].plot(t_, x, 'C0+') axes[i, 0].plot(t_, x_hat, 'C2.', ms=4)