From d43be08b10be339ee6cec3233950934c18d22049 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Tue, 1 Jul 2025 18:04:28 -0700 Subject: [PATCH] investigated #116 and added some comments to help clarify the code --- pynumdiff/kalman_smooth/_kalman_smooth.py | 2 +- pynumdiff/tests/test_diff_methods.py | 4 ++-- .../_total_variation_regularization.py | 21 +++++++------------ 3 files changed, 11 insertions(+), 16 deletions(-) diff --git a/pynumdiff/kalman_smooth/_kalman_smooth.py b/pynumdiff/kalman_smooth/_kalman_smooth.py index f68a923..dbc09ad 100644 --- a/pynumdiff/kalman_smooth/_kalman_smooth.py +++ b/pynumdiff/kalman_smooth/_kalman_smooth.py @@ -36,7 +36,7 @@ def _kalman_forward_filter(xhat0, P0, y, A, C, Q, R, u=None, B=None): xhat_pre.append(xhat_) # the [n]th index holds _{n|n-1} info P_pre.append(P_) - xhat = xhat_.copy() + xhat = xhat_.copy() # copies very important here. See #122 P = P_.copy() if not np.isnan(y[n]): # handle missing data K = P_ @ C.T @ np.linalg.inv(C @ P_ @ C.T + R) diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index 6aaf8a9..846c7a1 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -185,10 +185,10 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) [(0, 0), (1, 0), (0, -1), (1, 0)], [(1, 1), (2, 2), (1, 1), (2, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]], - jerk_sliding: [[(-15, -15), (-16, -17), (0, -1), (1, 0)], + jerk_sliding: [[(-15, -15), (-16, -16), (0, -1), (1, 0)], [(-14, -14), (-14, -14), (0, -1), (0, 0)], [(-14, -14), (-14, -14), (0, -1), (0, 0)], - [(-1, -1), (0, 0), (0, -1), (0, 0)], + [(-1, -1), (0, 0), (0, -1), (1, 0)], [(0, 0), (2, 2), (0, 0), (2, 2)], [(1, 1), (3, 3), (1, 1), (3, 3)]] } diff --git a/pynumdiff/total_variation_regularization/_total_variation_regularization.py b/pynumdiff/total_variation_regularization/_total_variation_regularization.py index 56ba776..540df9e 100644 --- a/pynumdiff/total_variation_regularization/_total_variation_regularization.py +++ b/pynumdiff/total_variation_regularization/_total_variation_regularization.py @@ -69,7 +69,7 @@ def _total_variation_regularized_derivative(x, dt, order, gamma, solver=None): - **x_hat** -- estimated (smoothed) x - **dxdt_hat** -- estimated derivative of x """ - # Normalize + # Normalize for numerical consistency with convex solver mean = np.mean(x) std = np.std(x) if std == 0: std = 1 # safety guard @@ -79,7 +79,8 @@ def _total_variation_regularized_derivative(x, dt, order, gamma, solver=None): deriv_values = cvxpy.Variable(len(x)) # values of the order^th derivative, in which we're penalizing variation integration_constants = cvxpy.Variable(order) # constants of integration that help get us back to x - # Recursively integrate the highest order derivative to get back to the position + # Recursively integrate the highest order derivative to get back to the position. This is a first- + # order scheme, but it's very fast and tends to not do markedly worse than 2nd order. See #116 y = deriv_values for i in range(order): y = cvxpy.cumsum(y) + integration_constants[i] @@ -87,7 +88,7 @@ def _total_variation_regularized_derivative(x, dt, order, gamma, solver=None): # Set up and solve the optimization problem prob = cvxpy.Problem(cvxpy.Minimize( # Compare the recursively integrated position to the noisy position, and add TVR penalty - cvxpy.sum_squares(y - x) + cvxpy.sum(gamma*cvxpy.tv(deriv_values)) )) + cvxpy.sum_squares(y - x) + gamma*cvxpy.sum(cvxpy.tv(deriv_values)) )) prob.solve(solver=solver) # Recursively integrate the final derivative values to get back to the function and derivative values @@ -97,14 +98,10 @@ def _total_variation_regularized_derivative(x, dt, order, gamma, solver=None): dxdt_hat = y/dt # y only holds the dx values; to get deriv scale by dt x_hat = np.cumsum(y) + integration_constants.value[order-1] # smoothed data - dxdt_hat = (dxdt_hat[:-1] + dxdt_hat[1:])/2 # take first order FD to smooth a touch - ddxdt_hat_f = dxdt_hat[-1] - dxdt_hat[-2] - dxdt_hat_f = dxdt_hat[-1] + ddxdt_hat_f # What is this doing? Could we use a 2nd order FD above natively? - dxdt_hat = np.hstack((dxdt_hat, dxdt_hat_f)) - - # fix first point - d = dxdt_hat[2] - dxdt_hat[1] - dxdt_hat[0] = dxdt_hat[1] - d + # Due to the first-order nature of the derivative, it has a slight lag. Average together every two values + # to better center the answer. But this leaves us one-short, so devise a good last value. + dxdt_hat = (dxdt_hat[:-1] + dxdt_hat[1:])/2 + dxdt_hat = np.hstack((dxdt_hat, 2*dxdt_hat[-1] - dxdt_hat[-2])) # last value = penultimate value [-1] + diff between [-1] and [-2] return x_hat*std+mean, dxdt_hat*std # derivative is linear, so scale derivative by std @@ -262,5 +259,3 @@ def jerk_sliding(x, dt, params=None, options=None, gamma=None, solver=None, wind ramp = window_size//5 kernel = np.hstack((np.arange(1, ramp+1)/ramp, np.ones(window_size - 2*ramp), (np.arange(1, ramp+1)/ramp)[::-1])) return utility.slide_function(_total_variation_regularized_derivative, x, dt, kernel, 3, gamma, stride=ramp, solver=solver) - -