Merge pull request #149 from florisvb/add-rbfdiff

Making the house even bigger

Author: pavelkomarov

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
+from .linear_model import spectraldiff, lineardiff, rbfdiff

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

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

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

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)