finishing touches post claude

Author: pavelkomarov

pynumdiff/kalman_smooth.py

@@ -111,7 +111,8 @@ def rtsdiff(x, dt_or_t, order, log_qr_ratio, forwardbackward):
     """Perform Rauch-Tung-Striebel smoothing with a naive constant derivative model. Makes use of :code:`kalman_filter`
     and :code:`rts_smooth`, which are made public. :code:`constant_X` methods in this module call this function.
 
-    :param np.array[float] x: data series to differentiate
+    :param np.array[float] x: data series to differentiate. May contain NaN values (missing data); NaNs are excluded from
+        fitting and imputed by dynamical model evolution. 
     :param float or array[float] dt_or_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 int order: which derivative to stabilize in the constant-derivative model

pynumdiff/polynomial_fit.py

@@ -51,11 +51,10 @@ def polydiff(x, dt_or_t, params=None, options=None, degree=None, window_size=Non
     kernel='friedrichs'):
     """Fit polynomials to the data, and differentiate the polynomials.
 
-    :param np.array[float] x: data to differentiate. May contain NaN values (missing data);
-        NaNs are excluded from fitting and imputed by polynomial interpolation.
-    :param float or np.array[float] dt_or_t: This function supports variable step size. This
-        parameter is either the constant :math:`\\Delta t` if given as a single float, or data
-        locations if given as an array of same length as :code:`x`.
+    :param np.array[float] x: data to differentiate. May contain NaN values (missing data); NaNs are excluded from
+        fitting and imputed by polynomial interpolation.
+    :param float or array[float] dt_or_t: This function supports variable step size. This parameter is either the constant
+        :math:`\\Delta t` if given as a single float, or data locations if given as an array of same length as :code:`x`.
     :param list[int] params: (**deprecated**, prefer :code:`degree` 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)}
