23028417751

Committed 12 Mar 2026 11:16PM UTC coverage: 91.932% (+0.07%) from 91.863%

Build # 23028417751

Build Type

Pull #188

github

Committed by

web-flow

Commit Message

Merge 882a9023d into dff1f9c95

Pull Request Pull Request #188: I added multidimensionality functions to the tvrdiff method by callin…

Coverage Stats

26 of 26 new or added lines in 1 file covered. (100.0%)

4 existing lines in 1 file now uncovered.

866 of 942 relevant lines covered (91.93%)

0.92 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

93.9

/pynumdiff/total_variation_regularization.py

"""This module implements some common total variation regularization methods."""
from warnings import warn
import numpy as np
from scipy.stats import median_abs_deviation
try: import cvxpy
except ImportError: pass

from pynumdiff.utils import _chartrand_tvregdiff, utility


def iterative_velocity(x, dt, params=None, 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: data to differentiate
    :param float dt: 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 and 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: - **x_hat** (np.array) -- estimated (smoothed) x
             - **dxdt_hat** (np.array) -- estimated derivative of x
    """
    if params is not 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 is not None:
            if 'cg_maxiter' in options: cg_maxiter = options['cg_maxiter']
            if 'scale' in options: scale = options['scale']
    elif num_iterations is None or gamma is 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)
    x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt)
    x0 = utility.estimate_integration_constant(x, x_hat)
    x_hat = x_hat + x0

    return x_hat, dxdt_hat

def tvrdiff(x, dt, order, gamma, huberM=float('inf'), solver=None, axis=0):
    """Generalized total variation regularized derivatives. Use convex optimization (cvxpy) to solve for a
    total variation regularized derivative. Other convex-solver-based methods in this module call this function.

    :param np.array[float] x: data to differentiate
    :param float dt: step size
    :param int order: 1, 2, or 3, the derivative to regularize
    :param float gamma: regularization parameter
    :param float huberM: Huber loss parameter, in units of scaled median absolute deviation of input data.
                    :math:`M = \\infty` reduces to :math:`\\ell_2` loss squared on first, fidelity cost term, and
                    :math:`M = 0` reduces to :math:`\\ell_1` loss, which seeks sparse residuals.
    :param str solver: Solver to use. Solver options include: 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS'.
                    If not given, fall back to CVXPY's default.
    :param int axis: data dimension along which differentiation is performed

    :return: - **x_hat** (np.array) -- estimated (smoothed) x
             - **dxdt_hat** (np.array) -- estimated derivative of x
    """
    x = np.asarray(x, dtype=float)

    # had to loop in python :(
    x_hat = np.empty_like(x)
    dxdt_hat = np.empty_like(x)
    for vec_idx in np.ndindex(x.shape[:axis] + x.shape[axis+1:]):
        s = vec_idx[:axis] + (slice(None),) + vec_idx[axis:]
        xv = x[s]

        # Normalize for numerical consistency with convex solver
        mu = np.mean(xv)
        sigma = median_abs_deviation(xv, scale='normal') # robust alternative to std()
        if sigma == 0: sigma = 1 # safety guard
        y = (xv-mu)/sigma

        # Define the variables for the highest order derivative and the integration constants
        deriv_values = cvxpy.Variable(len(y)) # 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. This is a first-
        # order scheme, but it's very fast and tends to do not markedly worse than 2nd order. See #116
        # I also tried a trapezoidal integration rule here, and it works no better. See #116 too.
        hx = deriv_values # variables are integrated to produce the signal estimate variables, \hat{x} in the math
        for i in range(order):
            hx = cvxpy.cumsum(hx) + integration_constants[i] # cumsum is like integration assuming dt = 1

        # Compare the recursively integrated position to the noisy position. \ell_2 doesn't get scaled by 1/2 here,
        # so cvxpy's doubled Huber is already the right scale, and \ell_1 should be scaled by 2\sqrt{2} to match.
        fidelity_cost = cvxpy.sum_squares(y - hx) if huberM == float('inf') \
                else np.sqrt(8)*cvxpy.norm(y - hx, 1) if huberM == 0 \
                else utility.huber_const(huberM)*cvxpy.sum(cvxpy.huber(y - hx, huberM)) # data is already scaled, so M rather than M*sigma
        # Set up and solve the optimization problem
        prob = cvxpy.Problem(cvxpy.Minimize(fidelity_cost + gamma*cvxpy.sum(cvxpy.tv(deriv_values)) ))
        prob.solve(solver=solver)

        # Recursively integrate the final derivative values to get back to the function and derivative values
        v = deriv_values.value
        for i in range(order-1): # stop one short to get the first derivative
            v = np.cumsum(v) + integration_constants.value[i]
        dxdt_v = v/dt # v only holds the dx values; to get deriv scale by dt
        x_v = np.cumsum(v) + integration_constants.value[order-1] # smoothed data

        # Due to the first-order nature of the derivative, it has a slight lag. Average together every two values
        # to better center the answer. But this leaves us one-short, so devise a good last value.
        dxdt_v = (dxdt_v[:-1] + dxdt_v[1:])/2
        dxdt_v = np.hstack((dxdt_v, 2*dxdt_v[-1] - dxdt_v[-2])) # last value = penultimate value [-1] + diff between [-1] and [-2]

        x_hat[s] = x_v*sigma+mu
        dxdt_hat[s] = dxdt_v*sigma # derivative is linear, so scale derivative by scatter

    return x_hat, dxdt_hat


def velocity(x, dt, params=None, options=None, gamma=None, solver=None):
    """Use convex optimization (cvxpy) to solve for the velocity total variation regularized derivative.\n
    **Deprecated**, prefer :code:`tvrdiff` with order 1 instead.

    :param np.array[float] x: data to differentiate
    :param float dt: 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.
                If not given, fall back to CVXPY's default.

    :return: - **x_hat** (np.array) -- estimated (smoothed) x
             - **dxdt_hat** (np.array) -- estimated derivative of x
    """
    if params is not 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 is not None:
            if 'solver' in options: solver = options['solver']
    elif gamma is None:
        raise ValueError("`gamma` must be given.")

    warn("`velocity` is deprecated. Call `tvrdiff` with order 1 instead.", DeprecationWarning)
    return tvrdiff(x, dt, 1, gamma, solver=solver)


def acceleration(x, dt, params=None, options=None, gamma=None, solver=None):
    """Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative.\n
    **Deprecated**, prefer :code:`tvrdiff` with order 2 instead.
    
    :param np.array[float] x: data to differentiate
    :param float dt: 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. If not given, fall back to CVXPY's default.

    :return: - **x_hat** (np.array) -- estimated (smoothed) x
             - **dxdt_hat** (np.array) -- estimated derivative of x
    """
    if params is not 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 is not None:
            if 'solver' in options: solver = options['solver']
    elif gamma is None:
        raise ValueError("`gamma` must be given.")

    warn("`acceleration` is deprecated. Call `tvrdiff` with order 2 instead.", DeprecationWarning)
    return tvrdiff(x, dt, 2, gamma, solver=solver)


def jerk(x, dt, params=None, options=None, gamma=None, solver=None):
    """Use convex optimization (cvxpy) to solve for the jerk total variation regularized derivative.\n
    **Deprecated**, prefer :code:`tvrdiff` with order 3 instead.

    :param np.array[float] x: data to differentiate
    :param float dt: 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. If not given, fall back to CVXPY's default.

    :return: - **x_hat** (np.array) -- estimated (smoothed) x
             - **dxdt_hat** (np.array) -- estimated derivative of x
    """
    if params is not 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 is not None:
            if 'solver' in options: solver = options['solver']
    elif gamma is None:
        raise ValueError("`gamma` must be given.")

    warn("`jerk` is deprecated. Call `tvrdiff` with order 3 instead.", DeprecationWarning)
    return tvrdiff(x, dt, 3, gamma, solver=solver)


def smooth_acceleration(x, dt, params=None, options=None, gamma=None, window_size=None, solver=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.
    The end result is similar to the jerk method, but can be more time-efficient.

    :param np.array[float] x: data to differentiate
    :param float dt: 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. If not given, fall back to CVXPY's default.

    :return: - **x_hat** (np.array) -- estimated (smoothed) x
             - **dxdt_hat** (np.array) -- estimated derivative of x
    """
    if params is not 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 is not None:
            if 'solver' in options: solver = options['solver']
    elif gamma is None or window_size is None:
        raise ValueError("`gamma` and `window_size` must be given.")

    _, dxdt_hat = tvrdiff(x, dt, 2, gamma, solver=solver)

    kernel = utility.gaussian_kernel(window_size)
    dxdt_hat = utility.convolutional_smoother(dxdt_hat, kernel, 1)

    x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt)
    x0 = utility.estimate_integration_constant(x, x_hat)
    x_hat = x_hat + x0

    return x_hat, dxdt_hat

