From aa252ba8308de189a410ee14885c1c96d3a14dbd Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Mon, 23 Jun 2025 19:55:48 -0700 Subject: [PATCH 1/5] beginning to touch TVR --- pynumdiff/tests/test_diff_methods.py | 18 ++++++++++-------- .../__chartrand_tvregdiff__.py | 4 ++-- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index 5e9796d..7178bef 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -37,6 +37,7 @@ (savgoldiff, {'polynomial_order':2, 'window_size':4, 'smoothing_win':4}), (savgoldiff, [2, 4, 4]), (spectraldiff, {'high_freq_cutoff':0.1}), (spectraldiff, [0.1]) ] +diff_methods_and_params = [(velocity, [0.5])] # 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 @@ -78,7 +79,8 @@ [(1, 0), (1, 1), (1, 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)]] + [(1, 1), (3, 3), (1, 1), (3, 3)]], + velocity: None } @mark.filterwarnings("ignore::DeprecationWarning") # I want to test the old and new functionality intentionally @@ -104,17 +106,17 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re else diff_method(x_noisy, dt, params) if isinstance(params, list) else diff_method(x_noisy, dt) # 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="") + 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)) - #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"({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 {diff_method.__name__}") if request.config.getoption("--plot"): # Get the plot flag from pytest configuration fig, axes = request.config.plots[diff_method] diff --git a/pynumdiff/total_variation_regularization/__chartrand_tvregdiff__.py b/pynumdiff/total_variation_regularization/__chartrand_tvregdiff__.py index 2e226a6..d64d075 100644 --- a/pynumdiff/total_variation_regularization/__chartrand_tvregdiff__.py +++ b/pynumdiff/total_variation_regularization/__chartrand_tvregdiff__.py @@ -133,9 +133,9 @@ # # # (1) # # scipy.sparse.linalg.cg seems to require more # -# iterations to reach the same result as # +# iterations to reach the same result as # # MATLAB's pcg. (scipy version 1.1.0) # -# A good guess is the length of the data # +# A good guess is the length of the data # # # # (2) # # Drop last entry of derivative (u) for small scale # From 5aa7cab2440a85a4bb447b9e951e0daa66b52f33 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Thu, 26 Jun 2025 14:32:44 -0700 Subject: [PATCH 2/5] went around removing triple quotes from warnings, because they preserve newlines and print weird --- .../finite_difference/_finite_difference.py | 2 +- pynumdiff/linear_model/__init__.py | 4 +- pynumdiff/linear_model/_linear_model.py | 16 ++-- pynumdiff/tests/test_diff_methods.py | 18 ++-- .../__init__.py | 8 +- ...tvregdiff__.py => _chartrand_tvregdiff.py} | 0 .../_total_variation_regularization.py | 86 +++++++++---------- 7 files changed, 64 insertions(+), 70 deletions(-) rename pynumdiff/total_variation_regularization/{__chartrand_tvregdiff__.py => _chartrand_tvregdiff.py} (100%) diff --git a/pynumdiff/finite_difference/_finite_difference.py b/pynumdiff/finite_difference/_finite_difference.py index 8232736..3397e14 100644 --- a/pynumdiff/finite_difference/_finite_difference.py +++ b/pynumdiff/finite_difference/_finite_difference.py @@ -18,7 +18,7 @@ def first_order(x, dt, params=None, options={}, num_iterations=None): - **dxdt_hat** -- estimated derivative of x """ if params != None and 'iterate' in options: - warn("""`params` and `options` parameters will be removed in a future version. Use `num_iterations` instead.""", DeprecationWarning) + warn("`params` and `options` parameters will be removed in a future version. Use `num_iterations` instead.", DeprecationWarning) if isinstance(params, list): params = params[0] return _iterate_first_order(x, dt, params) elif num_iterations: diff --git a/pynumdiff/linear_model/__init__.py b/pynumdiff/linear_model/__init__.py index 933d074..0118fd4 100644 --- a/pynumdiff/linear_model/__init__.py +++ b/pynumdiff/linear_model/__init__.py @@ -5,8 +5,8 @@ from ._linear_model import lineardiff except: from warnings import warn - warn("""Limited Linear Model Support Detected! CVXPY is not installed. - Install CVXPY to use lineardiff derivatives. You can still use other methods.""") + warn("Limited Linear Model Support Detected! CVXPY is not installed. " + + "Install CVXPY to use lineardiff derivatives. You can still use other methods.") from ._linear_model import savgoldiff, polydiff, spectraldiff diff --git a/pynumdiff/linear_model/_linear_model.py b/pynumdiff/linear_model/_linear_model.py index d67d11b..1fd73a9 100644 --- a/pynumdiff/linear_model/_linear_model.py +++ b/pynumdiff/linear_model/_linear_model.py @@ -112,8 +112,8 @@ def savgoldiff(x, dt, params=None, options=None, polynomial_order=None, window_s - **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 `polynomial_order`, - `window_size`, and `smoothing_win` instead.""", DeprecationWarning) + warn("`params` and `options` parameters will be removed in a future version. Use `polynomial_order`, " + + "`window_size`, and `smoothing_win` instead.", DeprecationWarning) polynomial_order, window_size, smoothing_win = params elif polynomial_order == None or window_size == None or smoothing_win == None: raise ValueError("`polynomial_order`, `window_size`, and `smoothing_win` must be given.") @@ -192,8 +192,8 @@ def polydiff(x, dt, params=None, options=None, polynomial_order=None, window_siz - **dxdt_hat** -- estimated derivative of x """ if params != None: - warn("""`params` and `options` parameters will be removed in a future version. Use `polynomial_order` - and `window_size` instead.""", DeprecationWarning) + warn("`params` and `options` parameters will be removed in a future version. Use `polynomial_order` " + + "and `window_size` instead.", DeprecationWarning) polynomial_order, window_size = params if options != None: if 'sliding' in options: sliding = options['sliding'] @@ -426,8 +426,8 @@ def lineardiff(x, dt, params=None, options=None, order=None, gamma=None, window_ - **dxdt_hat** -- estimated derivative of x """ if params != None: - warn("""`params` and `options` parameters will be removed in a future version. Use `order`, - `gamma`, and `window_size` instead.""", DeprecationWarning) + warn("`params` and `options` parameters will be removed in a future version. Use `order`, " + + "`gamma`, and `window_size` instead.", DeprecationWarning) order, gamma, window_size = params if options != None: if 'sliding' in options: sliding = options['sliding'] @@ -483,8 +483,8 @@ def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_e - **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 `high_freq_cutoff`, - `even_extension`, and `pad_to_zero_dxdt` instead.""", DeprecationWarning) + warn("`params` and `options` parameters will be removed in a future version. Use `high_freq_cutoff`, " + + "`even_extension`, and `pad_to_zero_dxdt` instead.", DeprecationWarning) high_freq_cutoff = params[0] if isinstance(params, list) else params if options != None: if 'even_extension' in options: even_extension = options['even_extension'] diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index 5ec93a6..8c6e150 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -44,7 +44,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) (butterdiff, {'filter_order':3, 'cutoff_freq':0.074}), (butterdiff, [3, 0.074]), (splinediff, {'order':5, 's':2}), (splinediff, [5, 2]) ] -diff_methods_and_params = [(velocity, [0.5])] +diff_methods_and_params = [(velocity, [0.5], {'solver': 'CVXOPT'})] # 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 @@ -155,16 +155,16 @@ def test_diff_method(diff_method_and_params, test_func_and_deriv, request): # re # 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 {diff_method.__name__}") + # 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 {diff_method.__name__}") 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 diff --git a/pynumdiff/total_variation_regularization/__init__.py b/pynumdiff/total_variation_regularization/__init__.py index 4729519..8037587 100644 --- a/pynumdiff/total_variation_regularization/__init__.py +++ b/pynumdiff/total_variation_regularization/__init__.py @@ -5,10 +5,10 @@ from ._total_variation_regularization import velocity, acceleration, jerk, jerk_sliding, smooth_acceleration except: from warnings import warn - warn("""Limited Total Variation Regularization Support Detected! CVXPY is not installed. - Many Total Variation Methods require CVXPY including: velocity, acceleration, jerk, jerk_sliding, smooth_acceleration. - Please install CVXPY to use these methods. Recommended to also install MOSEK and obtain a MOSEK license. - You can still use: total_variation_regularization.iterative_velocity.""") + warn("Limited Total Variation Regularization Support Detected! CVXPY is not installed. " + + "Many Total Variation Methods require CVXPY including: velocity, acceleration, jerk, jerk_sliding, smooth_acceleration. " + + "Please install CVXPY to use these methods. Recommended to also install MOSEK and obtain a MOSEK license. " + + "You can still use: total_variation_regularization.iterative_velocity.") from ._total_variation_regularization import iterative_velocity diff --git a/pynumdiff/total_variation_regularization/__chartrand_tvregdiff__.py b/pynumdiff/total_variation_regularization/_chartrand_tvregdiff.py similarity index 100% rename from pynumdiff/total_variation_regularization/__chartrand_tvregdiff__.py rename to pynumdiff/total_variation_regularization/_chartrand_tvregdiff.py diff --git a/pynumdiff/total_variation_regularization/_total_variation_regularization.py b/pynumdiff/total_variation_regularization/_total_variation_regularization.py index 04bff1a..96e0baf 100644 --- a/pynumdiff/total_variation_regularization/_total_variation_regularization.py +++ b/pynumdiff/total_variation_regularization/_total_variation_regularization.py @@ -1,8 +1,7 @@ import logging import numpy as np -from pynumdiff.total_variation_regularization import __chartrand_tvregdiff__ -import pynumdiff.smooth_finite_difference +from pynumdiff.total_variation_regularization import _chartrand_tvregdiff from pynumdiff.utils import utility try: @@ -10,50 +9,46 @@ except ImportError: pass -# Iterative total variation regularization -def iterative_velocity(x, dt, params, options=None): - """ - Use an iterative solver to find the total variation regularized 1st derivative. - See __chartrand_tvregdiff__.py for details, author info, and license - Methods described in: Rick Chartrand, "Numerical differentiation of noisy, nonsmooth data," - ISRN Applied Mathematics, Vol. 2011, Article ID 164564, 2011. - Original code (MATLAB and python): https://sites.google.com/site/dnartrahckcir/home/tvdiff-code - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: a list consisting of: - - - iterations: Number of iterations to run the solver. More iterations results in blockier derivatives, which approach the convex result - - gamma: Regularization parameter. - - :type params: list (int, float) - - :param options: a dictionary with 2 key value pairs - - - 'cg_maxiter': Max number of iterations to use in scipy.sparse.linalg.cg. Default is None, results in maxiter = len(x). This works well in our test examples. - - 'scale': This method has two different numerical options. From __chartrand_tvregdiff__.py: 'large' or 'small' (case insensitive). Default is 'small'. 'small' has somewhat better boundary behavior, but becomes unwieldly for data larger than 1000 entries or so. 'large' has simpler numerics but is more efficient for large-scale problems. 'large' is more readily modified for higher-order derivatives, since the implicit differentiation matrix is square. - :type options: dict {'cg_maxiter': (int), 'scale': (string)}, optional - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x - - :rtype: tuple -> (np.array, np.array) +def iterative_velocity(x, dt, params, options=None, num_iterations=None, gamma=None, cg_maxiter=1000, scale='small'): + """Use an iterative solver to find the total variation regularized 1st derivative. See + _chartrand_tvregdiff.py for details, author info, and license. Methods described in: + Rick Chartrand, "Numerical differentiation of noisy, nonsmooth data," ISRN Applied Mathematics, + Vol. 2011, Article ID 164564, 2011. Original code at https://sites.google.com/site/dnartrahckcir/home/tvdiff-code + + :param np.array[float] x: array of time series to differentiate + :param float dt: time step size + :param list params: (**deprecated**, prefer :code:`num_iterations` and :code:`gamma`) + :param dict options: (**deprecated**, prefer :code:`cg_maxiter` and :code:`scale`) + a dictionary consisting of {'cg_maxiter': (int), 'scale': (str)} + :param int num_iterations: number of iterations to run the solver. More iterations results in + blockier derivatives, which approach the convex result + :param float gamma: regularization parameter + :param int cg_maxiter: Max number of iterations to use in :code:`scipy.sparse.linalg.cg`. Default + :code:`None` results in maxiter = len(x). This works well in our test examples. + :param str scale: This method has two different numerical options. From :code:`_chartrand_tvregdiff.py`: + :code:`'large'` or :code:`'small'` (case insensitive). Default is :code:`'small'`. :code:`'small'` + has somewhat better boundary behavior, but becomes unwieldly for data larger than 1000 entries or so. + :code:`'large'` has simpler numerics but is more efficient for large-scale problems. :code:`'large'` + is more readily modified for higher-order derivatives, since the implicit differentiation matrix is square. + + :return: tuple[np.array, np.array] of\n + - **x_hat** -- estimated (smoothed) x + - **dxdt_hat** -- estimated derivative of x """ - - if options is None: - options = {'cg_maxiter': 1000, 'scale': 'small'} - - iterations, gamma = params - dxdt_hat = __chartrand_tvregdiff__.TVRegDiff(x, iterations, gamma, dx=dt, - maxit=options['cg_maxiter'], scale=options['scale'], - ep=1e-6, u0=None, plotflag=False, diagflag=1) + 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 `num_iterations`, " + + "`gamma`, `cg_maxiter`, and `scale` instead.", DeprecationWarning) + num_iterations, gamma = params + if options != None: + if 'cg_maxiter' in options: cg_maxiter = options['cg_maxiter'] + if 'scale' in options: scale = options['scale'] + elif num_iterations == None or gamma == None: + raise ValueError("`num_iterations` and `gamma` must be given.") + + dxdt_hat = _chartrand_tvregdiff.TVRegDiff(x, num_iterations, gamma, dx=dt, + maxit=cg_maxiter, scale=scale, + ep=1e-6, u0=None, plotflag=False, diagflag=1) x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt) x0 = utility.estimate_initial_condition(x, x_hat) x_hat = x_hat + x0 @@ -61,9 +56,8 @@ def iterative_velocity(x, dt, params, options=None): return x_hat, dxdt_hat -# Generalized total variation regularized derivatives def __total_variation_regularized_derivative__(x, dt, N, gamma, solver='MOSEK'): - """ + """Generalized total variation regularized derivatives Use convex optimization (cvxpy) to solve for the Nth total variation regularized derivative. Default solver is MOSEK: https://www.mosek.com/ From d902e7fa5d7a51eee80f7ad0bb3119cab322c055 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Thu, 26 Jun 2025 16:01:06 -0700 Subject: [PATCH 3/5] decided to finish #68 for TVR module instead of moving tests --- .../_total_variation_regularization.py | 361 ++++++++---------- 1 file changed, 151 insertions(+), 210 deletions(-) diff --git a/pynumdiff/total_variation_regularization/_total_variation_regularization.py b/pynumdiff/total_variation_regularization/_total_variation_regularization.py index 96e0baf..c649e44 100644 --- a/pynumdiff/total_variation_regularization/_total_variation_regularization.py +++ b/pynumdiff/total_variation_regularization/_total_variation_regularization.py @@ -1,13 +1,18 @@ -import logging import numpy as np +from warnings import warn from pynumdiff.total_variation_regularization import _chartrand_tvregdiff from pynumdiff.utils import utility +try:import cvxpy +except ImportError: pass + try: - import cvxpy + import mosek + solver = 'MOSEK' # https://www.mosek.com/ except ImportError: - pass + warn("MOSEK not installed, falling back to CVXPY's defaults") + solver = None # passing this to solve() allows CVXPY to use whatever it deems best def iterative_velocity(x, dt, params, options=None, num_iterations=None, gamma=None, cg_maxiter=1000, scale='small'): @@ -56,62 +61,52 @@ def iterative_velocity(x, dt, params, options=None, num_iterations=None, gamma=N return x_hat, dxdt_hat -def __total_variation_regularized_derivative__(x, dt, N, gamma, solver='MOSEK'): - """Generalized total variation regularized derivatives - Use convex optimization (cvxpy) to solve for the Nth total variation regularized derivative. - Default solver is MOSEK: https://www.mosek.com/ +def _total_variation_regularized_derivative(x, dt, order, gamma, solver=solver): + """Generalized total variation regularized derivatives. Use convex optimization (cvxpy) to solve for a + total variation regularized derivative. Default solver is MOSEK: , which is not + always available. Fall back to CVXPY's default. - :param x: (np.array of floats, 1xN) time series to differentiate - :param dt: (float) time step - :param N: (int) 1, 2, or 3, the Nth derivative to regularize - :param gamma: (float) regularization parameter - :param solver: (string) Solver to use. Solver options include: 'MOSEK' and 'CVXOPT', - in testing, 'MOSEK' was the most robust. - :return: x_hat : estimated (smoothed) x - dxdt_hat : estimated derivative of x + :param np.array[float] x: data to differentiate + :param float dt: time step + :param int order: 1, 2, or 3, the derivative to regularize + :param float gamma: regularization parameter + :param str solver: Solver to use. Solver options include: 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS'. + In testing, 'MOSEK' was the most robust. + + :return: tuple[np.array, np.array] of\n + - **x_hat** -- estimated (smoothed) x + - **dxdt_hat** -- estimated derivative of x """ - # Normalize mean = np.mean(x) std = np.std(x) x = (x-mean)/std # Define the variables for the highest order derivative and the integration constants - var = cvxpy.Variable(len(x) + N) + 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 - derivatives = [var[N:]] - for i in range(N): - d = cvxpy.cumsum(derivatives[-1]) + var[i] - derivatives.append(d) - - # Compare the recursively integration position to the noisy position - sum_squared_error = cvxpy.sum_squares(derivatives[-1] - x) - - # Total variation regularization on the highest order derivative - r = cvxpy.sum(gamma*cvxpy.tv(derivatives[0])) + y = deriv_values + for i in range(order): + y = cvxpy.cumsum(y) + integration_constants[i] # Set up and solve the optimization problem - obj = cvxpy.Minimize(sum_squared_error + r) - prob = cvxpy.Problem(obj) + 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)) )) prob.solve(solver=solver) - # Recursively calculate the value of each derivative - final_derivative = var.value[N:] - derivative_values = [final_derivative] - for i in range(N): - d = np.cumsum(derivative_values[-1]) + var.value[i] - derivative_values.append(d) - for i, _ in enumerate(derivative_values): - derivative_values[i] = derivative_values[i]/(dt**(N-i)) - - # Extract the velocity and smoothed position - dxdt_hat = derivative_values[-2] - x_hat = derivative_values[-1] + # Recursively integrate the final derivative values to get back to the function and derivative values + y = deriv_values.value + for i in range(order-1): # stop one short to get the first derivative + y = np.cumsum(y) + integration_constants.value[i] + 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[0:-1] + dxdt_hat[1:])/2 + dxdt_hat = (dxdt_hat[0:-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 + 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 @@ -121,158 +116,116 @@ def __total_variation_regularized_derivative__(x, dt, N, gamma, solver='MOSEK'): return x_hat*std+mean, dxdt_hat*std -def velocity(x, dt, params, options=None): - """ - Use convex optimization (cvxpy) to solve for the velocity total variation regularized derivative. - Default solver is MOSEK: https://www.mosek.com/ - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: [gamma], where gamma (float) is the regularization parameter - or if 'iterate' in options: [gamma, num_iterations] - - :type params: list (float) or float - - :param options: {'solver': SOLVER} SOLVER options include: 'MOSEK' and 'CVXOPT', - in testing, 'MOSEK' was the most robust. - - :type options: dict {'solver': SOLVER}, optional - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x - - :rtype: tuple -> (np.array, np.array) - """ - - if options is None: - options = {'solver': 'MOSEK'} - - if isinstance(params, list): - gamma = params[0] - else: - gamma = params +def velocity(x, dt, params, options=None, gamma=None, solver=solver): + """Use convex optimization (cvxpy) to solve for the velocity total variation regularized derivative. - return __total_variation_regularized_derivative__(x, dt, 1, gamma, solver=options['solver']) - - -def acceleration(x, dt, params, options=None): - """ - Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative. - Default solver is MOSEK: https://www.mosek.com/ - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: [gamma], where gamma (float) is the regularization parameter - or if 'iterate' in options: [gamma, num_iterations] - - :type params: list (float) or float - - :param options: {'solver': SOLVER} SOLVER options include: 'MOSEK' and 'CVXOPT', - in testing, 'MOSEK' was the most robust. - - :type options: dict {'solver': SOLVER}, optional - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x + :param np.array[float] x: data to differentiate + :param float dt: time step size + :param params: (**deprecated**, prefer :code:`gamma`) + :param dict options: (**deprecated**, prefer :code:`solver`) a dictionary consisting of {'solver': (str)} + :param float gamma: the regularization parameter + :param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc. + In testing, 'MOSEK' was the most robust. Default is to use CVXPY's default. - :rtype: tuple -> (np.array, np.array) + :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 `gamma` " + + "and `solver` instead.", DeprecationWarning) + gamma = params[0] if isinstance(params, list) else params + if options != None: + if 'solver' in options: solver = options['solver'] + elif gamma == None: + raise ValueError("`gamma` must be given.") - if options is None: - options = {'solver': 'MOSEK'} - - if isinstance(params, list): - gamma = params[0] - else: - gamma = params + return _total_variation_regularized_derivative(x, dt, 1, gamma, solver=solver) - return __total_variation_regularized_derivative__(x, dt, 2, gamma, solver=options['solver']) +def acceleration(x, dt, params, options=None, gamma=None, solver=solver): + """Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative. + + :param np.array[float] x: data to differentiate + :param float dt: time step size + :param params: (**deprecated**, prefer :code:`gamma`) + :param dict options: (**deprecated**, prefer :code:`solver`) a dictionary consisting of {'solver': (str)} + :param float gamma: the regularization parameter + :param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc. + In testing, 'MOSEK' was the most robust. Default is to use CVXPY's default. -def jerk(x, dt, params, options=None): + :return: tuple[np.array, np.array] of\n + - **x_hat** -- estimated (smoothed) x + - **dxdt_hat** -- estimated derivative of x """ - Use convex optimization (cvxpy) to solve for the jerk total variation regularized derivative. - Default solver is MOSEK: https://www.mosek.com/ - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: [gamma], where gamma (float) is the regularization parameter - or if 'iterate' in options: [gamma, num_iterations] - - :type params: list (float) or float + 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 `gamma` " + + "and `solver` instead.", DeprecationWarning) + gamma = params[0] if isinstance(params, list) else params + if options != None: + if 'solver' in options: solver = options['solver'] + elif gamma == None: + raise ValueError("`gamma` must be given.") - :param options: {'solver': SOLVER} SOLVER options include: 'MOSEK' and 'CVXOPT', - in testing, 'MOSEK' was the most robust. + return _total_variation_regularized_derivative(x, dt, 2, gamma, solver=solver) - :type options: dict {'solver': SOLVER}, optional - :return: a tuple consisting of: +def jerk(x, dt, params, options=None, gamma=None, solver=solver): + """Use convex optimization (cvxpy) to solve for the jerk total variation regularized derivative. - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x + :param np.array[float] x: data to differentiate + :param float dt: time step size + :param params: (**deprecated**, prefer :code:`gamma`) + :param dict options: (**deprecated**, prefer :code:`solver`) a dictionary consisting of {'solver': (str)} + :param float gamma: the regularization parameter + :param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc. + In testing, 'MOSEK' was the most robust. Default is to use CVXPY's default. - :rtype: tuple -> (np.array, np.array) + :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 `gamma` " + + "and `solver` instead.", DeprecationWarning) + gamma = params[0] if isinstance(params, list) else params + if options != None: + if 'solver' in options: solver = options['solver'] + elif gamma == None: + raise ValueError("`gamma` must be given.") - if options is None: - options = {'solver': 'MOSEK'} - - if isinstance(params, list): - gamma = params[0] - else: - gamma = params - - return __total_variation_regularized_derivative__(x, dt, 3, gamma, solver=options['solver']) + return _total_variation_regularized_derivative(x, dt, 3, gamma, solver=solver) -def smooth_acceleration(x, dt, params, options=None): - """ - Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative. - And then apply a convolutional gaussian smoother to the resulting derivative to smooth out the peaks. +def smooth_acceleration(x, dt, params, options=None, gamma=None, window_size=None, solver=solver): + """Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative, + and then apply a convolutional gaussian smoother to the resulting derivative to smooth out the peaks. The end result is similar to the jerk method, but can be more time-efficient. - Default solver is MOSEK: https://www.mosek.com/ - - :param x: time series to differentiate - :type x: np.array of floats, 1xN - - :param dt: time step - :type dt: float - - :param params: list with values [gamma, window_size], where gamma (float) is the regularization parameter, window_size (int) is the window_size to use for the gaussian kernel - :type params: list -> [float, int] - - :param options: a dictionary indicating which SOLVER option to use, ie. 'MOSEK' or 'CVXOPT', in testing, 'MOSEK' was the most robust. - :type options: dict {'solver': SOLVER} - - :return: a tuple consisting of: - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x - :rtype: tuple -> (np.array, np.array) + :param np.array[float] x: data to differentiate + :param float dt: time step size + :param params: (**deprecated**, prefer :code:`gamma` and :code:`window_size`) + :param dict options: (**deprecated**, prefer :code:`solver`) a dictionary consisting of {'solver': (str)} + :param float gamma: the regularization parameter + :param int window_size: window size for gaussian kernel + :param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc. + In testing, 'MOSEK' was the most robust. Default is to use CVXPY's default. + :return: tuple[np.array, np.array] of\n + - **x_hat** -- estimated (smoothed) x + - **dxdt_hat** -- estimated derivative of x """ - if options is None: - options = {'solver': 'MOSEK'} + 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 `gamma` " + + "and `solver` instead.", DeprecationWarning) + gamma, window_size = params + if options != None: + if 'solver' in options: solver = options['solver'] + elif gamma == None or window_size == None: + raise ValueError("`gamma` and `window_size` must be given.") - gamma, window_size = params + _, dxdt_hat = _total_variation_regularized_derivative(x, dt, 2, gamma, solver=solver) - x_hat, dxdt_hat = acceleration(x, dt, [gamma], options=options) kernel = utility._gaussian_kernel(window_size) dxdt_hat = utility.convolutional_smoother(dxdt_hat, kernel, 1) @@ -283,48 +236,35 @@ def smooth_acceleration(x, dt, params, options=None): return x_hat, dxdt_hat -def jerk_sliding(x, dt, params, options=None): - """ - Use convex optimization (cvxpy) to solve for the jerk total variation regularized derivative. - Default solver is MOSEK: https://www.mosek.com/ - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: [gamma], where gamma (float) is the regularization parameter - or if 'iterate' in options: [gamma, num_iterations] - - :type params: list (float) or float - - :param options: {'solver': SOLVER} SOLVER options include: 'MOSEK' and 'CVXOPT', - in testing, 'MOSEK' was the most robust. +def jerk_sliding(x, dt, params, options=None, gamma=None, solver=solver): + """Use convex optimization (cvxpy) to solve for the jerk total variation regularized derivative. - :type options: dict {'solver': SOLVER}, optional - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x + :param np.array[float] x: data to differentiate + :param float dt: time step size + :param params: (**deprecated**, prefer :code:`gamma`) + :param dict options: (**deprecated**, prefer :code:`solver`) a dictionary consisting of {'solver': (str)} + :param float gamma: the regularization parameter + :param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc. + In testing, 'MOSEK' was the most robust. Default is to use CVXPY's default. - :rtype: tuple -> (np.array, np.array) + :return: tuple[np.array, np.array] of\n + - **x_hat** -- estimated (smoothed) x + - **dxdt_hat** -- estimated derivative of x """ - - if options is None: - options = {'solver': 'MOSEK'} - - if isinstance(params, list): - gamma = params[0] - else: - gamma = params + 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 `gamma` " + + "and `solver` instead.", DeprecationWarning) + gamma = params[0] if isinstance(params, list) else params + if options != None: + if 'solver' in options: solver = options['solver'] + elif gamma == None: + raise ValueError("`gamma` must be given.") window_size = 1000 stride = 200 if len(x) < window_size: - return jerk(x, dt, params, options=options) + return _total_variation_regularized_derivative(x, dt, 3, gamma, solver=solver) # slide the jerk final_xsmooth = [] @@ -337,8 +277,8 @@ def jerk_sliding(x, dt, params, options=None): try: while not last_loop: - xsmooth, xdot_hat = __total_variation_regularized_derivative__(x[first_idx:final_idx], dt, 3, - gamma, solver=options['solver']) + xsmooth, xdot_hat = _total_variation_regularized_derivative(x[first_idx:final_idx], dt, 3, + gamma, solver=solver) xsmooth = np.hstack(([0]*first_idx, xsmooth, [0]*(len(x)-final_idx))) final_xsmooth.append(xsmooth) @@ -379,5 +319,6 @@ def jerk_sliding(x, dt, params, options=None): return xsmooth, xdot_hat except ValueError: - print('Solver failed, returning finite difference instead') - return pynumdiff.utils.utility.finite_difference(x, dt) + warn('Solver failed, returning finite difference instead') + from pynumdiff.finite_difference import second_order + return second_order(x, dt) From 102ccd25821f1cf042dc84f7947d552ceacb5220 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Thu, 26 Jun 2025 16:15:37 -0700 Subject: [PATCH 4/5] excluding the test_diff_methods file --- pynumdiff/tests/test_diff_methods.py | 118 +++++------------- .../_total_variation_regularization.py | 2 +- 2 files changed, 35 insertions(+), 85 deletions(-) diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index 8c6e150..5e9796d 100644 --- 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 -def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) +iterated_first_order = lambda *args, **kwargs: first_order(*args, **kwargs) dt = 0.1 -t = np.arange(0, 3+dt, dt) # sample locations, including the endpoint +t = np.arange(0, 3, dt) # sample locations 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,7 +22,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) (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 @@ -30,109 +30,62 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) # 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, {}), # empty dictionary for the case of no parameters. no params -> no diff in new vs old - (iterated_first_order, {'num_iterations':5}), (iterated_first_order, [5], {'iterate':True}), - (second_order, {}), + (first_order, None), (iterated_first_order, {'num_iterations':5}), + (second_order, None), #(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]), - (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]) + (spectraldiff, {'high_freq_cutoff':0.1}), (spectraldiff, [0.1]) ] -diff_methods_and_params = [(velocity, [0.5], {'solver': 'CVXOPT'})] # 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), (-13, -14), (0, 0), (1, 1)], + [(-25, -25), (-14, -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: [[(-8, -9), (-11, -11), (0, -1), (0, 0)], - [(-6, -6), (-6, -7), (0, -1), (0, 0)], + iterated_first_order: [[(-7, -7), (-9, -10), (0, -1), (0, 0)], + [(-5, -5), (-5, -6), (0, -1), (0, 0)], [(-1, -1), (0, 0), (0, -1), (0, 0)], - [(0, 0), (1, 0), (0, 0), (1, 0)], + [(0, 0), (1, 1), (0, 0), (1, 1)], [(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), (-13, -13), (0, 0), (1, 1)], - [(-25, -25), (-13, -13), (0, 0), (1, 1)], + [(-25, -25), (-14, -14), (0, 0), (1, 1)], + [(-25, -25), (-13, -14), (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)], [(-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), (1, 1), (0, -1), (1, 1)], + [(0, 0), (2, 1), (0, 0), (2, 1)], [(0, 0), (3, 3), (0, 0), (3, 3)]], - savgoldiff: [[(-9, -10), (-13, -14), (0, -1), (0, 0)], - [(-9, -10), (-13, -13), (0, -1), (0, 0)], + savgoldiff: [[(-7, -8), (-13, -14), (0, -1), (0, 0)], + [(-7, -8), (-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: [[(-9, -10), (-14, -15), (-1, -1), (0, 0)], - [(0, 0), (1, 1), (0, 0), (1, 1)], - [(1, 0), (1, 1), (1, 1), (1, 1)], + spectraldiff: [[(-7, -8), (-14, -15), (-1, -1), (0, 0)], [(0, 0), (1, 1), (0, 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)], + [(1, 0), (1, 1), (1, 0), (1, 1)], [(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)]], - 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)]] + [(1, 1), (3, 3), (1, 1), (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 - # 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 + diff_method, params = diff_method_and_params # unpack 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 @@ -144,30 +97,27 @@ 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, 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) + # 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) 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) and len(diff_method_and_params) < 3) \ - else diff_method(x_noisy, dt, params, options) + else diff_method(x_noisy, dt, params) if isinstance(params, list) else diff_method(x_noisy, dt) # 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="") + #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)) - #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=", ") - # 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 {diff_method.__name__}") + #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)}") - 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 + if request.config.getoption("--plot"): # Get the plot flag from pytest configuration + fig, axes = request.config.plots[diff_method] 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/total_variation_regularization/_total_variation_regularization.py b/pynumdiff/total_variation_regularization/_total_variation_regularization.py index c649e44..1f22a9a 100644 --- a/pynumdiff/total_variation_regularization/_total_variation_regularization.py +++ b/pynumdiff/total_variation_regularization/_total_variation_regularization.py @@ -4,7 +4,7 @@ from pynumdiff.total_variation_regularization import _chartrand_tvregdiff from pynumdiff.utils import utility -try:import cvxpy +try: import cvxpy except ImportError: pass try: From d67270f20bb266931843588a71b14c9a04ff28f4 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Thu, 26 Jun 2025 16:21:15 -0700 Subject: [PATCH 5/5] pulled test_diff_methods from wrong branch before --- pynumdiff/tests/test_diff_methods.py | 111 +++++++++++++++++++-------- 1 file changed, 80 insertions(+), 31 deletions(-) diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index 5e9796d..def67bc 100644 --- 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,7 +22,7 @@ (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 @@ -30,12 +30,19 @@ # 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, {}), # empty dictionary for the case of no parameters. no params -> no diff in new vs old + (iterated_first_order, {'num_iterations':5}), (iterated_first_order, [5], {'iterate':True}), + (second_order, {}), #(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 @@ -43,49 +50,88 @@ # 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), (-11, -11), (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)], [(-14, -14), (-13, -13), (0, -1), (1, 1)], - [(-14, -15), (-13, -14), (0, -1), (1, 1)], + [(-14, -14), (-13, -13), (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}$")