smoothing over

Author: pavelkomarov

pynumdiff/basis_fit.py

@@ -43,9 +43,11 @@ def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_e
         x = np.concatenate((pre, x, post), axis=axis) # extend the edges
         kernel = utility.mean_kernel(padding//2)
         x_smoothed = utility.convolutional_smoother(x, kernel, axis=axis) # smooth the padded edges in
-        center = (slice(None),)*axis + (slice(padding, L+padding),) + (slice(None),)*(x.ndim-axis-1)
-        x_smoothed[center] = x[center] # restore original signal in the middle
+        original_signal = (slice(None),)*axis + (slice(padding, L+padding),) + (slice(None),)*(x.ndim-axis-1)
+        x_smoothed[original_signal] = x[original_signal] # restore original signal in the middle
         x = x_smoothed
+    else:
+        padding = 0
 
     # Do even extension (optional)
     if even_extension is True:
@@ -60,8 +62,8 @@ def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_e
 
     # Filter to zero out higher wavenumbers
     discrete_cutoff = int(high_freq_cutoff * N / 2) # Nyquist is at N/2 location, and we're cutting off as a fraction of that
-    filt = np.ones(k.shape)  # start with all frequencies passing
-    filt[discrete_cutoff:-discrete_cutoff] = 0  # zero out high-frequency components
+    filt = np.ones(k.shape) # start with all frequencies passing
+    filt[discrete_cutoff:-discrete_cutoff] = 0 # zero out high-frequency components
 
     # Smoothed signal
     X = np.fft.fft(x, axis=axis)
@@ -71,7 +73,8 @@ def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_e
     omega = 2*np.pi/(dt*N) # factor of 2pi/T turns wavenumbers into frequencies in radians/s
     dxdt = np.real(np.fft.ifft(1j * k[s] * omega * filt[s] * X, axis=axis))
 
-    return (x_hat[center], dxdt[center]) if pad_to_zero_dxdt else (x_hat, dxdt)
+    original_signal = (slice(None),)*axis + (slice(padding, L+padding),) + (slice(None),)*(x_hat.ndim-axis-1)
+    return x_hat[original_signal], dxdt[original_signal]
 
 
 def rbfdiff(x, dt_or_t, sigma=1, lmbd=0.01, axis=0):
@@ -90,19 +93,19 @@ def rbfdiff(x, dt_or_t, sigma=1, lmbd=0.01, axis=0):
     :return: - **x_hat** (np.array) -- estimated (smoothed) x
              - **dxdt_hat** (np.array) -- estimated derivative of x
     """
-    x = np.asarray(x)
-    x = np.moveaxis(x, axis, 0)   # bring target axis to front
-    orig_shape = x.shape
-    N = orig_shape[0]
-    x_2d = x.reshape(N, -1)       # (N, M) — build matrix once, solve for all M columns
+    N = x.shape[axis]
+    x = np.moveaxis(x, axis, 0) # bring axis of differentiation to front so each N repeats comprise vector
+    plump = x.shape
+    x_flattened = x.reshape(N, -1) # (N, M) matrix where each column is a vector along the original axis
 
     if np.isscalar(dt_or_t):
         t = np.arange(N)*dt_or_t
     else: # support variable step size for this function
         if N != len(dt_or_t): raise ValueError("If `dt_or_t` is given as array-like, must have same length as `x`.")
         t = dt_or_t
 
-    # The below does the approximate equivalent of this code, but sparsely in O(N sigma^2), since the rbf falls off rapidly
+    # For each vector along the axis of differentiation, the below does the approximate equivalent of this code,
+    # but sparsely in O(N sigma^2), 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, O(N^2) entries
@@ -125,8 +128,9 @@ def rbfdiff(x, dt_or_t, sigma=1, lmbd=0.01, axis=0):
     rbf = sparse.csr_matrix((vals, (rows, cols)), shape=(N, N)) # Build sparse kernels, O(N sigma) entries
     drbfdt = sparse.csr_matrix((dvals, (rows, cols)), shape=(N, N))
     rbf_regularized = rbf + lmbd*sparse.eye(N, format="csr") # identity matrix gives a little extra height at the centers
-    alpha = sparse.linalg.spsolve(rbf_regularized, x_2d) # solve sparse system targeting the noisy data, O(N sigma^2)
+    alpha = sparse.linalg.spsolve(rbf_regularized, x_flattened) # solve sparse system targeting the noisy data,
+                                                                # can take matrix target, O(N sigma^2) for each vector
+    x_hat_flattened = rbf @ alpha # find samples of reconstructions using the smooth bases
+    dxdt_hat_flattened = drbfdt @ alpha
 
-    x_hat = np.moveaxis((rbf @ alpha).reshape(orig_shape), 0, axis)     # find samples of reconstructions using the smooth bases
-    dxdt_hat = np.moveaxis((drbfdt @ alpha).reshape(orig_shape), 0, axis)
-    return x_hat, dxdt_hat
+    return np.moveaxis(x_hat_flattened.reshape(plump), 0, axis), np.moveaxis(dxdt_hat_flattened.reshape(plump), 0, axis)

pynumdiff/tests/test_diff_methods.py

@@ -328,7 +328,8 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re
     (polydiff, {'degree': 2, 'window_size': 5}),
     (savgoldiff, {'degree': 3, 'window_size': 11, 'smoothing_win': 3}),
     (rtsdiff, {'order':2, 'log_qr_ratio':7, 'forwardbackward':True}),
-    (rbfdiff, {'sigma': 0.5, 'lmbd': 0.001}),
+    (spectraldiff, {'high_freq_cutoff': 0.25, 'pad_to_zero_dxdt': False}),
+    (rbfdiff, {'sigma': 0.5, 'lmbd': 1e-6}),
     (splinediff, {'degree': 9, 's': 1e-6}),
     (robustdiff, {'order':2, 'log_q':7, 'log_r':2})
 ]
@@ -344,7 +345,8 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re
     polydiff: [(1, -1), (1, 0)],
     savgoldiff: [(0, -1), (1, 1)],
     rtsdiff: [(1, -1), (1, 0)],
-    rbfdiff: [(2, 1), (2, 1)],
+    spectraldiff: [(2, 1), (3, 2)], # lot of Gibbs ringing in 2nd order derivatives along t1 with t_1^2 sin(3 pi t_2 / 2)
+    rbfdiff: [(0, -1), (1, 0)],
     splinediff: [(0, -1), (1, 0)],
     robustdiff: [(-2, -3), (0, -1)]
 }