new polynomial_fit module

docs/source/code.rst

@@ -6,6 +6,7 @@ API documentation
 
    finite_difference
    smooth_finite_difference
+   polynomial_fit
    total_variation_regularization
    kalman_smooth
    linear_model

docs/source/polynomial_fit.rst

@@ -0,0 +1,5 @@
+polynomial_fit
+==============
+
+.. automodule:: pynumdiff.polynomial_fit
+    :members:

docs/source/utils/evaluate.rst

@@ -1,7 +1,5 @@
 evaluate
 ========
 
-.. currentmodule:: pynumdiff.utils.evaluate
-
 .. automodule:: pynumdiff.utils.evaluate
     :members:

docs/source/utils/simulate.rst

@@ -1,7 +1,5 @@
 simulate
 ========
 
-.. currentmodule:: pynumdiff.utils.simulate
-
 .. automodule:: pynumdiff.utils.simulate
     :members:

docs/source/utils/utility.rst

@@ -1,7 +1,5 @@
 utility
 =======
 
-.. currentmodule:: pynumdiff.utils.utility
-
 .. automodule:: pynumdiff.utils.utility
     :members:

pynumdiff/__init__.py

@@ -1,11 +1,10 @@
 """Import useful functions from all modules
 """
-from pynumdiff._version import __version__
-from pynumdiff.finite_difference import first_order, second_order, fourth_order
-from pynumdiff.smooth_finite_difference import mediandiff, meandiff, gaussiandiff,\
-    friedrichsdiff, butterdiff, splinediff
-from pynumdiff.total_variation_regularization import iterative_velocity, velocity,\
-    acceleration, jerk, smooth_acceleration, jerk_sliding
-from pynumdiff.linear_model import lineardiff, polydiff, spectraldiff, savgoldiff
-from pynumdiff.kalman_smooth import constant_velocity, constant_acceleration, constant_jerk,\
-    known_dynamics
+from ._version import __version__
+
+from .finite_difference import finitediff, first_order, second_order, fourth_order
+from .smooth_finite_difference import meandiff, mediandiff, gaussiandiff, friedrichsdiff, butterdiff
+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 rts_const_deriv, constant_velocity, constant_acceleration, constant_jerk, known_dynamics
+from .linear_model import spectraldiff, lineardiff

pynumdiff/finite_difference/__init__.py

@@ -1,3 +1,3 @@
 """This module implements some common finite difference schemes
 """
-from ._finite_difference import finite_difference, first_order, second_order, fourth_order
+from ._finite_difference import finitediff, first_order, second_order, fourth_order

pynumdiff/finite_difference/_finite_difference.py

@@ -4,7 +4,7 @@
 from warnings import warn
 
 
-def finite_difference(x, dt, num_iterations, order):
+def finitediff(x, dt, num_iterations, order):
     """Perform iterated finite difference of a given order. This serves as the common backing function for
     all other methods in this module.
 
@@ -85,7 +85,7 @@ def first_order(x, dt, params=None, options={}, num_iterations=1):
         warn("`params` and `options` parameters will be removed in a future version. Use `num_iterations` instead.", DeprecationWarning)
         num_iterations = params[0] if isinstance(params, list) else params
 
-    return finite_difference(x, dt, num_iterations, 1)
+    return finitediff(x, dt, num_iterations, 1)
 
 
 def second_order(x, dt, num_iterations=1):
@@ -100,7 +100,7 @@ def second_order(x, dt, num_iterations=1):
              - **x_hat** -- original x if :code:`num_iterations=1`, else smoothed x that yielded dxdt_hat
              - **dxdt_hat** -- estimated derivative of x
     """
-    return finite_difference(x, dt, num_iterations, 2)
+    return finitediff(x, dt, num_iterations, 2)
 
 
 def fourth_order(x, dt, num_iterations=1):
@@ -115,4 +115,4 @@ def fourth_order(x, dt, num_iterations=1):
              - **x_hat** -- original x if :code:`num_iterations=1`, else smoothed x that yielded dxdt_hat
              - **dxdt_hat** -- estimated derivative of x
     """
-    return finite_difference(x, dt, num_iterations, 4)
+    return finitediff(x, dt, num_iterations, 4)

pynumdiff/linear_model/__init__.py

@@ -10,4 +10,6 @@
 
 from ._linear_model import savgoldiff, polydiff, spectraldiff
 
-__all__ = ['lineardiff', 'savgoldiff', 'polydiff', 'spectraldiff'] # So these get treated as direct members of the module by sphinx
+__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

pynumdiff/linear_model/_linear_model.py

@@ -1,115 +1,25 @@
-import copy, math, logging, scipy
+import math, scipy
 import numpy as np
 from warnings import warn
 
-from pynumdiff.finite_difference import first_order as finite_difference
+from pynumdiff.finite_difference import second_order as finite_difference
+from pynumdiff.polynomial_fit import savgoldiff as _savgoldiff # patch through
+from pynumdiff.polynomial_fit import polydiff as _polydiff # patch through
 from pynumdiff.utils import utility
 
 try: import cvxpy
 except ImportError: pass
 
 
-#########################
-# Savitzky-Golay filter #
-#########################
-def savgoldiff(x, dt, params=None, options=None, poly_order=None, window_size=None, smoothing_win=None):
-    """Use the Savitzky-Golay to smooth the data and calculate the first derivative. It uses
-    scipy.signal.savgol_filter. The Savitzky-Golay is very similar to the sliding polynomial fit,
-    but slightly noisier, and much faster.
+def savgoldiff(*args, **kwargs):
+    warn("`savgoldiff` has moved to `polynomial_fit.savgoldiff` and will be removed from "
+        + "`linear_model` in a future release.", DeprecationWarning)
+    return _savgoldiff(*args, **kwargs)
 
-    :param np.array[float] x: data to differentiate
-    :param float dt: step size
-    :param list params: (**deprecated**, prefer :code:`poly_order`, :code:`window_size`, and :code:`smoothing_win`)
-    :param dict options: (**deprecated**)
-    :param int poly_order: order of the polynomial
-    :param int window_size: size of the sliding window, must be odd (if not, 1 is added)
-    :param int smoothing_win: size of the window used for gaussian smoothing, a good default is
-        window_size, but smaller for high frequnecy data
-
-    :return: tuple[np.array, np.array] of\n
-             - **x_hat** -- estimated (smoothed) x
-             - **dxdt_hat** -- estimated derivative of x
-    """
-    if params != None: # Warning to support old interface for a while. Remove these lines along with params in a future release.
-        warn("`params` and `options` parameters will be removed in a future version. Use `poly_order`, " +
-            "`window_size`, and `smoothing_win` instead.", DeprecationWarning)
-        poly_order, window_size, smoothing_win = params
-    elif poly_order == None or window_size == None or smoothing_win == None:
-        raise ValueError("`poly_order`, `window_size`, and `smoothing_win` must be given.")
-
-    window_size = np.clip(window_size, poly_order + 1, len(x)-1)
-    if window_size % 2 == 0:
-        window_size += 1 # window_size needs to be odd
-        warn("Kernel window size should be odd. Added 1 to length.")
-    smoothing_win = min(smoothing_win, len(x)-1)
-
-    dxdt_hat = scipy.signal.savgol_filter(x, window_size, poly_order, deriv=1)/dt
-
-    kernel = utility.gaussian_kernel(smoothing_win)
-    dxdt_hat = utility.convolutional_smoother(dxdt_hat, kernel)
-
-    x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt)
-    x0 = utility.estimate_integration_constant(x, x_hat)
-    x_hat = x_hat + x0
-
-    return x_hat, dxdt_hat
-
-
-######################
-# Polynomial fitting #
-######################
-def polydiff(x, dt, params=None, options=None, poly_order=None, window_size=None, step_size=1,
-    kernel='friedrichs'):
-    """Fit polynomials to the data, and differentiate the polynomials.
-
-    :param np.array[float] x: data to differentiate
-    :param float dt: step size
-    :param list[int] params: (**deprecated**, prefer :code:`poly_order` and :code:`window_size`)
-    :param dict options: (**deprecated**, prefer :code:`step_size` and :code:`kernel`)
-            a dictionary consisting of {'sliding': (bool), 'step_size': (int), 'kernel_name': (str)}
-    :param int poly_order: order of the polynomial
-    :param int window_size: size of the sliding window, if not given no sliding
-    :param int step_size: step size for sliding
-    :param str kernel: name of kernel to use for weighting and smoothing windows ('gaussian' or 'friedrichs')
-
-    :return: tuple[np.array, np.array] of\n
-             - **x_hat** -- estimated (smoothed) x
-             - **dxdt_hat** -- estimated derivative of x
-    """
-    if params != None:
-        warn("`params` and `options` parameters will be removed in a future version. Use `poly_order` " +
-            "and `window_size` instead.", DeprecationWarning)
-        poly_order = params[0]
-        if len(params) > 1: window_size = params[1]
-        if options != None:
-            if 'sliding' in options and not options['sliding']: window_size = None
-            if 'step_size' in options: step_size = options['step_size']
-            if 'kernel_name' in options: kernel = options['kernel_name']
-    elif poly_order == None or window_size == None:
-        raise ValueError("`poly_order` and `window_size` must be given.")
-
-    if window_size < poly_order*3:
-        window_size = poly_order*3+1
-    if window_size % 2 == 0:
-        window_size += 1
-        warn("Kernel window size should be odd. Added 1 to length.")
-
-    def _polydiff(x, dt, poly_order, weights=None):
-        t = np.arange(len(x))*dt
-
-        r = np.polyfit(t, x, poly_order, w=weights) # polyfit returns highest order first
-        dr = np.polyder(r) # power rule already implemented for us
-
-        dxdt_hat = np.polyval(dr, t) # evaluate the derivative and original polynomials at points t
-        x_hat = np.polyval(r, t) # smoothed x
-
-        return x_hat, dxdt_hat
-
-    if not window_size:
-        return _polydiff(x, dt, poly_order)
-
-    kernel = {'gaussian':utility.gaussian_kernel, 'friedrichs':utility.friedrichs_kernel}[kernel](window_size)
-    return utility.slide_function(_polydiff, x, dt, kernel, poly_order, stride=step_size, pass_weights=True)
+def polydiff(*args, **kwargs):
+    warn("`polydiff` has moved to `polynomial_fit.polydiff` and will be removed from "
+        + "`linear_model` in a future release.", DeprecationWarning)
+    return _polydiff(*args, **kwargs)
 
 
 ###############

pynumdiff/optimize/_optimize.py

@@ -8,9 +8,10 @@
 from tqdm import tqdm
 
 from ..utils import evaluate
-from ..finite_difference import finite_difference, first_order, second_order, fourth_order
-from ..smooth_finite_difference import mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff, splinediff
-from ..linear_model import spectraldiff, polydiff, savgoldiff, lineardiff
+from ..finite_difference import finitediff, first_order, second_order, fourth_order
+from ..smooth_finite_difference import mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff
+from ..linear_model import spectraldiff, lineardiff
+from ..polynomial_fit import polydiff, savgoldiff, splinediff
 from ..total_variation_regularization import tvrdiff, velocity, acceleration, jerk, iterative_velocity, smooth_acceleration, jerk_sliding
 from ..kalman_smooth import rts_const_deriv, constant_velocity, constant_acceleration, constant_jerk
 
@@ -41,9 +42,9 @@
                  {'order': (1, 5),
                   'gamma': (1e-3, 1000),
                   'window_size': (15, 1000)}),
-    finite_difference: ({'num_iterations': [5, 10, 30, 50],
-                         'order': {2, 4}}, # order is categorical here, because it can't be 3
-                        {'num_iterations': (1, 1000)}),
+    finitediff: ({'num_iterations': [5, 10, 30, 50],
+                  'order': {2, 4}}, # order is categorical here, because it can't be 3
+                 {'num_iterations': (1, 1000)}),
     first_order: ({'num_iterations': [5, 10, 30, 50]},
                   {'num_iterations': (1, 1000)}),
     mediandiff: ({'window_size': [5, 15, 30, 50],
@@ -236,7 +237,7 @@ def suggest_method(x, dt, dxdt_truth=None, cutoff_frequency=None):
             raise ValueError('Either dxdt_truth or cutoff_frequency must be provided.')
         tvgamma = np.exp(-1.6*np.log(cutoff_frequency) -0.71*np.log(dt) - 5.1) # See https://ieeexplore.ieee.org/document/9241009
 
-    methods = [finite_difference, mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff,
+    methods = [finitediff, mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff,
         splinediff, spectraldiff, polydiff, savgoldiff, rts_const_deriv]
     try: # optionally skip some methods
         import cvxpy

pynumdiff/polynomial_fit/__init__.py

@@ -0,0 +1,5 @@
+"""Find a smooth fit with some kind of polynomials, and then take derivative.
+"""
+from ._polynomial_fit import splinediff, polydiff, savgoldiff
+
+__all__ = ['splinediff', 'polydiff', 'savgoldiff'] # So automodule from the .rst finds them