utility

Simple, short, reusable helper functions

pynumdiff.utils.utility.convolutional_smoother(x, kernel, num_iterations=1, axis=0)

Perform smoothing by convolving x with a kernel.

Parameters:

  • x (np.array[float]) – 1D data

  • kernel (np.array[float]) – kernel to use in convolution

  • num_iterations (int) – number of iterations, >=1

  • axis (int) – data dimension along which convolution is performed

Returns:

x_hat (np.array[float]) – smoothed x

pynumdiff.utils.utility.estimate_integration_constant(x, x_hat, M=6, axis=0)

Integration leaves an unknown integration constant. This function finds a best fit integration constant to correct the DC of x_hat (the integral of dxdt_hat) by optimizing \(\min_c J(x - \hat{x} + c)\), where \(J\) is the Huber loss function or the \(\ell_1\) or \(\ell_2\) norm.

Parameters:

  • x (np.array[float]) – timeseries of measurements

  • x_hat (np.array[float]) – smoothed estimate of x

  • M (float) – constant estimation is robustified with the Huber loss. M here is in units of scaled mean absolute deviation of residuals, so scatter can be calculated and used to normalize without being thrown off by outliers. The default is intended to capture the idea of “six sigma”: Assuming Gaussian x - xhat errors, the portion of inliers beyond the Huber loss’ transition is only about 1.97e-9.

  • axis (int) – data dimension along which integration was performed

Returns:

integration constant (float or np.array[float]) – initial condition(s) to best align \(\mathbf{\hat{x}}\) with \(\mathbf{x}\)

pynumdiff.utils.utility.friedrichs_kernel(window_size)

A bump function :param int window_size: the width of the return value :return: kernel (np.array[float]) – samples of the kernel function

pynumdiff.utils.utility.gaussian_kernel(window_size)

A truncated gaussian :param int window_size: the width of the return value :return: kernel (np.array[float]) – samples of the kernel function

pynumdiff.utils.utility.huber_const(M)

Scale that makes sum(huber()) interpolate \(\sqrt{2}\|\cdot\|_1\) and \(\frac{1}{2}\|\cdot\|_2^2\), from https://jmlr.org/papers/volume14/aravkin13a/aravkin13a.pdf, with correction for missing sqrt. Here huber refers to scipy.special.huber.

Parameters:

M (float) – Huber parameter, where the function turns from quadratic to linear

Returns:

(float) – appropriate scale factor to normalize the Huber function

pynumdiff.utils.utility.integrate_dxdt_hat(dxdt_hat, dt_or_t, axis=0)

Wrapper for scipy.integrate.cumulative_trapezoid. Use 0 as first value so lengths match, see #88.

Parameters:

  • dxdt_hat (np.array[float]) – estimate derivative of timeseries

  • dt_or_t (float) – step size if given as a scalar or a vector of sample locations

  • axis (int) – data dimension along which to integrate

Returns:

x_hat (np.array[float]) – integral of dxdt_hat

pynumdiff.utils.utility.mean_kernel(window_size)

A uniform boxcar of total integral 1 :param int window_size: the width of the return value :return: kernel (np.array[float]) – samples of the kernel function

pynumdiff.utils.utility.peakdet(x, delta, t=None)

Find peaks and valleys of 1D array. A point is considered a maximum peak if it has the maximal value, and was preceded (to the left) by a value lower by delta. Converted from MATLAB script at http://billauer.co.il/peakdet.html Eli Billauer, 3.4.05 (Explicitly not copyrighted). This function is released to the public domain; Any use is allowed.

Parameters:

  • x (np.array[float]) – array for which to find peaks and valleys

  • delta (float) – threshold for finding peaks and valleys. A point is considered a maximum peak if it has the maximal value, and was preceded (to the left) by a value lower by delta.

  • t (np.array[float]) – optional domain points where data comes from, to make indices into locations

Returns:

  • maxtab – indices or locations (column 1) and values (column 2) of maxima

  • mintab – indices or locations (column 1) and values (column 2) of minima

pynumdiff.utils.utility.slide_function(func, x, dt_or_t, kernel, *args, stride=1, pass_weights=False, **kwargs)

Slide a smoothing derivative function across a timeseries with specified window size, and combine the results according to kernel weights.

Parameters:

  • func (callable) – name of the function to slide

  • x (np.array[float]) – data to differentiate

  • dt_or_t (float or np.array[float]) – constant step size (scalar) or array of sample locations (same length as x). When given as an array, the actual time values for each window are passed to func, enabling nonuniform spacing.

  • kernel (np.array[float]) – values to weight the sliding window

  • args (list) – passed to func

  • stride (int) – step size for slide (e.g. 1 means slide by 1 index location)

  • pass_weights (bool) – whether weights should be passed to func via update to kwargs

  • kwargs (dict) – passed to func

Returns:

  • x_hat – estimated (smoothed) x

  • dxdt_hat – estimated derivative of x