moved smooth_finite_difference_tests and in the process found a way to enable options dictionaries in my old-style invocations

pavelkomarov · pavelkomarov · commit 7e6519e68713 · 2025-06-25T13:32:39.000-07:00
diff --git a/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py b/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py
@@ -170,14 +170,14 @@ def butterdiff(x, dt, params=None, options={}, filter_order=2, cutoff_freq=0.5,
     return x_hat, dxdt_hat
 
 
-def splinediff(x, dt, params=None, options={}, k=3, s=None, num_iterations=1):
+def splinediff(x, dt, params=None, options={}, order=3, s=None, num_iterations=1):
     """Perform spline smoothing on x with scipy.interpolate.UnivariateSpline followed by first order finite difference
 
     :param np.array[float] x: array of time series to differentiate
     :param float dt: time step size
     :param list params: (**deprecated**, prefer :code:`order`, :code:`cutoff_freq`, and :code:`num_iterations`)
     :param dict options: (**deprecated**, prefer :code:`num_iterations`) an empty dictionary or {'iterate': (bool)}
-    :param int k: polynomial order of the spline. A kth order spline can be differentiated k times.
+    :param int order: polynomial order of the spline. A kth order spline can be differentiated k times.
     :param float s: positive smoothing factor used to choose the number of knots. Number of knots will be increased
         until the smoothing condition is satisfied: :math:`\\sum_t (x[t] - \\text{spline}[t])^2 \\leq s`
     :param int num_iterations: how many times to apply smoothing
@@ -187,17 +187,17 @@ def splinediff(x, dt, params=None, options={}, k=3, s=None, num_iterations=1):
              - **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 `k`, `s`, and " +
+        warn("`params` and `options` parameters will be removed in a future version. Use `order`, `s`, and " +
             "`num_iterations` instead.", DeprecationWarning)
-        k, s = params[0:2]
+        order, s = params[0:2]
         if 'iterate' in options and options['iterate']:
             num_iterations = params[2]
 
     t = np.arange(0, len(x)*dt, dt)
 
     x_hat = x
     for _ in range(num_iterations):
-        spline = scipy.interpolate.UnivariateSpline(t, x_hat, k=k, s=s)
+        spline = scipy.interpolate.UnivariateSpline(t, x_hat, k=order, s=s)
         x_hat = spline(t)
 
     x_hat, dxdt_hat = finite_difference(x_hat, dt)
diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py
@@ -6,13 +6,13 @@
 from ..linear_model import lineardiff, polydiff, savgoldiff, spectraldiff
 from ..total_variation_regularization import velocity, acceleration, jerk, iterative_velocity
 from ..kalman_smooth import * # constant_velocity, constant_acceleration, constant_jerk, known_dynamics
-from ..smooth_finite_difference import * # mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff, splinediff
+from ..smooth_finite_difference import mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff, splinediff
 from ..finite_difference import first_order, second_order
 # Function aliases for testing cases where parameters change the behavior in a big way
-iterated_first_order = lambda *args, **kwargs: first_order(*args, **kwargs)
+def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs)
 
 dt = 0.1
-t = np.arange(0, 3, dt) # sample locations
+t = np.arange(0, 3+dt, dt) # sample locations, including the endpoint
 tt = np.linspace(0, 3) # full domain, for visualizing denser plots
 np.random.seed(7) # for repeatability of the test, so we don't get random failures
 noise = 0.05*np.random.randn(*t.shape)
@@ -22,70 +22,116 @@
     (0, r"$x(t)=1$",            lambda t: np.ones(t.shape), lambda t: np.zeros(t.shape)),   # constant
     (1, r"$x(t)=2t+1$",         lambda t: 2*t + 1,          lambda t: 2*np.ones(t.shape)),  # affine
     (2, r"$x(t)=t^2-t+1$",      lambda t: t**2 - t + 1,     lambda t: 2*t - 1),             # quadratic
-    (3, r"$x(t)=\sin(3t)+1/2$", lambda t: np.sin(3*t) + 1/2, lambda t: 3*np.cos(3*t)),     # sinuoidal
+    (3, r"$x(t)=\sin(3t)+1/2$", lambda t: np.sin(3*t) + 1/2, lambda t: 3*np.cos(3*t)),      # sinuoidal
     (4, r"$x(t)=e^t\sin(5t)$",  lambda t: np.exp(t)*np.sin(5*t),                            # growing sinusoidal
                                 lambda t: np.exp(t)*(5*np.cos(5*t) + np.sin(5*t))),
     (5, r"$x(t)=\frac{\sin(8t)}{(t+0.1)^{3/2}}$", lambda t: np.sin(8*t)/((t + 0.1)**(3/2)), # steep challenger
                                 lambda t: ((0.8 + 8*t)*np.cos(8*t) - 1.5*np.sin(8*t))/(0.1 + t)**(5/2))]
 
 # Call both ways, with kwargs (new) and with params list with default options dict (legacy), to ensure both work
 diff_methods_and_params = [
-    (first_order, None), (iterated_first_order, {'num_iterations':5}),
-    (second_order, None),
+    (first_order, {}),
+    (iterated_first_order, {'num_iterations':5}),
+    (second_order, {}), # empty dictionary for the case of no parameters
     #(lineardiff, {'order':3, 'gamma':5, 'window_size':10, 'solver':'CVXOPT'}),
     (polydiff, {'polynomial_order':2, 'window_size':3}), (polydiff, [2, 3]),
     (savgoldiff, {'polynomial_order':2, 'window_size':4, 'smoothing_win':4}), (savgoldiff, [2, 4, 4]),
-    (spectraldiff, {'high_freq_cutoff':0.1}), (spectraldiff, [0.1])
+    (spectraldiff, {'high_freq_cutoff':0.1}), (spectraldiff, [0.1]),
+    (mediandiff, {'window_size':3, 'num_iterations':2}), (mediandiff, [3, 2], {'iterate':True}),
+    (meandiff, {'window_size':3, 'num_iterations':2}), (meandiff, [3, 2], {'iterate':True}),
+    (gaussiandiff, {'window_size':5}), (gaussiandiff, [5]),
+    (friedrichsdiff, {'window_size':5}), (friedrichsdiff, [5]),
+    (butterdiff, {'filter_order':3, 'cutoff_freq':0.074}), (butterdiff, [3, 0.074]),
+    (splinediff, {'order':5, 's':2}), (splinediff, [5, 2])
     ]
 
 # All the testing methodology follows the exact same pattern; the only thing that changes is the
 # closeness to the right answer various methods achieve with the given parameterizations. So index a
 # big ol' table by the method, then the test function, then the pair of quantities we're comparing.
 error_bounds = {
     first_order: [[(-25, -25), (-25, -25), (0, 0), (1, 1)],
-                  [(-25, -25), (-14, -14), (0, 0), (1, 1)],
+                  [(-25, -25), (-13, -14), (0, 0), (1, 1)],
                   [(-25, -25), (0, 0), (0, 0), (1, 0)],
                   [(-25, -25), (0, 0), (0, 0), (1, 1)],
                   [(-25, -25), (2, 2), (0, 0), (2, 2)],
                   [(-25, -25), (3, 3), (0, 0), (3, 3)]],
-    iterated_first_order: [[(-7, -7), (-9, -10), (0, -1), (0, 0)],
-                           [(-5, -5), (-5, -6), (0, -1), (0, 0)],
+    iterated_first_order: [[(-8, -9), (-25, -25), (0, -1), (0, 0)],
+                           [(-6, -6), (-6, -7), (0, -1), (0, 0)],
                            [(-1, -1), (0, 0), (0, -1), (0, 0)],
-                           [(0, 0), (1, 1), (0, 0), (1, 1)],
+                           [(0, 0), (1, 0), (0, 0), (1, 0)],
                            [(1, 1), (2, 2), (1, 1), (2, 2)],
                            [(1, 1), (3, 3), (1, 1), (3, 3)]],
     second_order: [[(-25, -25), (-25, -25), (0, 0), (1, 1)],
-                   [(-25, -25), (-14, -14), (0, 0), (1, 1)],
-                   [(-25, -25), (-13, -14), (0, 0), (1, 1)],
+                   [(-25, -25), (-13, -13), (0, 0), (1, 1)],
+                   [(-25, -25), (-13, -13), (0, 0), (1, 1)],
                    [(-25, -25), (0, -1), (0, 0), (1, 1)],
                    [(-25, -25), (1, 1), (0, 0), (1, 1)],
                    [(-25, -25), (3, 3), (0, 0), (3, 3)]],
     #lineardiff: [TBD when #91 is solved],
-    polydiff: [[(-14, -15), (-14, -14), (0, -1), (1, 1)],
+    polydiff: [[(-15, -15), (-14, -14), (0, -1), (1, 1)],
+               [(-14, -14), (-13, -13), (0, -1), (1, 1)],
                [(-14, -14), (-13, -13), (0, -1), (1, 1)],
-               [(-14, -15), (-13, -14), (0, -1), (1, 1)],
                [(-2, -2), (0, 0), (0, -1), (1, 1)],
-               [(0, 0), (2, 1), (0, 0), (2, 1)],
+               [(0, 0), (1, 1), (0, -1), (1, 1)],
                [(0, 0), (3, 3), (0, 0), (3, 3)]],
-    savgoldiff: [[(-7, -8), (-13, -14), (0, -1), (0, 0)],
-                 [(-7, -8), (-13, -13), (0, -1), (0, 0)],
+    savgoldiff: [[(-9, -10), (-13, -14), (0, -1), (0, 0)],
+                 [(-9, -10), (-13, -13), (0, -1), (0, 0)],
                  [(-1, -1), (0, -1), (0, -1), (0, 0)],
                  [(0, -1), (0, 0), (0, -1), (1, 0)],
                  [(1, 1), (2, 2), (1, 1), (2, 2)],
                  [(1, 1), (3, 3), (1, 1), (3, 3)]],
-    spectraldiff: [[(-7, -8), (-14, -15), (-1, -1), (0, 0)],
+    spectraldiff: [[(-9, -10), (-14, -15), (-1, -1), (0, 0)],
+                   [(0, 0), (1, 1), (0, 0), (1, 1)],
+                   [(1, 0), (1, 1), (1, 1), (1, 1)],
                    [(0, 0), (1, 1), (0, 0), (1, 1)],
-                   [(1, 0), (1, 1), (1, 0), (1, 1)],
+                   [(1, 1), (2, 2), (1, 1), (2, 2)],
+                   [(1, 1), (3, 3), (1, 1), (3, 3)]],
+    mediandiff: [[(-25, -25), (-25, -25), (-1, -1), (0, 0)],
+                 [(0, 0), (1, 1), (0, 0), (1, 1)],
+                 [(0, 0), (1, 1), (0, 0), (1, 1)],
+                 [(-1, -1), (0, 0), (0, 0), (1, 1)],
+                 [(0, 0), (2, 2), (0, 0), (2, 2)],
+                 [(1, 1), (3, 3), (1, 1), (3, 3)]],
+    meandiff: [[(-25, -25), (-25, -25), (0, -1), (0, 0)],
+               [(0, 0), (1, 0), (0, 0), (1, 1)],
+               [(0, 0), (1, 1), (0, 0), (1, 1)],
+               [(0, 0), (1, 1), (0, 0), (1, 1)],
+               [(1, 1), (2, 2), (1, 1), (2, 2)],
+               [(1, 1), (3, 3), (1, 1), (3, 3)]],
+    gaussiandiff: [[(-14, -15), (-14, -14), (0, -1), (1, 0)],
+                   [(-1, -1), (0, 0), (0, 0), (1, 0)],
                    [(0, 0), (1, 1), (0, 0), (1, 1)],
+                   [(0, -1), (1, 0), (0, 0), (1, 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)]],
+    friedrichsdiff: [[(-25, -25), (-25, -25), (0, -1), (0, 0)],
+                     [(-1, -1), (0, 0), (0, 0), (1, 0)],
+                     [(0, 0), (1, 1), (0, 0), (1, 1)],
+                     [(0, -1), (1, 1), (0, 0), (1, 1)],
+                     [(1, 1), (2, 2), (1, 1), (2, 2)],
+                     [(1, 1), (3, 3), (1, 1), (3, 3)]],
+    butterdiff: [[(-13, -14), (-13, -13), (0, -1), (0, 0)],
+                 [(0, -1), (0, 0), (0, -1), (0, 0)],
+                 [(0, 0), (1, 1), (0, 0), (1, 1)],
+                 [(1, 0), (1, 1), (1, 0), (1, 1)],
+                 [(2, 2), (3, 2), (2, 2), (3, 2)],
+                 [(2, 1), (3, 3), (2, 1), (3, 3)]],
+    splinediff: [[(-14, -15), (-14, -14), (-1, -1), (0, 0)],
+                 [(-14, -14), (-13, -14), (-1, -1), (0, 0)],
+                 [(-14, -14), (0, 0), (-1, -1), (0, 0)],
+                 [(0, 0), (1, 1), (0, 0), (1, 1)],
+                 [(1, 0), (2, 2), (1, 0), (2, 2)],
+                 [(1, 0), (3, 3), (1, 0), (3, 3)]]
 }
 
+# Essentially run the cartesian product of [diff methods] x [test functions] through this one test
 @mark.filterwarnings("ignore::DeprecationWarning") # I want to test the old and new functionality intentionally
 @mark.parametrize("diff_method_and_params", diff_methods_and_params)
 @mark.parametrize("test_func_and_deriv", test_funcs_and_derivs)
 def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # request gives access to context
-    diff_method, params = diff_method_and_params # unpack
+    # unpack
+    diff_method, params = diff_method_and_params[:2]
+    if len(diff_method_and_params) == 3: options = diff_method_and_params[2] # optionally pass old-style `options` dict
     i, latex_name, f, df = test_func_and_deriv
 
     # some methods rely on cvxpy, and we'd like to allow use of pynumdiff without convex optimization
@@ -97,27 +143,30 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re
     x_noisy = x + noise # add a little noise
     dxdt = df(t) # true values of the derivative at samples
 
-    # differentiate without and with noise
-    x_hat, dxdt_hat = diff_method(x, dt, **params) if isinstance(params, dict) else diff_method(x, dt, params) \
-        if isinstance(params, list) else diff_method(x, dt)
+    # differentiate without and with noise, accounting for new and old styles of calling functions
+    x_hat, dxdt_hat = diff_method(x, dt, **params) if isinstance(params, dict) \
+        else diff_method(x, dt, params) if (isinstance(params, list) and len(diff_method_and_params) < 3) \
+        else diff_method(x, dt, params, options)
     x_hat_noisy, dxdt_hat_noisy = diff_method(x_noisy, dt, **params) if isinstance(params, dict) \
-        else diff_method(x_noisy, dt, params) if isinstance(params, list) else diff_method(x_noisy, dt)
+        else diff_method(x_noisy, dt, params) if (isinstance(params, list) and len(diff_method_and_params) < 3) \
+        else diff_method(x_noisy, dt, params, options)
     
     # check x_hat and x_hat_noisy are close to x and that dxdt_hat and dxdt_hat_noisy are close to dxdt
     #print("]\n[", end="")
     for j,(a,b) in enumerate([(x,x_hat), (dxdt,dxdt_hat), (x,x_hat_noisy), (dxdt,dxdt_hat_noisy)]):
-        l2_error = np.linalg.norm(a - b)
-        linf_error = np.max(np.abs(a - b))
+        #l2_error = np.linalg.norm(a - b)
+        #linf_error = np.max(np.abs(a - b))
+        #print(f"({l2_error},{linf_error})", end=", ")
+        #print(f"({int(np.ceil(np.log10(l2_error))) if l2_error > 0 else -25}, {int(np.ceil(np.log10(linf_error))) if linf_error > 0 else -25})", end=", ")
         
-        #print(f"({int(np.ceil(np.log10(l2_error))) if l2_error> 0 else -25}, {int(np.ceil(np.log10(linf_error))) if linf_error > 0 else -25})", end=", ")
         log_l2_bound, log_linf_bound = error_bounds[diff_method][i][j]
         assert np.linalg.norm(a - b) < 10**log_l2_bound
         assert np.max(np.abs(a - b)) < 10**log_linf_bound
         if 0 < np.linalg.norm(a - b) < 10**(log_l2_bound - 1) or 0 < np.max(np.abs(a - b)) < 10**(log_linf_bound - 1):
-            print(f"Improvement detected for method {str(diff_method)}")
+            print(f"Improvement detected for method {diff_method.__name__}")
 
-    if request.config.getoption("--plot"): # Get the plot flag from pytest configuration
-        fig, axes = request.config.plots[diff_method]
+    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
         axes[i, 0].plot(t, f(t), label=r"$x(t)$")
         axes[i, 0].plot(t, x, 'C0+')
         axes[i, 0].plot(tt, df(tt), label=r"$\frac{dx(t)}{dt}$")
diff --git a/pynumdiff/tests/test_smooth_finite_difference.py b/pynumdiff/tests/test_smooth_finite_difference.py