@@ -86,15 +85,9 @@ def polydiff(x, dt_or_t, params=None, options=None, degree=None, window_size=Non
         warn("Kernel window size should be odd. Added 1 to length.")
 
     def _polydiff(x, dt_or_t, degree, weights=None):
-        # Build time array: either from scalar dt or from passed array of locations
-        if np.isscalar(dt_or_t):
-            t = np.arange(len(x)) * dt_or_t
-        else:
-            t = np.asarray(dt_or_t)
-
-        # Filter out NaN values so polyfit doesn't fail on missing data
-        mask = ~np.isnan(x)
-        if not np.any(mask): return np.full_like(x, np.nan), np.full_like(x, np.nan) # all NaN window
+        t = dt_or_t if not np.isscalar(dt_or_t) else np.arange(len(x)) * dt_or_t # sample locations
+        mask = ~np.isnan(x) # Filter out any NaN values so polyfit doesn't lose its mind in the event of missing data
+        if not np.any(mask): warn("Window of all NaNs encountered. `polyfit` will fail. Choose a wider `window_size`?")
 
         r = np.polyfit(t[mask], x[mask], degree, w=weights[mask] if weights is not None else None) # polyfit returns highest order first
         dr = np.polyder(r) # power rule already implemented for us
@@ -104,8 +97,7 @@ def _polydiff(x, dt_or_t, degree, weights=None):
 
         return x_hat, dxdt_hat
 
-    if not window_size:
-        return _polydiff(x, dt_or_t, degree)
+    if not window_size: return _polydiff(x, dt_or_t, degree)
 
     kernel = {'gaussian':utility.gaussian_kernel, 'friedrichs':utility.friedrichs_kernel}[kernel](window_size)
     return utility.slide_function(_polydiff, x, dt_or_t, kernel, degree, stride=step_size, pass_weights=True)

pynumdiff/tests/test_diff_methods.py

@@ -14,6 +14,7 @@ def iterated_second_order(*args, **kwargs): return second_order(*args, **kwargs)
 def iterated_fourth_order(*args, **kwargs): return fourth_order(*args, **kwargs)
 def spline_irreg_step(*args, **kwargs): return splinediff(*args, **kwargs)
 def polydiff_irreg_step(*args, **kwargs): return polydiff(*args, **kwargs)
+irreg_list = [spline_irreg_step, polydiff_irreg_step, rbfdiff, rtsdiff] # methods to test with irregular time steps
 
 dt = 0.1
 t = np.linspace(0, 3, 31) # sample locations, including the endpoint
@@ -135,7 +136,7 @@ def polydiff_irreg_step(*args, **kwargs): return polydiff(*args, **kwargs)
                [(0, 0), (1, 1), (0, -1), (1, 1)],
                [(0, 0), (3, 3), (0, 0), (3, 3)]],
     polydiff_irreg_step: [[(-14, -15), (-14, -14), (0, -1), (1, 1)],
-                          [(-14, -14), (-13, -14), (0, -1), (1, 1)],
+                          [(-14, -14), (-13, -13), (0, -1), (1, 1)],
                           [(-14, -14), (-13, -13), (0, -1), (1, 1)],
                           [(-2, -2), (0, 0), (0, -1), (1, 1)],
                           [(0, 0), (1, 1), (0, 0), (1, 1)],
@@ -250,9 +251,9 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re
     i, latex_name, f, df = test_func_and_deriv
 
     # sample the true function and true derivative, and make noisy samples
-    x = f(t) if diff_method not in [spline_irreg_step, polydiff_irreg_step, rbfdiff, rtsdiff] else f(t_irreg)
-    dxdt = df(t) if diff_method not in [spline_irreg_step, polydiff_irreg_step, rbfdiff, rtsdiff] else df(t_irreg)
-    _t = dt if diff_method not in [spline_irreg_step, polydiff_irreg_step, rbfdiff, rtsdiff] else t_irreg
+    x = f(t) if diff_method not in irreg_list else f(t_irreg)
+    dxdt = df(t) if diff_method not in irreg_list else df(t_irreg)
+    _t = dt if diff_method not in irreg_list else t_irreg
     x_noisy = x + noise
 
     # differentiate without and with noise, accounting for new and old styles of calling functions
@@ -266,7 +267,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, polydiff_irreg_step, rtsdiff, rbfdiff] else t
+        t_ = t_irreg if diff_method in irreg_list 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)
@@ -381,6 +382,3 @@ def test_multidimensionality(multidim_method_and_params, request):
         ax3.plot_wireframe(T1, T2, computed_laplacian, label='computed')
         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)
-
-
-

pynumdiff/tests/test_utils.py

@@ -4,7 +4,7 @@
 
 from pynumdiff.utils import utility, evaluate
 from pynumdiff.utils.simulate import sine, triangle, pop_dyn, linear_autonomous, pi_cruise_control, lorenz_x
-np.random.seed(42) # The answer to life, the universe, and everything
+np.random.seed(7)
 
 
 def test_integrate_dxdt_hat():
@@ -55,9 +55,9 @@ def test_convolutional_smoother():
 
 def test_peakdet(request):
     """Verify peakdet finds peaks and valleys"""
-    rng = np.random.RandomState(42) # own seed so test order doesn't affect results
     t = np.arange(0, 10, 0.001)
-    x = 0.3*np.sin(t) + np.sin(1.3*t) + 0.9*np.sin(4.2*t) + 0.02*rng.randn(10000)
+    x = 0.3*np.sin(t) + np.sin(1.3*t) + 0.9*np.sin(4.2*t) + \
+        0.02*np.random.RandomState(42).randn(10000) # isolated source of randomness so test order doesn't affect results
     maxtab, mintab = utility.peakdet(x, 0.5, t)
 
     if request.config.getoption("--plot"):
@@ -88,13 +88,11 @@ def identity(x, dt_or_t): return x, 0 # should come back the same
     x = np.arange(100, dtype=float)
     kernel = utility.gaussian_kernel(9)
 
-    # Scalar dt: existing behavior should be unchanged
-    x_hat, dxdt_hat = utility.slide_function(identity, x, 0.1, kernel, stride=2)
-    assert np.allclose(x, x_hat)
+    x_hat_dt, _ = utility.slide_function(identity, x, 0.1, kernel, stride=2)
+    assert np.allclose(x, x_hat_dt)
 