1	"""This module implements some common total variation regularization methods."""
2	from warnings import warn	1✔
3	import numpy as np	1✔
4	from scipy.stats import median_abs_deviation	1✔
5	try: import cvxpy	1✔
6	except ImportError: pass
7
8	from pynumdiff.utils import _chartrand_tvregdiff, utility	1✔
9
10
11	def iterative_velocity(x, dt, params=None, options=None, num_iterations=None, gamma=None, cg_maxiter=1000, scale='small'):	1✔
12	"""Use an iterative solver to find the total variation regularized 1st derivative. See
13	_chartrand_tvregdiff.py for details, author info, and license. Methods described in:
14	Rick Chartrand, "Numerical differentiation of noisy, nonsmooth data," ISRN Applied Mathematics,
15	Vol. 2011, Article ID 164564, 2011. Original code at https://sites.google.com/site/dnartrahckcir/home/tvdiff-code
16
17	:param np.array[float] x: data to differentiate
18	:param float dt: step size
19	:param list params: (deprecated, prefer :code:`num_iterations` and :code:`gamma`)
20	:param dict options: (deprecated, prefer :code:`cg_maxiter` and :code:`scale`)
21	a dictionary consisting of {'cg_maxiter': (int), 'scale': (str)}
22	:param int num_iterations: number of iterations to run the solver. More iterations results in
23	blockier derivatives, which approach the convex result
24	:param float gamma: regularization parameter
25	:param int cg_maxiter: Max number of iterations to use in :code:`scipy.sparse.linalg.cg`. Default
26	:code:`None` results in maxiter = len(x). This works well in our test examples.
27	:param str scale: This method has two different numerical options. From :code:`_chartrand_tvregdiff.py`:
28	:code:`'large'` or :code:`'small'` (case insensitive). Default is :code:`'small'`. :code:`'small'`
29	has somewhat better boundary behavior, but becomes unwieldly for data larger than 1000 entries or so.
30	:code:`'large'` has simpler numerics and is more efficient for large-scale problems. :code:`'large'`
31	is more readily modified for higher-order derivatives, since the implicit differentiation matrix is square.
32
33	:return: - x_hat (np.array) -- estimated (smoothed) x
34	- dxdt_hat (np.array) -- estimated derivative of x
35	"""
36	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
37	warn("`params` and `options` parameters will be removed in a future version. Use `num_iterations`, " +	1✔
38	"`gamma`, `cg_maxiter`, and `scale` instead.", DeprecationWarning)
39	num_iterations, gamma = params	1✔
40	if options is not None:
41	if 'cg_maxiter' in options: cg_maxiter = options['cg_maxiter']
42	if 'scale' in options: scale = options['scale']
43	elif num_iterations is None or gamma is None:	1✔
44	raise ValueError("`num_iterations` and `gamma` must be given.")	×
45
46	dxdt_hat = _chartrand_tvregdiff.TVRegDiff(x, num_iterations, gamma, dx=dt,	1✔
47	maxit=cg_maxiter, scale=scale,
48	ep=1e-6, u0=None, plotflag=False)
49	x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt)	1✔
50	x0 = utility.estimate_integration_constant(x, x_hat)	1✔
51	x_hat = x_hat + x0	1✔
52
53	return x_hat, dxdt_hat	1✔
54
55	def tvrdiff(x, dt, order, gamma, huberM=float('inf'), solver=None, axis=0):	1✔
56	"""Generalized total variation regularized derivatives. Use convex optimization (cvxpy) to solve for a
57	total variation regularized derivative. Other convex-solver-based methods in this module call this function.
58
59	:param np.array[float] x: data to differentiate
60	:param float dt: step size
61	:param int order: 1, 2, or 3, the derivative to regularize
62	:param float gamma: regularization parameter
63	:param float huberM: Huber loss parameter, in units of scaled median absolute deviation of input data.
64	:math:`M = \\infty` reduces to :math:`\\ell_2` loss squared on first, fidelity cost term, and
65	:math:`M = 0` reduces to :math:`\\ell_1` loss, which seeks sparse residuals.
66	:param str solver: Solver to use. Solver options include: 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS'.
67	If not given, fall back to CVXPY's default.
68	:param int axis: data dimension along which differentiation is performed
69
70	:return: - x_hat (np.array) -- estimated (smoothed) x
71	- dxdt_hat (np.array) -- estimated derivative of x
72	"""
73	x = np.asarray(x, dtype=float)	1✔
74
75	# had to loop in python :(
76	x_hat = np.empty_like(x)	1✔
77	dxdt_hat = np.empty_like(x)	1✔
78	for vec_idx in np.ndindex(x.shape[:axis] + x.shape[axis+1:]):	1✔
79	s = vec_idx[:axis] + (slice(None),) + vec_idx[axis:]	1✔
80	xv = x[s]	1✔
81
82	# Normalize for numerical consistency with convex solver
83	mu = np.mean(xv)	1✔
84	sigma = median_abs_deviation(xv, scale='normal') # robust alternative to std()	1✔
85	if sigma == 0: sigma = 1 # safety guard	1✔
86	y = (xv-mu)/sigma	1✔
87
88	# Define the variables for the highest order derivative and the integration constants
89	deriv_values = cvxpy.Variable(len(y)) # values of the order^th derivative, in which we're penalizing variation	1✔
90	integration_constants = cvxpy.Variable(order) # constants of integration that help get us back to x	1✔
91
92	# Recursively integrate the highest order derivative to get back to the position. This is a first-
93	# order scheme, but it's very fast and tends to do not markedly worse than 2nd order. See #116
94	# I also tried a trapezoidal integration rule here, and it works no better. See #116 too.
95	hx = deriv_values # variables are integrated to produce the signal estimate variables, \hat{x} in the math	1✔
96	for i in range(order):	1✔
97	hx = cvxpy.cumsum(hx) + integration_constants[i] # cumsum is like integration assuming dt = 1	1✔
98
99	# Compare the recursively integrated position to the noisy position. \ell_2 doesn't get scaled by 1/2 here,
100	# so cvxpy's doubled Huber is already the right scale, and \ell_1 should be scaled by 2\sqrt{2} to match.
101	fidelity_cost = cvxpy.sum_squares(y - hx) if huberM == float('inf') \	1✔
102	else np.sqrt(8)*cvxpy.norm(y - hx, 1) if huberM == 0 \
103	else utility.huber_const(huberM)cvxpy.sum(cvxpy.huber(y - hx, huberM)) # data is already scaled, so M rather than Msigma
104	# Set up and solve the optimization problem
105	prob = cvxpy.Problem(cvxpy.Minimize(fidelity_cost + gamma*cvxpy.sum(cvxpy.tv(deriv_values)) ))	1✔
106	prob.solve(solver=solver)	1✔
107
108	# Recursively integrate the final derivative values to get back to the function and derivative values
109	v = deriv_values.value	1✔
110	for i in range(order-1): # stop one short to get the first derivative	1✔
111	v = np.cumsum(v) + integration_constants.value[i]	1✔
112	dxdt_v = v/dt # v only holds the dx values; to get deriv scale by dt	1✔
113	x_v = np.cumsum(v) + integration_constants.value[order-1] # smoothed data	1✔
114
115	# Due to the first-order nature of the derivative, it has a slight lag. Average together every two values
116	# to better center the answer. But this leaves us one-short, so devise a good last value.
117	dxdt_v = (dxdt_v[:-1] + dxdt_v[1:])/2	1✔
118	dxdt_v = np.hstack((dxdt_v, 2*dxdt_v[-1] - dxdt_v[-2])) # last value = penultimate value [-1] + diff between [-1] and [-2]	1✔
119
120	x_hat[s] = x_v*sigma+mu	1✔
121	dxdt_hat[s] = dxdt_v*sigma # derivative is linear, so scale derivative by scatter	1✔
122
123	return x_hat, dxdt_hat	1✔
124
125
126	def velocity(x, dt, params=None, options=None, gamma=None, solver=None):	1✔
127	"""Use convex optimization (cvxpy) to solve for the velocity total variation regularized derivative.\n
128	Deprecated, prefer :code:`tvrdiff` with order 1 instead.
129
130	:param np.array[float] x: data to differentiate
131	:param float dt: step size
132	:param params: (deprecated, prefer :code:`gamma`)
133	:param dict options: (deprecated, prefer :code:`solver`) a dictionary consisting of {'solver': (str)}
134	:param float gamma: the regularization parameter
135	:param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc.
136	If not given, fall back to CVXPY's default.
137
138	:return: - x_hat (np.array) -- estimated (smoothed) x
139	- dxdt_hat (np.array) -- estimated derivative of x
140	"""
141	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
142	warn("`params` and `options` parameters will be removed in a future version. Use `gamma` " +	1✔
143	"and `solver` instead.", DeprecationWarning)
144	gamma = params[0] if isinstance(params, list) else params	1✔
145	if options is not None:
146	if 'solver' in options: solver = options['solver']
147	elif gamma is None:	1✔
UNCOV 148	raise ValueError("`gamma` must be given.")	×
149
150	warn("`velocity` is deprecated. Call `tvrdiff` with order 1 instead.", DeprecationWarning)	1✔
151	return tvrdiff(x, dt, 1, gamma, solver=solver)	1✔
152
153
154	def acceleration(x, dt, params=None, options=None, gamma=None, solver=None):	1✔
155	"""Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative.\n
156	Deprecated, prefer :code:`tvrdiff` with order 2 instead.
157
158	:param np.array[float] x: data to differentiate
159	:param float dt: step size
160	:param params: (deprecated, prefer :code:`gamma`)
161	:param dict options: (deprecated, prefer :code:`solver`) a dictionary consisting of {'solver': (str)}
162	:param float gamma: the regularization parameter
163	:param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc.
164	In testing, 'MOSEK' was the most robust. If not given, fall back to CVXPY's default.
165
166	:return: - x_hat (np.array) -- estimated (smoothed) x
167	- dxdt_hat (np.array) -- estimated derivative of x
168	"""
169	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
170	warn("`params` and `options` parameters will be removed in a future version. Use `gamma` " +	1✔
171	"and `solver` instead.", DeprecationWarning)
172	gamma = params[0] if isinstance(params, list) else params	1✔
173	if options is not None:
174	if 'solver' in options: solver = options['solver']
175	elif gamma is None:	1✔
UNCOV 176	raise ValueError("`gamma` must be given.")	×
177
178	warn("`acceleration` is deprecated. Call `tvrdiff` with order 2 instead.", DeprecationWarning)	1✔
179	return tvrdiff(x, dt, 2, gamma, solver=solver)	1✔
180
181
182	def jerk(x, dt, params=None, options=None, gamma=None, solver=None):	1✔
183	"""Use convex optimization (cvxpy) to solve for the jerk total variation regularized derivative.\n
184	Deprecated, prefer :code:`tvrdiff` with order 3 instead.
185
186	:param np.array[float] x: data to differentiate
187	:param float dt: step size
188	:param params: (deprecated, prefer :code:`gamma`)
189	:param dict options: (deprecated, prefer :code:`solver`) a dictionary consisting of {'solver': (str)}
190	:param float gamma: the regularization parameter
191	:param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc.
192	In testing, 'MOSEK' was the most robust. If not given, fall back to CVXPY's default.
193
194	:return: - x_hat (np.array) -- estimated (smoothed) x
195	- dxdt_hat (np.array) -- estimated derivative of x
196	"""
197	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
198	warn("`params` and `options` parameters will be removed in a future version. Use `gamma` " +	1✔
199	"and `solver` instead.", DeprecationWarning)
200	gamma = params[0] if isinstance(params, list) else params	1✔
201	if options is not None:
202	if 'solver' in options: solver = options['solver']
203	elif gamma is None:	1✔
UNCOV 204	raise ValueError("`gamma` must be given.")	×
205
206	warn("`jerk` is deprecated. Call `tvrdiff` with order 3 instead.", DeprecationWarning)	1✔
207	return tvrdiff(x, dt, 3, gamma, solver=solver)	1✔
208
209
210	def smooth_acceleration(x, dt, params=None, options=None, gamma=None, window_size=None, solver=None):	1✔
211	"""Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative,
212	and then apply a convolutional gaussian smoother to the resulting derivative to smooth out the peaks.
213	The end result is similar to the jerk method, but can be more time-efficient.
214
215	:param np.array[float] x: data to differentiate
216	:param float dt: step size
217	:param params: (deprecated, prefer :code:`gamma` and :code:`window_size`)
218	:param dict options: (deprecated, prefer :code:`solver`) a dictionary consisting of {'solver': (str)}
219	:param float gamma: the regularization parameter
220	:param int window_size: window size for gaussian kernel
221	:param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc.
222	In testing, 'MOSEK' was the most robust. If not given, fall back to CVXPY's default.
223
224	:return: - x_hat (np.array) -- estimated (smoothed) x
225	- dxdt_hat (np.array) -- estimated derivative of x
226	"""
227	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
228	warn("`params` and `options` parameters will be removed in a future version. Use `gamma` " +	1✔
229	"and `solver` instead.", DeprecationWarning)
230	gamma, window_size = params	1✔
231	if options is not None:
232	if 'solver' in options: solver = options['solver']
233	elif gamma is None or window_size is None:	1✔
UNCOV 234	raise ValueError("`gamma` and `window_size` must be given.")	×
235
236	_, dxdt_hat = tvrdiff(x, dt, 2, gamma, solver=solver)	1✔
237
238	kernel = utility.gaussian_kernel(window_size)	1✔
239	dxdt_hat = utility.convolutional_smoother(dxdt_hat, kernel, 1)	1✔
240
241	x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt)	1✔
242	x0 = utility.estimate_integration_constant(x, x_hat)	1✔
243	x_hat = x_hat + x0	1✔
244
245	return x_hat, dxdt_hat	1✔

florisvb / PyNumDiff / 23028417751

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous