23027887759

Committed 12 Mar 2026 10:58PM UTC coverage: 90.516% (-1.3%) from 91.863%

Build # 23027887759

Build Type

Pull #188

github

Committed by

web-flow

Commit Message

Merge b9d73b18d into dff1f9c95

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

Coverage Stats

1 of 15 new or added lines in 1 file covered. (6.67%)

4 existing lines in 1 file now uncovered.

859 of 949 relevant lines covered (90.52%)

0.91 hits per line

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

78.65

/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

#N-d case:
def tvrdiff(x, dt, order, gamma, huberM=float('inf'), solver=None, axis=0):
    """
    Generalized total variation regularized derivatives (cvxpy). Supports multidimensionality by differentiating along
    'axis', independently for each vector obtained by fixing all other indices.
    
    :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.

    :return: - **x_hat** (np.array) -- estimated (smoothed) x
             - **dxdt_hat** (np.array) -- estimated derivative of x
    """

    x0 = np.moveaxis(x, axis, 0)

    # end quick if it's just 1d case
    if x0.ndim == 1:
        x_hat0, dxdt0 = tvrdiff(x0, dt, order, gamma, huberM, solver)
        return x_hat0, dxdt0
    
    x_hat0 = np.empty_like(x0, dtype=float)
    dxdt0 = np.empty_like(x0, dtype=float)
    rest = x0.shape[1:]
    print(rest)

    # had to loop in python:(
    for i in np.ndindex(rest):
        slice = (slice(None),) + i
        x_hat0[slice], dxdt0[slice] = tvrdiff(x0[slice], dt, order, gamma, huberM, solver)

    x_hat = np.moveaxis(x_hat0, 0, axis)
    dxdt_hat = np.moveaxis(dxdt0, 0, axis)

    return x_hat, dxdt_hat

# 1-d case:
def tvrdiff(x, dt, order, gamma, huberM=float('inf'), solver=None):
    """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.

    :return: - **x_hat** (np.array) -- estimated (smoothed) x
             - **dxdt_hat** (np.array) -- estimated derivative of x
    """

    # Normalize for numerical consistency with convex solver
    mu = np.mean(x)
    sigma = median_abs_deviation(x, scale='normal') # robust alternative to std()
    if sigma == 0: sigma = 1 # safety guard
    y = (x-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_hat = v/dt # v only holds the dx values; to get deriv scale by dt
    x_hat = 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_hat = (dxdt_hat[:-1] + dxdt_hat[1:])/2
    dxdt_hat = np.hstack((dxdt_hat, 2*dxdt_hat[-1] - dxdt_hat[-2])) # last value = penultimate value [-1] + diff between [-1] and [-2]

    return x_hat*sigma+mu, dxdt_hat*sigma # derivative is linear, so scale derivative by scatter


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	#N-d case:
56	def tvrdiff(x, dt, order, gamma, huberM=float('inf'), solver=None, axis=0):	1✔
57	"""
58	Generalized total variation regularized derivatives (cvxpy). Supports multidimensionality by differentiating along
59	'axis', independently for each vector obtained by fixing all other indices.
60
61	:param np.array[float] x: data to differentiate
62	:param float dt: step size
63	:param int order: 1, 2, or 3, the derivative to regularize
64	:param float gamma: regularization parameter
65	:param float huberM: Huber loss parameter, in units of scaled median absolute deviation of input data.
66	:math:`M = \\infty` reduces to :math:`\\ell_2` loss squared on first, fidelity cost term, and
67	:math:`M = 0` reduces to :math:`\\ell_1` loss, which seeks sparse residuals.
68	:param str solver: Solver to use. Solver options include: 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS'.
69	If not given, fall back to CVXPY's default.
70
71	:return: - x_hat (np.array) -- estimated (smoothed) x
72	- dxdt_hat (np.array) -- estimated derivative of x
73	"""
74
NEW 75	x0 = np.moveaxis(x, axis, 0)	×
76
77	# end quick if it's just 1d case
NEW 78	if x0.ndim == 1:	×
NEW 79	x_hat0, dxdt0 = tvrdiff(x0, dt, order, gamma, huberM, solver)	×
NEW 80	return x_hat0, dxdt0	×
81
NEW 82	x_hat0 = np.empty_like(x0, dtype=float)	×
NEW 83	dxdt0 = np.empty_like(x0, dtype=float)	×
NEW 84	rest = x0.shape[1:]	×
NEW 85	print(rest)	×
86
87	# had to loop in python:(
NEW 88	for i in np.ndindex(rest):	×
NEW 89	slice = (slice(None),) + i	×
NEW 90	x_hat0[slice], dxdt0[slice] = tvrdiff(x0[slice], dt, order, gamma, huberM, solver)	×
91
NEW 92	x_hat = np.moveaxis(x_hat0, 0, axis)	×
NEW 93	dxdt_hat = np.moveaxis(dxdt0, 0, axis)	×
94
NEW 95	return x_hat, dxdt_hat	×
96
97	# 1-d case:
98	def tvrdiff(x, dt, order, gamma, huberM=float('inf'), solver=None):	1✔
99	"""Generalized total variation regularized derivatives. Use convex optimization (cvxpy) to solve for a
100	total variation regularized derivative. Other convex-solver-based methods in this module call this function.
101
102	:param np.array[float] x: data to differentiate
103	:param float dt: step size
104	:param int order: 1, 2, or 3, the derivative to regularize
105	:param float gamma: regularization parameter
106	:param float huberM: Huber loss parameter, in units of scaled median absolute deviation of input data.
107	:math:`M = \\infty` reduces to :math:`\\ell_2` loss squared on first, fidelity cost term, and
108	:math:`M = 0` reduces to :math:`\\ell_1` loss, which seeks sparse residuals.
109	:param str solver: Solver to use. Solver options include: 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS'.
110	If not given, fall back to CVXPY's default.
111
112	:return: - x_hat (np.array) -- estimated (smoothed) x
113	- dxdt_hat (np.array) -- estimated derivative of x
114	"""
115
116	# Normalize for numerical consistency with convex solver
117	mu = np.mean(x)	1✔
118	sigma = median_abs_deviation(x, scale='normal') # robust alternative to std()	1✔
119	if sigma == 0: sigma = 1 # safety guard	1✔
120	y = (x-mu)/sigma	1✔
121
122	# Define the variables for the highest order derivative and the integration constants
123	deriv_values = cvxpy.Variable(len(y)) # values of the order^th derivative, in which we're penalizing variation	1✔
124	integration_constants = cvxpy.Variable(order) # constants of integration that help get us back to x	1✔
125
126	# Recursively integrate the highest order derivative to get back to the position. This is a first-
127	# order scheme, but it's very fast and tends to do not markedly worse than 2nd order. See #116
128	# I also tried a trapezoidal integration rule here, and it works no better. See #116 too.
129	hx = deriv_values # variables are integrated to produce the signal estimate variables, \hat{x} in the math	1✔
130	for i in range(order):	1✔
131	hx = cvxpy.cumsum(hx) + integration_constants[i] # cumsum is like integration assuming dt = 1	1✔
132
133	# Compare the recursively integrated position to the noisy position. \ell_2 doesn't get scaled by 1/2 here,
134	# so cvxpy's doubled Huber is already the right scale, and \ell_1 should be scaled by 2\sqrt{2} to match.
135	fidelity_cost = cvxpy.sum_squares(y - hx) if huberM == float('inf') \	1✔
136	else np.sqrt(8)*cvxpy.norm(y - hx, 1) if huberM == 0 \
137	else utility.huber_const(huberM)cvxpy.sum(cvxpy.huber(y - hx, huberM)) # data is already scaled, so M rather than Msigma
138	# Set up and solve the optimization problem
139	prob = cvxpy.Problem(cvxpy.Minimize(fidelity_cost + gamma*cvxpy.sum(cvxpy.tv(deriv_values)) ))	1✔
140	prob.solve(solver=solver)	1✔
141
142	# Recursively integrate the final derivative values to get back to the function and derivative values
143	v = deriv_values.value	1✔
144	for i in range(order-1): # stop one short to get the first derivative	1✔
145	v = np.cumsum(v) + integration_constants.value[i]	1✔
146	dxdt_hat = v/dt # v only holds the dx values; to get deriv scale by dt	1✔
147	x_hat = np.cumsum(v) + integration_constants.value[order-1] # smoothed data	1✔
148
149	# Due to the first-order nature of the derivative, it has a slight lag. Average together every two values
150	# to better center the answer. But this leaves us one-short, so devise a good last value.
151	dxdt_hat = (dxdt_hat[:-1] + dxdt_hat[1:])/2	1✔
152	dxdt_hat = np.hstack((dxdt_hat, 2*dxdt_hat[-1] - dxdt_hat[-2])) # last value = penultimate value [-1] + diff between [-1] and [-2]	1✔
153
154	return x_hatsigma+mu, dxdt_hatsigma # derivative is linear, so scale derivative by scatter	1✔
155
156
157	def velocity(x, dt, params=None, options=None, gamma=None, solver=None):	1✔
158	"""Use convex optimization (cvxpy) to solve for the velocity total variation regularized derivative.\n
159	Deprecated, prefer :code:`tvrdiff` with order 1 instead.
160
161	:param np.array[float] x: data to differentiate
162	:param float dt: step size
163	:param params: (deprecated, prefer :code:`gamma`)
164	:param dict options: (deprecated, prefer :code:`solver`) a dictionary consisting of {'solver': (str)}
165	:param float gamma: the regularization parameter
166	:param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc.
167	If not given, fall back to CVXPY's default.
168
169	:return: - x_hat (np.array) -- estimated (smoothed) x
170	- dxdt_hat (np.array) -- estimated derivative of x
171	"""
172	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
173	warn("`params` and `options` parameters will be removed in a future version. Use `gamma` " +	1✔
174	"and `solver` instead.", DeprecationWarning)
175	gamma = params[0] if isinstance(params, list) else params	1✔
176	if options is not None:
177	if 'solver' in options: solver = options['solver']
178	elif gamma is None:	1✔
UNCOV 179	raise ValueError("`gamma` must be given.")	×
180
181	warn("`velocity` is deprecated. Call `tvrdiff` with order 1 instead.", DeprecationWarning)	1✔
182	return tvrdiff(x, dt, 1, gamma, solver=solver)	1✔
183
184
185	def acceleration(x, dt, params=None, options=None, gamma=None, solver=None):	1✔
186	"""Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative.\n
187	Deprecated, prefer :code:`tvrdiff` with order 2 instead.
188
189	:param np.array[float] x: data to differentiate
190	:param float dt: step size
191	:param params: (deprecated, prefer :code:`gamma`)
192	:param dict options: (deprecated, prefer :code:`solver`) a dictionary consisting of {'solver': (str)}
193	:param float gamma: the regularization parameter
194	:param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc.
195	In testing, 'MOSEK' was the most robust. If not given, fall back to CVXPY's default.
196
197	:return: - x_hat (np.array) -- estimated (smoothed) x
198	- dxdt_hat (np.array) -- estimated derivative of x
199	"""
200	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
201	warn("`params` and `options` parameters will be removed in a future version. Use `gamma` " +	1✔
202	"and `solver` instead.", DeprecationWarning)
203	gamma = params[0] if isinstance(params, list) else params	1✔
204	if options is not None:
205	if 'solver' in options: solver = options['solver']
206	elif gamma is None:	1✔
UNCOV 207	raise ValueError("`gamma` must be given.")	×
208
209	warn("`acceleration` is deprecated. Call `tvrdiff` with order 2 instead.", DeprecationWarning)	1✔
210	return tvrdiff(x, dt, 2, gamma, solver=solver)	1✔
211
212
213	def jerk(x, dt, params=None, options=None, gamma=None, solver=None):	1✔
214	"""Use convex optimization (cvxpy) to solve for the jerk total variation regularized derivative.\n
215	Deprecated, prefer :code:`tvrdiff` with order 3 instead.
216
217	:param np.array[float] x: data to differentiate
218	:param float dt: step size
219	:param params: (deprecated, prefer :code:`gamma`)
220	:param dict options: (deprecated, prefer :code:`solver`) a dictionary consisting of {'solver': (str)}
221	:param float gamma: the regularization parameter
222	:param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc.
223	In testing, 'MOSEK' was the most robust. If not given, fall back to CVXPY's default.
224
225	:return: - x_hat (np.array) -- estimated (smoothed) x
226	- dxdt_hat (np.array) -- estimated derivative of x
227	"""
228	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
229	warn("`params` and `options` parameters will be removed in a future version. Use `gamma` " +	1✔
230	"and `solver` instead.", DeprecationWarning)
231	gamma = params[0] if isinstance(params, list) else params	1✔
232	if options is not None:
233	if 'solver' in options: solver = options['solver']
234	elif gamma is None:	1✔
UNCOV 235	raise ValueError("`gamma` must be given.")	×
236
237	warn("`jerk` is deprecated. Call `tvrdiff` with order 3 instead.", DeprecationWarning)	1✔
238	return tvrdiff(x, dt, 3, gamma, solver=solver)	1✔
239
240
241	def smooth_acceleration(x, dt, params=None, options=None, gamma=None, window_size=None, solver=None):	1✔
242	"""Use convex optimization (cvxpy) to solve for the acceleration total variation regularized derivative,
243	and then apply a convolutional gaussian smoother to the resulting derivative to smooth out the peaks.
244	The end result is similar to the jerk method, but can be more time-efficient.
245
246	:param np.array[float] x: data to differentiate
247	:param float dt: step size
248	:param params: (deprecated, prefer :code:`gamma` and :code:`window_size`)
249	:param dict options: (deprecated, prefer :code:`solver`) a dictionary consisting of {'solver': (str)}
250	:param float gamma: the regularization parameter
251	:param int window_size: window size for gaussian kernel
252	:param str solver: the solver CVXPY should use, 'MOSEK', 'CVXOPT', 'CLARABEL', 'ECOS', etc.
253	In testing, 'MOSEK' was the most robust. If not given, fall back to CVXPY's default.
254
255	:return: - x_hat (np.array) -- estimated (smoothed) x
256	- dxdt_hat (np.array) -- estimated derivative of x
257	"""
258	if params is not None: # Warning to support old interface for a while. Remove these lines along with params in a future release.	1✔
259	warn("`params` and `options` parameters will be removed in a future version. Use `gamma` " +	1✔
260	"and `solver` instead.", DeprecationWarning)
261	gamma, window_size = params	1✔
262	if options is not None:
263	if 'solver' in options: solver = options['solver']
264	elif gamma is None or window_size is None:	1✔
UNCOV 265	raise ValueError("`gamma` and `window_size` must be given.")	×
266
267	_, dxdt_hat = tvrdiff(x, dt, 2, gamma, solver=solver)	1✔
268
269	kernel = utility.gaussian_kernel(window_size)	1✔
270	dxdt_hat = utility.convolutional_smoother(dxdt_hat, kernel, 1)	1✔
271
272	x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt)	1✔
273	x0 = utility.estimate_integration_constant(x, x_hat)	1✔
274	x_hat = x_hat + x0	1✔
275
276	return x_hat, dxdt_hat	1✔

florisvb / PyNumDiff / 23027887759

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