From ab43990902b4d2652c597c20860cb86ee3bdd83d Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Mon, 23 Jun 2025 20:31:29 -0700 Subject: [PATCH 1/6] halfway through the docstrings, and then gotta add them to the new unittest structure. Should be easy. --- .../_smooth_finite_difference.py | 160 +++++++----------- 1 file changed, 57 insertions(+), 103 deletions(-) diff --git a/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py b/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py index 58a4cde..fcfb834 100644 --- a/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py +++ b/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py @@ -1,5 +1,6 @@ import numpy as np import scipy.signal +from warnings import warn # included code from pynumdiff.finite_difference import first_order as finite_difference @@ -9,140 +10,93 @@ ################################ # Smoothing finite differences # ################################ -def mediandiff(x, dt, params, options={}): - """ - Perform median smoothing using scipy.signal.medfilt - followed by first order finite difference - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: filter window size - :type params: list (int) or int - - :param options: an empty dictionary or a dictionary with 1 key value pair - - - 'iterate': whether to run multiple iterations of the smoother. Note: iterate does nothing for median smoother. - - :type options: dict {'iterate': (boolean)} +def mediandiff(x, dt, params=None, options={}, window_size=5, num_iterations=1): + """Perform median smoothing using scipy.signal.medfilt followed by first order finite difference - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x + :param np.array[float] x: array of time series to differentiate + :param float dt: time step size + :param list[int] params: (**deprecated**, prefer :code:`window_size` and :code:`num_iterations`) + :param dict options: (**deprecated**, prefer :code:`num_iterations`) an empty dictionary or {'iterate': (bool)} + :param int window_size: filter window size + :param int num_iterations: how many times to apply median smoothing - :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 'iterate' in options.keys() and options['iterate'] is True: - window_size, iterations = params - else: - iterations = 1 - if isinstance(params, list): - window_size = params[0] - else: - window_size = 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 `window_size` + and `num_iterations` instead.""", DeprecationWarning) + window_size = params[0] if isinstance(params, list) else params + if 'iterate' in options and options['iterate']: + num_iterations = params[1] if not window_size % 2: - window_size += 1 # assert window_size % 2 == 1 # is odd + window_size += 1 # make sure window_size is odd x_hat = x - for _ in range(iterations): + for _ in range(num_iterations): x_hat = scipy.signal.medfilt(x_hat, window_size) x_hat, dxdt_hat = finite_difference(x_hat, dt) return x_hat, dxdt_hat -def meandiff(x, dt, params, options={}): - """ - Perform mean smoothing by convolving mean kernel with x - followed by first order finite difference +def meandiff(x, dt, params=None, options={}, window_size=5, num_iterations=1): + """Perform mean smoothing by convolving mean kernel with x followed by first order finite difference :param np.ndarray[float] x: array of time series to differentiate :param float dt: time step size - :param params: [filter_window_size] or if 'iterate' in options: - [filter_window_size, num_iterations] - - :type params: list (int) - - :param options: an empty dictionary or a dictionary with 1 key value pair - - - 'iterate': whether to run multiple iterations of the smoother. Note: iterate does nothing for median smoother. - - :type options: dict {'iterate': (boolean)} - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x - + :param list[int] params: (**deprecated**, prefer :code:`window_size` and :code:`num_iterations`) + :code:`[window_size]` or, :code:`if 'iterate' in options`, :code:`[window_size, num_iterations]` + :param dict options: (**deprecated**, prefer :code:`num_iterations`) an empty dictionary or {'iterate': (bool)} + :param int window_size: filter window size + :param int num_iterations: how many times to apply mean smoothing - :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 'iterate' in options.keys() and options['iterate'] is True: - window_size, iterations = params - else: - iterations = 1 - if isinstance(params, list): - window_size = params[0] - else: - window_size = 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 `window_size` + and `num_iterations` instead.""", DeprecationWarning) + window_size = params[0] if isinstance(params, list) else params + if 'iterate' in options and options['iterate']: + num_iterations = params[1] kernel = utility._mean_kernel(window_size) - x_hat = utility.convolutional_smoother(x, kernel, iterations) + x_hat = utility.convolutional_smoother(x, kernel, num_iterations) x_hat, dxdt_hat = finite_difference(x_hat, dt) return x_hat, dxdt_hat -def gaussiandiff(x, dt, params, options={}): - """ - Perform gaussian smoothing by convolving gaussian kernel with x - followed by first order finite difference - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: [filter_window_size] or if 'iterate' in options: - [filter_window_size, num_iterations] - - :type params: list (int) - - :param options: an empty dictionary or a dictionary with 1 key value pair - - - 'iterate': whether to run multiple iterations of the smoother. Note: iterate does nothing for median smoother. - - :type options: dict {'iterate': (boolean)} - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x - +def gaussiandiff(x, dt, params=None, options={}, window_size=5, num_iterations=1): + """Perform gaussian smoothing by convolving gaussian kernel with x followed by first order finite difference - :rtype: tuple -> (np.array, np.array) + :param np.array[float] x: array of time series to differentiate + :param float dt: time step size + :param list[int] params: (**deprecated**, prefer :code:`window_size` and :code:`num_iterations`) + :code:`[window_size]` or, :code:`if 'iterate' in options`, :code:`[window_size, num_iterations]` + :param dict options: (**deprecated**, prefer :code:`num_iterations`) an empty dictionary or {'iterate': (bool)} + :param int window_size: filter window size + :param int num_iterations: how many times to apply gaussian smoothing + + :return: tuple[np.array, np.array] of\n + - **x_hat** -- estimated (smoothed) x + - **dxdt_hat** -- estimated derivative of x """ - if 'iterate' in options.keys() and options['iterate'] is True: - window_size, iterations = params - else: - iterations = 1 - if isinstance(params, list): - window_size = params[0] - else: - window_size = 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 `window_size` + and `num_iterations` instead.""", DeprecationWarning) + window_size = params[0] if isinstance(params, list) else params + if 'iterate' in options and options['iterate']: + num_iterations = params[1] kernel = utility._gaussian_kernel(window_size) - x_hat = utility.convolutional_smoother(x, kernel, iterations) + x_hat = utility.convolutional_smoother(x, kernel, num_iterations) x_hat, dxdt_hat = finite_difference(x_hat, dt) return x_hat, dxdt_hat From 581515880af66c6446bd9457a01f1713e8a4b12a Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Wed, 25 Jun 2025 11:48:59 -0700 Subject: [PATCH 2/6] finished #68 for smooth_finite_difference module --- .../_smooth_finite_difference.py | 187 +++++++----------- 1 file changed, 69 insertions(+), 118 deletions(-) diff --git a/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py b/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py index fcfb834..0d924e9 100644 --- a/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py +++ b/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py @@ -26,8 +26,8 @@ def mediandiff(x, dt, params=None, options={}, window_size=5, 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 `window_size` - and `num_iterations` instead.""", DeprecationWarning) + warn("`params` and `options` parameters will be removed in a future version. Use `window_size` " + + "and `num_iterations` instead.", DeprecationWarning) window_size = params[0] if isinstance(params, list) else params if 'iterate' in options and options['iterate']: num_iterations = params[1] @@ -60,8 +60,8 @@ def meandiff(x, dt, params=None, options={}, window_size=5, 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 `window_size` - and `num_iterations` instead.""", DeprecationWarning) + warn("`params` and `options` parameters will be removed in a future version. Use `window_size` " + + "and `num_iterations` instead.", DeprecationWarning) window_size = params[0] if isinstance(params, list) else params if 'iterate' in options and options['iterate']: num_iterations = params[1] @@ -89,8 +89,8 @@ def gaussiandiff(x, dt, params=None, options={}, window_size=5, 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 `window_size` - and `num_iterations` instead.""", DeprecationWarning) + warn("`params` and `options` parameters will be removed in a future version. Use `window_size` " + + "and `num_iterations` instead.", DeprecationWarning) window_size = params[0] if isinstance(params, list) else params if 'iterate' in options and options['iterate']: num_iterations = params[1] @@ -102,99 +102,65 @@ def gaussiandiff(x, dt, params=None, options={}, window_size=5, num_iterations=1 return x_hat, dxdt_hat -def friedrichsdiff(x, dt, params, options={}): - """ - Perform friedrichs smoothing by convolving friedrichs kernel with x - followed by first order finite difference - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: [filter_window_size] or if 'iterate' in options: - [filter_window_size, num_iterations] - - :type params: list (int) - - :param options: an empty dictionary or a dictionary with 1 key value pair - - - 'iterate': whether to run multiple iterations of the smoother. Note: iterate does nothing for median smoother. - - :type options: dict {'iterate': (boolean)} - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x +def friedrichsdiff(x, dt, params=None, options={}, window_size=5, num_iterations=1): + """Perform friedrichs smoothing by convolving friedrichs kernel with x 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[int] params: (**deprecated**, prefer :code:`window_size` and :code:`num_iterations`) + :code:`[window_size]` or, :code:`if 'iterate' in options`, :code:`[window_size, num_iterations]` + :param dict options: (**deprecated**, prefer :code:`num_iterations`) an empty dictionary or {'iterate': (bool)} + :param int window_size: filter window size + :param int num_iterations: how many times to apply smoothing - :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 'iterate' in options.keys() and options['iterate'] is True: - window_size, iterations = params - else: - iterations = 1 - if isinstance(params, list): - window_size = params[0] - else: - window_size = 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 `window_size` " + + "and `num_iterations` instead.", DeprecationWarning) + window_size = params[0] if isinstance(params, list) else params + if 'iterate' in options and options['iterate']: + num_iterations = params[1] kernel = utility._friedrichs_kernel(window_size) - x_hat = utility.convolutional_smoother(x, kernel, iterations) + x_hat = utility.convolutional_smoother(x, kernel, num_iterations) x_hat, dxdt_hat = finite_difference(x_hat, dt) return x_hat, dxdt_hat -def butterdiff(x, dt, params, options={}): - """ - Perform butterworth smoothing on x with scipy.signal.filtfilt - followed by first order finite difference - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: [n, wn], n = order of the filter; wn = Cutoff frequency. - For a discrete timeseries, the value is normalized to the range 0-1, - where 1 is the Nyquist frequency. - - :type params: list (int) - - :param options: an empty dictionary or a dictionary with 2 key value pair - - - 'iterate': whether to run multiple iterations of the smoother. Note: iterate does nothing for median smoother. - - 'padmethod': "pad" or "gust", see scipy.signal.filtfilt - - :type options: dict {'iterate': (boolean), 'padmethod': string} - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x +def butterdiff(x, dt, params=None, options={}, filter_order=2, cutoff_freq=0.5, num_iterations=1): + """Perform butterworth smoothing on x with scipy.signal.filtfilt 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[int] params: (**deprecated**, prefer :code:`filter_order`, :code:`cutoff_freq`, + and :code:`num_iterations`) + :param dict options: (**deprecated**, prefer :code:`num_iterations`) an empty dictionary or {'iterate': (bool)} + :param int filter_order: order of the filter + :param float cutoff_freq: cutoff frequency :math:`\\in [0, 1]`. For a discrete timeseries, + the value is normalized to the range 0-1, where 1 is the Nyquist frequency. + :param int num_iterations: how many times to apply smoothing - :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 'iterate' in options.keys() and options['iterate'] is True: - n, wn, iterations = params - else: - iterations = 1 - n, wn = 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 `filter_order`, " + + "`cutoff_freq`, and `num_iterations` instead.", DeprecationWarning) + filter_order, cutoff_freq = params[0:2] + if 'iterate' in options and options['iterate']: + num_iterations = params[2] - b, a = scipy.signal.butter(n, wn) + b, a = scipy.signal.butter(filter_order, cutoff_freq) x_hat = x - for _ in range(iterations): - if len(x) < 9: - x_hat = scipy.signal.filtfilt(b, a, x_hat, method="pad", padlen=len(x)-1) - else: - x_hat = scipy.signal.filtfilt(b, a, x_hat, method="pad") + padlen = len(x)-1 if len(x) < 9 else None + for _ in range(num_iterations): + x_hat = scipy.signal.filtfilt(b, a, x_hat, method="pad", padlen=padlen) x_hat, dxdt_hat = finite_difference(x_hat, dt) @@ -204,48 +170,33 @@ def butterdiff(x, dt, params, options={}): return x_hat, dxdt_hat -def splinediff(x, dt, params, options={}): - """ - Perform spline smoothing on x with scipy.interpolate.UnivariateSpline - followed by first order finite difference - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: [k, s], k: Order of the spline. A kth order spline can be differentiated k times. - s: Positive smoothing factor used to choose the number of knots. - Number of knots will be increased until the smoothing condition is satisfied: - sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s - - :type params: list (int) - - :param options: an empty dictionary or a dictionary with 1 key value pair - - - 'iterate': whether to run multiple iterations of the smoother. Note: iterate does nothing for median smoother. - - :type options: dict {'iterate': (boolean)} - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x +def splinediff(x, dt, params=None, options={}, k=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 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 - :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 'iterate' in options.keys() and options['iterate'] is True: - k, s, iterations = params - else: - iterations = 1 - k, s = 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 `k`, `s`, and " + + "`num_iterations` instead.", DeprecationWarning) + k, 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(iterations): + for _ in range(num_iterations): spline = scipy.interpolate.UnivariateSpline(t, x_hat, k=k, s=s) x_hat = spline(t) From 7e6519e6871338127d825b0a59841caba3304a95 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Wed, 25 Jun 2025 13:32:39 -0700 Subject: [PATCH 3/6] moved smooth_finite_difference_tests and in the process found a way to enable options dictionaries in my old-style invocations --- .../_smooth_finite_difference.py | 10 +- pynumdiff/tests/test_diff_methods.py | 113 +++++++++++++----- .../tests/test_smooth_finite_difference.py | 78 ------------ 3 files changed, 86 insertions(+), 115 deletions(-) delete mode 100644 pynumdiff/tests/test_smooth_finite_difference.py diff --git a/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py b/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py index 0d924e9..f66c924 100644 --- 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,9 +187,9 @@ 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] @@ -197,7 +197,7 @@ def splinediff(x, dt, params=None, options={}, k=3, s=None, num_iterations=1): 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 index 5e9796d..d9e945d 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, {}), + (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 @@ -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), (-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 deleted file mode 100644 index 5429990..0000000 --- a/pynumdiff/tests/test_smooth_finite_difference.py +++ /dev/null @@ -1,78 +0,0 @@ -""" -Unit tests for smoothing + finite difference methods -""" -# pylint: skip-file -import numpy as np -from pynumdiff.smooth_finite_difference import mediandiff, meandiff, gaussiandiff, \ - friedrichsdiff, butterdiff, splinediff - - -x = np.array([1., 4., 9., 3., 20., - 8., 16., 2., 15., 10., - 15., 3., 5., 7., 4.]) -dt = 0.01 - - -def test_mediandiff(): - params = [3, 2] - x_hat, dxdt_hat = mediandiff(x, dt, params, options={'iterate': True}) - x_smooth = np.array([1., 4., 4., 8., 9., 8., 15., 10., 15., 10., 10., 5., 5., 5., 4.]) - dxdt = np.array([300., 150., 200., 250., 0., 300., 100., 0., - 0., -250., -250., -250., 0., -50., -100.]) - np.testing.assert_array_equal(x_smooth, x_hat) - np.testing.assert_array_equal(dxdt, dxdt_hat) - -def test_meandiff(): - params = [3, 2] - x_hat, dxdt_hat = meandiff(x, dt, params, options={'iterate': True}) - x_smooth = np.array([2.889, 4., 6.889, 8.778, 11.889, 11.222, 11.444, 9.556, - 11.111, 10.556, 10.111, 7.333, 6., 5.111, 5.111]) - dxdt = np.array([111.111, 200., 238.889, 250., 122.222, -22.222, - -83.333, -16.667, 50., -50., -161.111, -205.556, - -111.111, -44.444, 0.]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) - -def test_gaussiandiff(): - params = [5] - x_hat, dxdt_hat = gaussiandiff(x, dt, params, options={'iterate': False}) - x_smooth = np.array([1.805, 4.377, 6.66, 8.066, 13.508, 12.177, 11.278, 8.044, - 11.116, 11.955, 11.178, 6.187, 5.127, 5.819, 4.706]) - dxdt = np.array([257.235, 242.77, 184.438, 342.42, 205.553, -111.535, - -206.61, -8.093, 195.509, 3.089, -288.392, -302.545, - -18.409, -21.032, -111.263]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) - -def test_friedrichsdiff(): - params = [5] - x_hat, dxdt_hat = friedrichsdiff(x, dt, params, options={'iterate': False}) - x_smooth = np.array([1.884, 4.589, 5.759, 9.776, 11.456, 13.892, 9.519, 9.954, - 9.697, 12.946, 9.992, 7.124, 5., 5.527, 4.884]) - dxdt = np.array([270.539, 193.776, 259.335, 284.855, 205.809, -96.888, - -196.888, 8.921, 149.586, 14.73, -291.079, -249.586, - -79.875, -5.809, -64.316]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) - -def test_butterdiff(): - params = [3, 0.074] - x_hat, dxdt_hat = butterdiff(x, dt, params, options={'iterate': False}) - x_smooth = np.array([3.445, 4.753, 5.997, 7.131, 8.114, 8.914, 9.51, 9.891, - 10.058, 10.02, 9.798, 9.42, 8.919, 8.332, 7.699]) - dxdt = np.array([130.832, 127.617, 118.881, 105.827, 89.169, 69.833, 48.871, - 27.381, 6.431, -12.992, -30.023, -43.972, -54.368, -60.98, - -63.326]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) - -def test_splinediff(): - params = [5, 2] - x_hat, dxdt_hat = splinediff(x, dt, params, options={'iterate': False}) - x_smooth = np.array([0.995, 4.035, 8.874, 3.279, 19.555, 8.564, 15.386, 2.603, - 14.455, 10.45, 14.674, 3.193, 4.916, 7.023, 3.997]) - dxdt = np.array([303.996, 393.932, -37.815, 534.063, 264.225, -208.442, - -298.051, -46.561, 392.365, 10.93, -362.858, -487.87, - 191.508, -45.968, -302.579]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) From 4c734c46c96bce1f275a97546e4ef07571673e7b Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Wed, 25 Jun 2025 13:41:11 -0700 Subject: [PATCH 4/6] adjusting error bounds to get this to work in the CI --- pynumdiff/tests/test_diff_methods.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index d9e945d..36acf76 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -55,7 +55,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) [(-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), (-25, -25), (0, -1), (0, 0)], + iterated_first_order: [[(-8, -9), (-11, -25), (0, -1), (0, 0)], [(-6, -6), (-6, -7), (0, -1), (0, 0)], [(-1, -1), (0, 0), (0, -1), (0, 0)], [(0, 0), (1, 0), (0, 0), (1, 0)], @@ -68,7 +68,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) [(-25, -25), (1, 1), (0, 0), (1, 1)], [(-25, -25), (3, 3), (0, 0), (3, 3)]], #lineardiff: [TBD when #91 is solved], - polydiff: [[(-15, -15), (-14, -14), (0, -1), (1, 1)], + 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)], [(-2, -2), (0, 0), (0, -1), (1, 1)], From 8409c97a024f35eb7effebd21a79398097bc0745 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Wed, 25 Jun 2025 13:44:24 -0700 Subject: [PATCH 5/6] adjusting error bounds more. Last failure hid the present failure, though I suspected it would fail. This should work now --- pynumdiff/tests/test_diff_methods.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index 36acf76..f66adf0 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -55,7 +55,7 @@ def iterated_first_order(*args, **kwargs): return first_order(*args, **kwargs) [(-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, -25), (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, 0), (0, 0), (1, 0)], From 96f4631dbdb0128557bd620c3adcfc21be54e086 Mon Sep 17 00:00:00 2001 From: pavelkomarov Date: Wed, 25 Jun 2025 14:00:21 -0700 Subject: [PATCH 6/6] added test for iterated first order, called by old method signature --- pynumdiff/tests/test_diff_methods.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pynumdiff/tests/test_diff_methods.py b/pynumdiff/tests/test_diff_methods.py index f66adf0..def67bc 100644 --- a/pynumdiff/tests/test_diff_methods.py +++ b/pynumdiff/tests/test_diff_methods.py @@ -30,9 +30,9 @@ 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, {}), - (iterated_first_order, {'num_iterations':5}), - (second_order, {}), # empty dictionary for the case of no parameters + (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]),