-    # Array t: func receives a t-slice instead of scalar dt; identity still returns x unchanged
-    t = np.linspace(0, 10, 100)
-    x_hat_t, _ = utility.slide_function(identity, x, t, kernel, stride=2)
+    # time array: func receives a kernel-length slice of times instead of scalar dt; identity still returns x unchanged
+    x_hat_t, _ = utility.slide_function(identity, x, np.linspace(0, 10, 100, endpoint=False), kernel, stride=2)
     assert np.allclose(x, x_hat_t)
 
 

pynumdiff/utils/utility.py

@@ -120,23 +120,20 @@ def slide_function(func, x, dt_or_t, kernel, *args, stride=1, pass_weights=False
 
     :param callable func: name of the function to slide
     :param np.array[float] x: data to differentiate
-    :param float or np.array[float] dt_or_t: constant step size (scalar) or array of sample
-        locations (same length as x). When given as an array, the actual time values for each
-        window are passed to func, enabling support for nonuniform spacing.
+    :param float or np.array[float] dt_or_t: constant step size (scalar) or array of sample locations (same length as x).
+        When given as an array, the actual time values for each window are passed to :code:`func`, enabling nonuniform spacing.
     :param np.array[float] kernel: values to weight the sliding window
-    :param list args: passed to func
-    :param int stride: step size for slide (e.g. 1 means slide by 1 step)
+    :param list args: passed to :code:`func`
+    :param int stride: step size for slide (e.g. 1 means slide by 1 index location)
     :param bool pass_weights: whether weights should be passed to func via update to kwargs
-    :param dict kwargs: passed to func
+    :param dict kwargs: passed to :code:`func`
 
     :return: - **x_hat** -- estimated (smoothed) x
              - **dxdt_hat** -- estimated derivative of x
     """
     if len(kernel) % 2 == 0: raise ValueError("Kernel window size should be odd.")
     half_window_size = (len(kernel) - 1)//2 # int because len(kernel) is always odd
-
-    scalar_dt = np.isscalar(dt_or_t)
-    t = None if scalar_dt else np.asarray(dt_or_t)
+    equispaced = np.isscalar(dt_or_t)
 
     x_hat = np.zeros(x.shape)
     dxdt_hat = np.zeros(x.shape)
@@ -146,7 +143,7 @@ def slide_function(func, x, dt_or_t, kernel, *args, stride=1, pass_weights=False
         # find where to index data and kernel, taking care at edges
         start = max(0, midpoint - half_window_size)
         end = min(len(x), midpoint + half_window_size + 1) # +1 because slicing is exclusive of end
-        window = slice(start, end)
+        window = slice(start, end) # This is in terms of indices, not true time in the event of nonuniform spacing
 
         kstart = max(0, half_window_size - midpoint)
         kend = kstart + (end - start)
@@ -155,12 +152,8 @@ def slide_function(func, x, dt_or_t, kernel, *args, stride=1, pass_weights=False
         w = kernel if (end-start) == len(kernel) else kernel[kslice]/np.sum(kernel[kslice])
         if pass_weights: kwargs['weights'] = w
 
-        # When t is an array, pass the actual time values for the window so func can handle
-        # nonuniform spacing; otherwise pass dt as before (backward compatible).
-        dt_or_t_window = dt_or_t if scalar_dt else t[window]
-
-        # run the function on the window and add weighted results to cumulative answers
-        x_window_hat, dxdt_window_hat = func(x[window], dt_or_t_window, *args, **kwargs)
+        # Run the function on the window and add weighted results to cumulative answers. If not equispaced, pass times for window.
+        x_window_hat, dxdt_window_hat = func(x[window], dt_or_t if equispaced else dt_or_t[window], *args, **kwargs)
         x_hat[window] += w * x_window_hat
         dxdt_hat[window] += w * dxdt_window_hat
         weight_sum[window] += w # save sum of weights for normalization at the end