diff --git a/.gitignore b/.gitignore index 72cabbb..7ebc547 100644 --- a/.gitignore +++ b/.gitignore @@ -116,3 +116,4 @@ __pycache__ dist *.whl .vscode +docs/build diff --git a/README.md b/README.md index 2fc7894..8f1b481 100644 --- a/README.md +++ b/README.md @@ -9,8 +9,6 @@ Python methods for numerical differentiation of noisy data, including multi-obje

- - Documentation Status @@ -86,9 +84,8 @@ For more details, refer to [this paper](https://doi.org/10.1109/ACCESS.2020.3034 |- index.rst |- ... |- .gitignore - |- .travis.yml |- LICENSE.txt - |- requirements.txt + |- pyproject.toml ## Citation @@ -125,7 +122,7 @@ See CITATION.cff file as well as the following references. ### Prerequisite PyNumDiff requires common packages like `numpy`, `scipy`, `matplotlib`, `pytest` (for unittests), `pylint` -(for PEP8 style check). For a full list, you can check the file [requirements.txt](requirements.txt) +(for PEP8 style check). For a full list, you can check the file [pyproject.toml](pyproject.toml) In addition, it also requires certain additional packages for select functions, though these are not required for a successful install of PyNumDiff: * Total Variation Regularization methods: [`cvxpy`](http://www.cvxpy.org/install/index.html) diff --git a/docs/requirements.txt b/docs/requirements.txt deleted file mode 100644 index 9ec5929..0000000 --- a/docs/requirements.txt +++ /dev/null @@ -1,75 +0,0 @@ -appnope==0.1.2 -argon2-cffi==20.1.0 -astroid==2.6.2 -async-generator==1.10 -attrs==21.2.0 -backcall==0.2.0 -bleach==3.3.0 -certifi==2023.7.22 -cffi==1.14.5 -cvxopt==1.2.7 -cvxpy==1.1.13 -cycler==0.10.0 -debugpy==1.3.0 -decorator==5.0.9 -defusedxml==0.7.1 -ecos==2.0.7.post1 -entrypoints==0.3 -iniconfig==1.1.1 -ipykernel==6.0.1 -ipython==8.10.0 -ipython-genutils==0.2.0 -isort==5.9.1 -jedi==0.18.0 -Jinja2==3.0.1 -jsonschema==3.2.0 -jupyter-client==6.1.12 -jupyter-core==4.11.2 -jupyterlab-pygments==0.1.2 -kiwisolver==1.3.1 -lazy-object-proxy==1.6.0 -MarkupSafe==2.0.1 -matplotlib==3.4.2 -matplotlib-inline==0.1.2 -mccabe==0.6.1 -mistune==2.0.3 -Mosek==9.2.47 -nbclient==0.5.3 -nbconvert==6.5.1 -nbformat==5.1.3 -nest-asyncio==1.5.1 -notebook==6.4.12 -numpy==1.22.0 -osqp==0.6.2.post0 -packaging==21.0 -pandocfilters==1.4.3 -parso==0.8.2 -pexpect==4.8.0 -pickleshare==0.7.5 -Pillow==9.3.0 -pluggy==0.13.1 -prometheus-client==0.11.0 -prompt-toolkit==3.0.19 -ptyprocess==0.7.0 -py==1.10.0 -pycparser==2.20 -Pygments==2.15.0 -pylint==2.9.3 -pyparsing==2.4.7 -pyrsistent==0.18.0 -pytest==6.2.4 -python-dateutil==2.8.1 -pyzmq==22.1.0 -qdldl==0.1.5.post0 -scipy==1.7.0 -scs==2.1.4 -Send2Trash==1.7.1 -six==1.16.0 -terminado==0.10.1 -testpath==0.5.0 -toml==0.10.2 -tornado==6.3.3 -traitlets==5.0.5 -wcwidth==0.2.5 -webencodings==0.5.1 -wrapt==1.12.1 diff --git a/docs/source/LICENSE.rst b/docs/source/LICENSE.rst index b8349a0..c2124a3 100644 --- a/docs/source/LICENSE.rst +++ b/docs/source/LICENSE.rst @@ -6,6 +6,7 @@ The MIT License (MIT) Copyright (c) 2021 Floris van Breugel, The University of Nevada, Reno, Mechanical Engineering Department Yuying Liu, The University of Washington Applied Math Department +Pavel Komarov, Univerity of Washington Electrical and Computer Engineering Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/docs/source/conf.py b/docs/source/conf.py index ed50290..e441264 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -20,8 +20,8 @@ # -- Project information ----------------------------------------------------- project = 'PyNumDiff' -copyright = '2021, Floris van Breugel, Yuying Liu' -author = 'Floris van Breugel, Yuying Liu' +copyright = '2025, Floris van Breugel, Yuying Liu, Pavel Komarov' +author = 'Floris van Breugel, Yuying Liu, Pavel Komarov' # The full version, including alpha/beta/rc tags release = '0.1.2' diff --git a/docs/source/finite_difference.rst b/docs/source/finite_difference.rst index 638d5e7..c1e4747 100644 --- a/docs/source/finite_difference.rst +++ b/docs/source/finite_difference.rst @@ -1,7 +1,5 @@ finite_difference ================= -.. toctree:: - :maxdepth: 1 - - finite_difference/_finite_difference \ No newline at end of file +.. automodule:: pynumdiff.finite_difference + :members: diff --git a/docs/source/finite_difference/_finite_difference.rst b/docs/source/finite_difference/_finite_difference.rst deleted file mode 100644 index 8001c5c..0000000 --- a/docs/source/finite_difference/_finite_difference.rst +++ /dev/null @@ -1,7 +0,0 @@ -_finite_difference -===================== - -.. currentmodule:: pynumdiff.finite_difference._finite_difference - -.. automodule:: pynumdiff.finite_difference._finite_difference - :members: diff --git a/docs/source/index.rst b/docs/source/index.rst index 43629b6..0e486da 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -81,9 +81,8 @@ Structure |- ... |- setup.py |- .gitignore - |- .travis.yml |- LICENSE.txt - |- requirements.txt + |- pyproject.toml Getting Started --------------- diff --git a/docs/source/kalman_smooth.rst b/docs/source/kalman_smooth.rst index 7916e85..246ca27 100644 --- a/docs/source/kalman_smooth.rst +++ b/docs/source/kalman_smooth.rst @@ -1,7 +1,5 @@ kalman_smooth ============= -.. toctree:: - :maxdepth: 1 - - kalman_smooth/_kalman_smooth \ No newline at end of file +.. automodule:: pynumdiff.kalman_smooth + :members: diff --git a/docs/source/kalman_smooth/_kalman_smooth.rst b/docs/source/kalman_smooth/_kalman_smooth.rst deleted file mode 100644 index 9f3d391..0000000 --- a/docs/source/kalman_smooth/_kalman_smooth.rst +++ /dev/null @@ -1,7 +0,0 @@ -_kalman_smooth -================ - -.. currentmodule:: pynumdiff.kalman_smooth._kalman_smooth - -.. automodule:: pynumdiff.kalman_smooth._kalman_smooth - :members: \ No newline at end of file diff --git a/docs/source/linear_model.rst b/docs/source/linear_model.rst index fec816f..1a53152 100644 --- a/docs/source/linear_model.rst +++ b/docs/source/linear_model.rst @@ -1,7 +1,5 @@ linear_model ============ -.. toctree:: - :maxdepth: 1 - - linear_model/_linear_model \ No newline at end of file +.. automodule:: pynumdiff.linear_model + :members: \ No newline at end of file diff --git a/docs/source/linear_model/_linear_model.rst b/docs/source/linear_model/_linear_model.rst deleted file mode 100644 index 703453f..0000000 --- a/docs/source/linear_model/_linear_model.rst +++ /dev/null @@ -1,7 +0,0 @@ -_linear_model -================ - -.. currentmodule:: pynumdiff.linear_model._linear_model - -.. automodule:: pynumdiff.linear_model._linear_model - :members: diff --git a/docs/source/smooth_finite_difference.rst b/docs/source/smooth_finite_difference.rst index a841751..b48c9f3 100644 --- a/docs/source/smooth_finite_difference.rst +++ b/docs/source/smooth_finite_difference.rst @@ -1,7 +1,5 @@ smooth_finite_difference ======================== -.. toctree:: - :maxdepth: 1 - - smooth_finite_difference/_smooth_finite_difference \ No newline at end of file +.. automodule:: pynumdiff.smooth_finite_difference + :members: \ No newline at end of file diff --git a/docs/source/smooth_finite_difference/_smooth_finite_difference.rst b/docs/source/smooth_finite_difference/_smooth_finite_difference.rst deleted file mode 100644 index 5497e17..0000000 --- a/docs/source/smooth_finite_difference/_smooth_finite_difference.rst +++ /dev/null @@ -1,7 +0,0 @@ -_smooth_finite_difference -=========================== - -.. currentmodule:: pynumdiff.smooth_finite_difference._smooth_finite_difference - -.. automodule:: pynumdiff.smooth_finite_difference._smooth_finite_difference - :members: \ No newline at end of file diff --git a/docs/source/total_variation_regularization.rst b/docs/source/total_variation_regularization.rst index 0a57f81..ed75798 100644 --- a/docs/source/total_variation_regularization.rst +++ b/docs/source/total_variation_regularization.rst @@ -1,7 +1,5 @@ total_variation_regularization ============================== -.. toctree:: - :maxdepth: 1 - - total_variation_regularization/_total_variation_regularization \ No newline at end of file +.. automodule:: pynumdiff.total_variation_regularization + :members: \ No newline at end of file diff --git a/docs/source/total_variation_regularization/_total_variation_regularization.rst b/docs/source/total_variation_regularization/_total_variation_regularization.rst deleted file mode 100644 index bc6982e..0000000 --- a/docs/source/total_variation_regularization/_total_variation_regularization.rst +++ /dev/null @@ -1,7 +0,0 @@ -_total_variation_regularization -================================= - -.. currentmodule:: pynumdiff.total_variation_regularization._total_variation_regularization - -.. automodule:: pynumdiff.total_variation_regularization._total_variation_regularization - :members: \ No newline at end of file diff --git a/paper.bib b/paper.bib index 5fe0140..bf66c9a 100644 --- a/paper.bib +++ b/paper.bib @@ -225,7 +225,3 @@ @Misc{findiff year = {2018}, url = "https://github.com/maroba/findiff" } - - - - diff --git a/paper.md b/paper.md index 6e8d8cf..142ca8f 100644 --- a/paper.md +++ b/paper.md @@ -29,7 +29,7 @@ bibliography: paper.bib # Statement of need -The numerical computation of derivatives is ubiquitous in every scientific discipline and engineering application because derivatives express fundamental relationships among many quantities of interest. As a result, a large number of diverse algorithms have been developed to differentiate numerical data. These efforts are challenging because, in reality, practitioners often have sparse and noisy measurements and data, which undermine the ability to estimate accurate derivatives. Among the diversity of mathematical approaches that have been formulated, many are ad hoc in nature and require significant bespoke tuning of multiple parameters to produce reasonable results. Thus, at a practical level, it is often unclear which method should be used, how to choose parameters, and how to compare results from different methods. +The numerical computation of derivatives is ubiquitous in every scientific discipline and engineering application because derivatives express fundamental relationships among many quantities of interest. As a result, a large number of diverse algorithms have been developed to differentiate numerical data. These efforts are challenging because, in reality, practitioners often have sparse and noisy measurements and data, which undermine the ability to estimate accurate derivatives. Among the diversity of mathematical approaches that have been formulated, many are ad hoc in nature and require significant bespoke tuning of multiple parameters to produce reasonable results. Thus, at a practical level, it is often unclear which method should be used, how to choose parameters, and how to compare results from different methods. Regardless of application domain, scientists of various levels of mathematical expertise would benefit from a unified toolbox for differentiation techniques and parameter tuning. To address these needs, we built the open-source package `PyNumDiff`, with two primary goals in mind: (1) to develop a unified source for a diversity of differentiation methods using a common API, and (2) to provide an objective approach for choosing optimal parameters with a single universal hyperparameter (`gamma`) that functions similarly for all differentiation methods [@van2020numerical]. By filling these needs, `PyNumdiff` facilitates easy computations of derivatives on diverse time-series data sets. diff --git a/pynumdiff/finite_difference/__init__.py b/pynumdiff/finite_difference/__init__.py index 33de64b..77a6cae 100644 --- a/pynumdiff/finite_difference/__init__.py +++ b/pynumdiff/finite_difference/__init__.py @@ -1,5 +1,5 @@ +"""This module implements some common finite difference schemes """ -Import useful functions from _finite_difference module -""" -from pynumdiff.finite_difference._finite_difference import first_order -from pynumdiff.finite_difference._finite_difference import second_order +from ._finite_difference import first_order, second_order + +__all__ = ['first_order', 'second_order'] # So these get treated as direct members of the module by sphinx diff --git a/pynumdiff/finite_difference/_finite_difference.py b/pynumdiff/finite_difference/_finite_difference.py index 5728370..90eba03 100644 --- a/pynumdiff/finite_difference/_finite_difference.py +++ b/pynumdiff/finite_difference/_finite_difference.py @@ -1,65 +1,45 @@ -""" -This module implements some common finite difference schemes -""" import numpy as np from pynumdiff.utils import utility +from warnings import warn -def first_order(x, dt, params=None, options={}): - """ - First-order centered difference method - - :param x: array of time series to differentiate - :type x: np.array (floats) - - :param dt: time step size - :type dt: float - - :param params: number of iterations (if 'iterate' is enabled in options) - :type params: list (int) or int, optional - - :param options: a dictionary indicating whether to iteratively apply the difference method to smooth the estimates - :type options: dict {'iterate': boolean}, optional +def first_order(x, dt, params=None, options={}, num_iterations=None): + """First-order centered difference method - :return: a tuple consisting of: + :param np.array[float] x: array of time series to differentiate + :param float dt: time step size + :param list[float] or float params: (**deprecated**, prefer :code:`num_iterations`) + :param dict options: (**deprecated**, prefer :code:`num_iterations`) a dictionary consisting of {'iterate': (bool)} + :param int num_iterations: If performing iterated FD to smooth the estimates, give the number of iterations. + If ungiven, FD will not be iterated. - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x - - - :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: - assert params and isinstance(params, list), "params should be a non-empty list!" - return __iterate_first_order__(x, dt, params) + if params != None and 'iterate' in options: + warn("""`params` and `options` parameters will be removed in a future version. Use `num_iterations` instead.""") + if isinstance(params, list): params = params[0] + return _iterate_first_order(x, dt, params) + elif num_iterations: + return _iterate_first_order(x, dt, num_iterations) - # Calculate the finite difference - dxdt_hat = np.diff(x) / dt - # Pad the data - dxdt_hat = np.hstack((dxdt_hat[0], dxdt_hat, dxdt_hat[-1])) - # Re-finite dxdt_hat using linear interpolation - dxdt_hat = np.mean((dxdt_hat[0:-1], dxdt_hat[1:]), axis=0) + dxdt_hat = np.diff(x) / dt # Calculate the finite difference + dxdt_hat = np.hstack((dxdt_hat[0], dxdt_hat, dxdt_hat[-1])) # Pad the data + dxdt_hat = np.mean((dxdt_hat[0:-1], dxdt_hat[1:]), axis=0) # Re-finite dxdt_hat using linear interpolation return x, dxdt_hat def second_order(x, dt): - """ - Second-order centered difference method - - :param x: array of time series to differentiate - :type x: np.array (floats) - - :param dt: time step size - :type dt: float + """Second-order centered difference method - :return: a tuple consisting of: + :param np.array[float] x: array of time series to differentiate + :param float dt: time step size - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x - - - :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 """ dxdt_hat = (x[2:] - x[0:-2]) / (2 * dt) first_dxdt_hat = (-3 * x[0] + 4 * x[1] - x[2]) / (2 * dt) @@ -68,11 +48,8 @@ def second_order(x, dt): return x, dxdt_hat -def __x_hat_using_finite_difference__(x, dt): - """ - :param x: - :param dt: - :return: +def _x_hat_using_finite_difference(x, dt): + """Find a smoothed estimate of the true function by taking FD and then integrating with trapezoids """ x_hat, dxdt_hat = first_order(x, dt) x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt) @@ -80,40 +57,23 @@ def __x_hat_using_finite_difference__(x, dt): return x_hat + x0 -def __iterate_first_order__(x, dt, params): - """ - Iterative first order centered finite difference. - - :param x: array of time series to differentiate - :type x: np.array (floats) - - :param dt: time step size - :type dt: float - - :param params: number of iterations (if 'iterate' is enabled in options) - :type params: list (int) or int, optional - - :return: a tuple consisting of: +def _iterate_first_order(x, dt, num_iterations): + """Iterative first order centered finite difference. - - 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 int num_iterations: number of iterations - - :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 isinstance(params, list): - iterations = params[0] - else: - iterations = params - - # set up weights - w = np.arange(0, len(x), 1) - w = w / np.max(w) + w = np.arange(len(x)) / (len(x) - 1) # set up weights, [0., ... 1.0] # forward backward passes - for _ in range(iterations): - xf = __x_hat_using_finite_difference__(x, dt) - xb = __x_hat_using_finite_difference__(x[::-1], dt) + for _ in range(num_iterations): + xf = _x_hat_using_finite_difference(x, dt) + xb = _x_hat_using_finite_difference(x[::-1], dt) x = xf * w + xb[::-1] * (1 - w) x_hat, dxdt_hat = first_order(x, dt) diff --git a/pynumdiff/kalman_smooth/__init__.py b/pynumdiff/kalman_smooth/__init__.py index be3f896..6cf9932 100644 --- a/pynumdiff/kalman_smooth/__init__.py +++ b/pynumdiff/kalman_smooth/__init__.py @@ -1,9 +1,5 @@ +"""This module implements Kalman filters """ -Import useful functions from _finite_difference module -""" -from pynumdiff.kalman_smooth._kalman_smooth import constant_velocity -from pynumdiff.kalman_smooth._kalman_smooth import constant_acceleration -from pynumdiff.kalman_smooth._kalman_smooth import constant_jerk -from pynumdiff.kalman_smooth._kalman_smooth import known_dynamics -from pynumdiff.kalman_smooth._kalman_smooth import __kalman_forward_filter__ -from pynumdiff.kalman_smooth._kalman_smooth import __kalman_backward_smooth__ \ No newline at end of file +from ._kalman_smooth import constant_velocity, constant_acceleration, constant_jerk, known_dynamics, savgol_const_accel + +__all__ = ['constant_velocity', 'constant_acceleration', 'constant_jerk', 'known_dynamics', 'savgol_const_accel'] # So these get treated as direct members of the module by sphinx diff --git a/pynumdiff/kalman_smooth/_kalman_smooth.py b/pynumdiff/kalman_smooth/_kalman_smooth.py index ecd545f..57c1aae 100644 --- a/pynumdiff/kalman_smooth/_kalman_smooth.py +++ b/pynumdiff/kalman_smooth/_kalman_smooth.py @@ -1,6 +1,3 @@ -""" -This module implements Kalman filters -""" import copy import numpy as np diff --git a/pynumdiff/linear_model/__init__.py b/pynumdiff/linear_model/__init__.py index 3e8e775..933d074 100644 --- a/pynumdiff/linear_model/__init__.py +++ b/pynumdiff/linear_model/__init__.py @@ -1,49 +1,13 @@ +"""This module implements interpolation-based differentiation schemes. """ -Import useful functions from _linear_model module -""" -import logging as _logging -_logging.basicConfig( - level=_logging.INFO, - format="%(asctime)s [%(levelname)s] %(message)s", - handlers=[ - _logging.FileHandler("debug.log"), - _logging.StreamHandler() - ] -) - -# try: -# import pychebfun -# __pychebfun_installed__ = True -# except: -# __pychebfun_installed__ = False - try: import cvxpy - __cvxpy_installed__ = True + from ._linear_model import lineardiff except: - __cvxpy_installed__ = False - - -# if __pychebfun_installed__: -# from pynumdiff.linear_model._linear_model import chebydiff -# else: -# __warning__ = '\nLimited Linear Model Support Detected! \n'\ -# '---> PYCHEBFUN is not installed. \n'\ -# '---> Install pychebfun to use chebfun derivatives (https://github.com/pychebfun/pychebfun/) \n'\ -# 'You can still use other methods \n' -# _logging.info(__warning__) + from warnings import warn + warn("""Limited Linear Model Support Detected! CVXPY is not installed. + Install CVXPY to use lineardiff derivatives. You can still use other methods.""") -if __cvxpy_installed__: - from pynumdiff.linear_model._linear_model import lineardiff - from pynumdiff.linear_model._linear_model import __integrate_dxdt_hat_matrix__ - from pynumdiff.linear_model._linear_model import __solve_for_A_and_C_given_X_and_Xdot__ -else: - __warning__ = '\nLimited Linear Model Support Detected! \n'\ - '---> CVXPY is not installed. \n'\ - '---> Install CVXPY to use lineardiff derivatives \n'\ - 'You can still use other methods \n' - _logging.info(__warning__) +from ._linear_model import savgoldiff, polydiff, spectraldiff -from pynumdiff.linear_model._linear_model import savgoldiff -from pynumdiff.linear_model._linear_model import spectraldiff -from pynumdiff.linear_model._linear_model import polydiff \ No newline at end of file +__all__ = ['lineardiff', 'savgoldiff', 'polydiff', 'spectraldiff'] # So these get treated as direct members of the module by sphinx diff --git a/pynumdiff/linear_model/_linear_model.py b/pynumdiff/linear_model/_linear_model.py index 1b35822..e736a7d 100644 --- a/pynumdiff/linear_model/_linear_model.py +++ b/pynumdiff/linear_model/_linear_model.py @@ -1,36 +1,22 @@ -""" -This module implements some common finite difference schemes -""" -import copy -import math -import logging -import scipy +import copy, math, logging, scipy import numpy as np +from warnings import warn from pynumdiff import smooth_finite_difference from pynumdiff.finite_difference import first_order as finite_difference from pynumdiff.utils import utility -# try: -# import pychebfun -# except ImportError: -# pass - try: import cvxpy except ImportError: pass -__friedrichs_kernel__ = utility.__friedrichs_kernel__ -__gaussian_kernel__ = utility.__gaussian_kernel__ -KERNELS = {'friedrichs': __friedrichs_kernel__, - 'gaussian': __gaussian_kernel__} +KERNELS = {'friedrichs': utility._friedrichs_kernel, + 'gaussian': utility._gaussian_kernel} #################### # Helper functions # #################### - - def __slide_function__(func, x, dt, params, window_size, step_size, kernel_name, solver=None): """ Slide a smoothing derivative function across a timeseries with specified window size. @@ -109,8 +95,6 @@ def __slide_function__(func, x, dt, params, window_size, step_size, kernel_name, ######################### # Savitzky-Golay filter # ######################### - - def savgoldiff(x, dt, params, options=None): """ Use the Savitzky-Golay to smooth the data and calculate the first derivative. It wses scipy.signal.savgol_filter. The Savitzky-Golay is very similar to the sliding polynomial fit, but slightly noisier, and much faster @@ -153,8 +137,8 @@ def savgoldiff(x, dt, params, options=None): dxdt_hat = scipy.signal.savgol_filter(x, window_size, n, deriv=1) / dt - kernel = __gaussian_kernel__(smoothing_win) - dxdt_hat = smooth_finite_difference.__convolutional_smoother__(dxdt_hat, kernel, 1) + kernel = utility._gaussian_kernel(smoothing_win) + dxdt_hat = utility.convolutional_smoother(dxdt_hat, kernel, 1) x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt) x0 = utility.estimate_initial_condition(x, x_hat) @@ -166,8 +150,6 @@ def savgoldiff(x, dt, params, options=None): ###################### # Polynomial fitting # ###################### - - def __polydiff__(x, dt, params, options=None): """ Fit polynomials to the timeseries, and differentiate the polynomials. @@ -439,7 +421,6 @@ def __lineardiff__(x, dt, params, options=None): x = x - mean - # Generate the matrix of integrals of x X = [x] for n in range(1, N): @@ -546,52 +527,40 @@ def lineardiff(x, dt, params, options=None): return __lineardiff__(x, dt, params, options) -####################### -# Spectral derivative # -####################### - - -def spectraldiff(x, dt, params, options=None): - """ - Take a derivative in the fourier domain, with high frequency attentuation. - - :param x: array of time series to differentiate - :type x: np.array (float) - - :param dt: time step size - :type dt: float - - :param params: the high frequency cut off - - :type params: list (float) or float - - :param options: a dictionary consisting of 2 key value pairs: - - - 'even_extension': if True, extend the time series with an even extension so signal starts and ends at the same value. - - 'pad_to_zero_dxdt': if True, extend the time series with extensions that smoothly force the derivative to zero. This allows the spectral derivative to fit data which does not start and end with derivatives equal to zero. - - :type options: dict {'even_extension': (bool), 'pad_to_zero_dxdt': (bool)}, optional - - :return: a tuple consisting of: - - - x_hat: estimated (smoothed) x - - dxdt_hat: estimated derivative of x - - :rtype: tuple -> (np.array, np.array) +############################### +# Fourier Spectral derivative # +############################### +def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_extension=True, pad_to_zero_dxdt=True): + """Take a derivative in the fourier domain, with high frequency attentuation. + + :param np.array[float] x: array of time series to differentiate + :param float dt: time step size + :param list[float] or float params: (**deprecated**, prefer :code:`high_freq_cutoff`) + :param dict options: (**deprecated**, prefer :code:`even_extension` + and :code:`pad_to_zero_dxdt`) a dictionary consisting of {'even_extension': (bool), 'pad_to_zero_dxdt': (bool)} + :param float high_freq_cutoff: The high frequency cutoff. Frequencies below this threshold will be kept, and above will be zeroed. + :param bool even_extension: if True, extend the time series with an even extension so signal starts and ends at the same value. + :param bool pad_to_zero_dxdt: if True, extend the time series with extensions that smoothly force the derivative to zero. This + allows the spectral derivative to fit data which does not start and end with derivatives equal to zero. + + :return: tuple[np.array, np.array] of\n + - **x_hat** -- estimated (smoothed) x + - **dxdt_hat** -- estimated derivative of x """ - - if options is None: - options = {'even_extension': True, 'pad_to_zero_dxdt': True} - - if isinstance(params, list): - wn = params[0] - else: - 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 `high_freq_cutoff`, + `even_extension`, and `pad_to_zero_dxdt` instead.""") + if options is None: + even_extension = True + pad_to_zero_dxdt = True + high_freq_cutoff = params[0] if isinstance(params, list) else params + elif high_freq_cutoff == None: + raise ValueError("`high_freq_cutoff` must be given.") original_L = len(x) # make derivative go to zero at ends (optional) - if options['pad_to_zero_dxdt']: + if pad_to_zero_dxdt: padding = 100 pre = x[0]*np.ones(padding) post = x[-1]*np.ones(padding) @@ -603,7 +572,7 @@ def spectraldiff(x, dt, params, options=None): padding = 0 # Do even extension (optional) - if options['even_extension'] is True: + if even_extension is True: x = np.hstack((x, x[::-1])) # If odd, make N even, and pad x @@ -615,12 +584,12 @@ def spectraldiff(x, dt, params, options=None): N = L # Define the frequency range. - k = np.asarray(list(range(0, int(N/2))) + [0] + list(range( int(-N/2) + 1,0))) + k = np.asarray(list(range(0, int(N/2))) + [0] + list(range(int(-N/2) + 1,0))) k = k*2*np.pi/(dt*N) - # Frequency based smoothing: remove signals with a frequency higher than wn - discrete_wn = int(wn*N) - k[discrete_wn:N-discrete_wn] = 0 + # Frequency based smoothing: remove signals with a frequency higher than high_freq_cutoff + discrete_high_freq_cutoff = int(high_freq_cutoff*N) + k[discrete_high_freq_cutoff:N-discrete_high_freq_cutoff] = 0 # Derivative = 90 deg phase shift dxdt_hat = np.real(np.fft.ifft(1.0j * k * np.fft.fft(x))) diff --git a/pynumdiff/smooth_finite_difference/__init__.py b/pynumdiff/smooth_finite_difference/__init__.py index 2b6f0e9..25ed6f6 100644 --- a/pynumdiff/smooth_finite_difference/__init__.py +++ b/pynumdiff/smooth_finite_difference/__init__.py @@ -1,10 +1,5 @@ +"""Apply smoothing method before finite difference. """ -Import useful functions from _smooth_finite_difference -""" -from pynumdiff.smooth_finite_difference._smooth_finite_difference import mediandiff -from pynumdiff.smooth_finite_difference._smooth_finite_difference import meandiff -from pynumdiff.smooth_finite_difference._smooth_finite_difference import gaussiandiff -from pynumdiff.smooth_finite_difference._smooth_finite_difference import friedrichsdiff -from pynumdiff.smooth_finite_difference._smooth_finite_difference import butterdiff -from pynumdiff.smooth_finite_difference._smooth_finite_difference import splinediff -from pynumdiff.smooth_finite_difference._smooth_finite_difference import __convolutional_smoother__ +from ._smooth_finite_difference import mediandiff, meandiff, gaussiandiff, friedrichsdiff, butterdiff, splinediff + +__all__ = ['mediandiff', 'meandiff', 'gaussiandiff', 'friedrichsdiff', 'butterdiff', 'splinediff'] # So these get treated as direct members of the module by sphinx diff --git a/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py b/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py index 6a86559..58a4cde 100644 --- a/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py +++ b/pynumdiff/smooth_finite_difference/_smooth_finite_difference.py @@ -1,61 +1,14 @@ -""" -Apply smoothing method before finite difference -""" import numpy as np import scipy.signal # included code from pynumdiff.finite_difference import first_order as finite_difference from pynumdiff.utils import utility -__friedrichs_kernel__ = utility.__friedrichs_kernel__ -__gaussian_kernel__ = utility.__gaussian_kernel__ -__mean_kernel__ = utility.__mean_kernel__ -##################### -# Smoothing methods # -##################### - - -def __median_smooth__(x, window_size): - """ - :param x: - :param window_size: - :return: - """ - assert window_size % 2 == 1 # is odd - x_hat = scipy.signal.medfilt(x, window_size) - return x_hat - - -# convolve kernels -def __convolutional_smoother__(x, kernel, iterations): - """ - Perform mean smoothing by convolving mean kernel with x - followed by first order finite difference - - :param x: (np.array of floats, 1xN) time series to differentiate - :param kernel: (np.array of floats, 1 x window_size) kernel to use in convolution - :param iterations: (int) number of iterations, >=1 - :return: x_hat : smoothed x - """ - # pad - x_hat = np.hstack((x[::-1], x, x[::-1])) - for _ in range(iterations): - x_hat_f = np.convolve(x_hat, kernel, 'same') - x_hat_b = np.convolve(x_hat[::-1], kernel, 'same')[::-1] - - w = np.arange(0, len(x_hat_f), 1) - w = w/np.max(w) - x_hat = x_hat_f*w + x_hat_b*(1-w) - - return x_hat[len(x):len(x)*2] - ################################ # Smoothing finite differences # ################################ - - def mediandiff(x, dt, params, options={}): """ Perform median smoothing using scipy.signal.medfilt @@ -95,11 +48,11 @@ def mediandiff(x, dt, params, options={}): window_size = params if not window_size % 2: - window_size += 1 + window_size += 1 # assert window_size % 2 == 1 # is odd x_hat = x for _ in range(iterations): - x_hat = __median_smooth__(x_hat, window_size) + x_hat = scipy.signal.medfilt(x_hat, window_size) x_hat, dxdt_hat = finite_difference(x_hat, dt) return x_hat, dxdt_hat @@ -110,11 +63,8 @@ def meandiff(x, dt, params, options={}): Perform mean smoothing by convolving mean 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 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] @@ -145,8 +95,8 @@ def meandiff(x, dt, params, options={}): else: window_size = params - kernel = __mean_kernel__(window_size) - x_hat = __convolutional_smoother__(x, kernel, iterations) + kernel = utility._mean_kernel(window_size) + x_hat = utility.convolutional_smoother(x, kernel, iterations) x_hat, dxdt_hat = finite_difference(x_hat, dt) return x_hat, dxdt_hat @@ -191,8 +141,8 @@ def gaussiandiff(x, dt, params, options={}): else: window_size = params - kernel = __gaussian_kernel__(window_size) - x_hat = __convolutional_smoother__(x, kernel, iterations) + kernel = utility._gaussian_kernel(window_size) + x_hat = utility.convolutional_smoother(x, kernel, iterations) x_hat, dxdt_hat = finite_difference(x_hat, dt) return x_hat, dxdt_hat @@ -238,8 +188,8 @@ def friedrichsdiff(x, dt, params, options={}): else: window_size = params - kernel = __friedrichs_kernel__(window_size) - x_hat = __convolutional_smoother__(x, kernel, iterations) + kernel = utility._friedrichs_kernel(window_size) + x_hat = utility.convolutional_smoother(x, kernel, iterations) x_hat, dxdt_hat = finite_difference(x_hat, dt) return x_hat, dxdt_hat diff --git a/pynumdiff/tests/test_finite_difference.py b/pynumdiff/tests/test_finite_difference.py index 9380672..99a6272 100644 --- a/pynumdiff/tests/test_finite_difference.py +++ b/pynumdiff/tests/test_finite_difference.py @@ -3,42 +3,40 @@ """ # pylint: skip-file import numpy as np -from unittest import TestCase from pynumdiff import first_order, second_order -class TestFD(TestCase): - def test_first_order_1(self): - x = np.array([1, 2, 3, 4, 5]) - dt = 0.01 - dxdt = np.array([100, 100, 100, 100, 100]) - _, dxdt_hat = first_order(x, dt) - np.testing.assert_array_equal(dxdt_hat, dxdt) +def test_first_order_1(): + x = np.array([1, 2, 3, 4, 5]) + dt = 0.01 + dxdt = np.array([100, 100, 100, 100, 100]) + _, dxdt_hat = first_order(x, dt) + np.testing.assert_array_equal(dxdt_hat, dxdt) - def test_first_order_2(self): - x = np.array([8, 3, 14, 0, 9]) - dt = 0.01 - dxdt = np.array([-500, 300, -150, -250, 900]) - _, dxdt_hat = first_order(x, dt) - np.testing.assert_array_equal(dxdt_hat, dxdt) +def test_first_order_2(): + x = np.array([8, 3, 14, 0, 9]) + dt = 0.01 + dxdt = np.array([-500, 300, -150, -250, 900]) + _, dxdt_hat = first_order(x, dt) + np.testing.assert_array_equal(dxdt_hat, dxdt) - def test_first_order_iterative(self): - x = np.random.rand(100) - dt = 0.01 - params = [100] - _, dxdt_hat = first_order(x, dt, params, options={'iterate': True}) - assert x.shape == dxdt_hat.shape +def test_first_order_iterative(): + x = np.random.rand(100) + dt = 0.01 + params = [100] + _, dxdt_hat = first_order(x, dt, params, options={'iterate': True}) + assert x.shape == dxdt_hat.shape - def test_second_order_1(self): - x = np.array([1, 2, 3, 4, 5]) - dt = 0.01 - dxdt = np.array([100, 100, 100, 100, 100]) - _, dxdt_hat = second_order(x, dt) - np.testing.assert_array_equal(dxdt_hat, dxdt) +def test_second_order_1(): + x = np.array([1, 2, 3, 4, 5]) + dt = 0.01 + dxdt = np.array([100, 100, 100, 100, 100]) + _, dxdt_hat = second_order(x, dt) + np.testing.assert_array_equal(dxdt_hat, dxdt) - def test_second_order_2(self): - x = np.array([8, 3, 14, 0, 9]) - dt = 0.01 - dxdt = np.array([-1300, 300, -150, -250, 2050]) - _, dxdt_hat = second_order(x, dt) - np.testing.assert_array_equal(dxdt_hat, dxdt) +def test_second_order_2(): + x = np.array([8, 3, 14, 0, 9]) + dt = 0.01 + dxdt = np.array([-1300, 300, -150, -250, 2050]) + _, dxdt_hat = second_order(x, dt) + np.testing.assert_array_equal(dxdt_hat, dxdt) diff --git a/pynumdiff/tests/test_kalman_smooth.py b/pynumdiff/tests/test_kalman_smooth.py index 6d444c8..df8756c 100644 --- a/pynumdiff/tests/test_kalman_smooth.py +++ b/pynumdiff/tests/test_kalman_smooth.py @@ -4,7 +4,6 @@ # pylint: skip-file import numpy as np -from unittest import TestCase from pynumdiff.kalman_smooth import constant_velocity, constant_acceleration, \ constant_jerk, known_dynamics @@ -15,44 +14,41 @@ dt = 0.01 -class TestKS(TestCase): - def test_constant_velocity(self): - params = [1e-4, 1e-5] - x_hat, dxdt_hat = constant_velocity(x, dt, params) - x_smooth = np.array([7.952849, 7.714494, 7.769948, 7.81768, 8.330625, 8.332996, - 8.46594, 8.243244, 8.458473, 8.367324, 8.284892, 7.947729, - 7.998362, 8.123646, 8.303191]) - dxdt = np.array([88.750804, 93.567378, 102.828004, 62.994815, 92.01605, - 60.395089, 47.494064, 27.626483, 21.537133, 14.105156, - 8.138253, 7.996629, 4.016616, -0.1122, 2.319358]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) - - def test_constant_acceleration(self): - params = [1e-4, 1e-2] - x_hat, dxdt_hat = constant_acceleration(x, dt, params) - x_smooth = np.array([5.069524, 6.137091, 7.191819, 8.062104, 9.112695, 9.616349, - 10.029422, 9.945811, 10.05048, 9.703503, 9.180588, 8.309991, - 7.546839, 6.581819, 5.421122]) - dxdt = np.array([170.225553, 164.483647, 156.524187, 103.452558, 113.776639, - 64.258467, 33.813842, -1.889904, -25.372839, -48.272303, - -69.60202, -81.885049, -101.379641, -122.551681, -140.214842]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) - - def test_constant_jerk(self): - params = [1e-4, 1e-4] - x_hat, dxdt_hat = constant_jerk(x, dt, params) - x_smooth = np.array([5.066536, 6.135826, 7.191131, 8.061294, 9.110784, 9.613802, - 10.026445, 9.943029, 10.047933, 9.701807, 9.179971, 8.310492, - 7.547672, 6.582594, 5.421728]) - dxdt = np.array([170.262874, 164.484367, 156.478206, 103.371112, 113.682324, - 64.169044, 33.742701, -1.935552, -25.398252, -48.273806, - -69.59001, -81.873115, -101.384521, -122.579907, -140.269899]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) - - def test_known_dynamics(self): - return - - +def test_constant_velocity(): + params = [1e-4, 1e-5] + x_hat, dxdt_hat = constant_velocity(x, dt, params) + x_smooth = np.array([7.952849, 7.714494, 7.769948, 7.81768, 8.330625, 8.332996, + 8.46594, 8.243244, 8.458473, 8.367324, 8.284892, 7.947729, + 7.998362, 8.123646, 8.303191]) + dxdt = np.array([88.750804, 93.567378, 102.828004, 62.994815, 92.01605, + 60.395089, 47.494064, 27.626483, 21.537133, 14.105156, + 8.138253, 7.996629, 4.016616, -0.1122, 2.319358]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) + +def test_constant_acceleration(): + params = [1e-4, 1e-2] + x_hat, dxdt_hat = constant_acceleration(x, dt, params) + x_smooth = np.array([5.069524, 6.137091, 7.191819, 8.062104, 9.112695, 9.616349, + 10.029422, 9.945811, 10.05048, 9.703503, 9.180588, 8.309991, + 7.546839, 6.581819, 5.421122]) + dxdt = np.array([170.225553, 164.483647, 156.524187, 103.452558, 113.776639, + 64.258467, 33.813842, -1.889904, -25.372839, -48.272303, + -69.60202, -81.885049, -101.379641, -122.551681, -140.214842]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) + +def test_constant_jerk(): + params = [1e-4, 1e-4] + x_hat, dxdt_hat = constant_jerk(x, dt, params) + x_smooth = np.array([5.066536, 6.135826, 7.191131, 8.061294, 9.110784, 9.613802, + 10.026445, 9.943029, 10.047933, 9.701807, 9.179971, 8.310492, + 7.547672, 6.582594, 5.421728]) + dxdt = np.array([170.262874, 164.484367, 156.478206, 103.371112, 113.682324, + 64.169044, 33.742701, -1.935552, -25.398252, -48.273806, + -69.59001, -81.873115, -101.384521, -122.579907, -140.269899]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=3) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=3) + +def test_known_dynamics(): + return#TODO. We shouldn't have empty tests hanging around diff --git a/pynumdiff/tests/test_linear_model.py b/pynumdiff/tests/test_linear_model.py index cf577c0..0c8ed29 100644 --- a/pynumdiff/tests/test_linear_model.py +++ b/pynumdiff/tests/test_linear_model.py @@ -4,7 +4,6 @@ # pylint: skip-file import numpy as np -from unittest import TestCase import logging as _logging _logging.basicConfig( @@ -26,76 +25,75 @@ dt = 0.01 -class TestLM(TestCase): - def test_savgoldiff(self): - params = [2, 4, 4] - x_hat, dxdt_hat = savgoldiff(x, dt, params) - x_smooth = np.array([4.669816, 4.374363, 6.46848, 8.899164, 10.606681, 11.059424, - 10.519283, 10.058375, 10.191014, 10.193343, 9.208019, 7.445465, - 5.880869, 5.49672, 6.930156]) - dxdt = np.array([-29.5453, 156.853147, 261.970245, 224.16657, 117.336993, - -26.788542, -81.239512, -10.942197, 37.470096, -37.004311, - -160.060586, -192.450136, -120.46908, 43.639278, 243.047964]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) +def test_savgoldiff(): + params = [2, 4, 4] + x_hat, dxdt_hat = savgoldiff(x, dt, params) + x_smooth = np.array([4.669816, 4.374363, 6.46848, 8.899164, 10.606681, 11.059424, + 10.519283, 10.058375, 10.191014, 10.193343, 9.208019, 7.445465, + 5.880869, 5.49672, 6.930156]) + dxdt = np.array([-29.5453, 156.853147, 261.970245, 224.16657, 117.336993, + -26.788542, -81.239512, -10.942197, 37.470096, -37.004311, + -160.060586, -192.450136, -120.46908, 43.639278, 243.047964]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - def test_spectraldiff(self): - params = [0.1] - x_hat, dxdt_hat = spectraldiff(x, dt, params) - x_smooth = np.array([3.99, 5.038, 6.635, 8.365, 9.971, 11.201, 11.86, 11.86, - 11.231, 10.113, 8.722, 7.296, 6.047, 5.116, 4.556]) - dxdt = np.array([104.803, 147., 172.464, 173.547, 147.67, 98.194, - 33.754, -33.769, -92.105, -131.479, -146.761, -138.333, - -111.508, -74.752, -37.276]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) +def test_spectraldiff(): + params = [0.1] + x_hat, dxdt_hat = spectraldiff(x, dt, params) + x_smooth = np.array([3.99, 5.038, 6.635, 8.365, 9.971, 11.201, 11.86, 11.86, + 11.231, 10.113, 8.722, 7.296, 6.047, 5.116, 4.556]) + dxdt = np.array([104.803, 147., 172.464, 173.547, 147.67, 98.194, + 33.754, -33.769, -92.105, -131.479, -146.761, -138.333, + -111.508, -74.752, -37.276]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - def test_polydiff(self): - params = [2, 3] - x_hat, dxdt_hat = polydiff(x, dt, params) - x_smooth = np.array([1.16153, 4.506877, 6.407802, 8.544663, 13.431766, 14.051294, - 10.115687, 7.674865, 10.471466, 13.612046, 11.363571, 5.68407, - 4.443968, 6.213507, 4.695931]) - dxdt = np.array([330.730385, 284.267456, 299.891801, 305.441626, 205.475727, - -145.229037, -279.41178, 15.428548, 244.252341, 20.343789, - -326.727498, -288.988297, 33.647456, 27.861175, -344.695033]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) +def test_polydiff(): + params = [2, 3] + x_hat, dxdt_hat = polydiff(x, dt, params) + x_smooth = np.array([1.16153, 4.506877, 6.407802, 8.544663, 13.431766, 14.051294, + 10.115687, 7.674865, 10.471466, 13.612046, 11.363571, 5.68407, + 4.443968, 6.213507, 4.695931]) + dxdt = np.array([330.730385, 284.267456, 299.891801, 305.441626, 205.475727, + -145.229037, -279.41178, 15.428548, 244.252341, 20.343789, + -326.727498, -288.988297, 33.647456, 27.861175, -344.695033]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - # def test_chebydiff(self): - # try: - # import pychebfun - # except: - # __warning__ = '\nCannot import pychebfun, skipping chebydiff test.' - # _logging.info(__warning__) - # return +# def test_chebydiff(self): +# try: +# import pychebfun +# except: +# __warning__ = '\nCannot import pychebfun, skipping chebydiff test.' +# _logging.info(__warning__) +# return - # params = [2, 3] - # x_hat, dxdt_hat = chebydiff(x, dt, params) - # x_smooth = np.array([1., 4.638844, 7.184256, 6.644655, 15.614775, 11.60484, - # 12.284141, 6.082226, 12.000615, 12.058705, 12.462283, 5.018101, - # 4.674378, 6.431201, 4.]) - # dxdt = np.array([202.732652, 346.950235, -140.713336, 498.719617, 212.717775, - # -185.13847, -266.604056, -51.792587, 377.969849, -0.749768, - # -297.654931, -455.876155, 197.575692, -24.809441, -150.109487]) - # np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - # np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) +# params = [2, 3] +# x_hat, dxdt_hat = chebydiff(x, dt, params) +# x_smooth = np.array([1., 4.638844, 7.184256, 6.644655, 15.614775, 11.60484, +# 12.284141, 6.082226, 12.000615, 12.058705, 12.462283, 5.018101, +# 4.674378, 6.431201, 4.]) +# dxdt = np.array([202.732652, 346.950235, -140.713336, 498.719617, 212.717775, +# -185.13847, -266.604056, -51.792587, 377.969849, -0.749768, +# -297.654931, -455.876155, 197.575692, -24.809441, -150.109487]) +# np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) +# np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - def test_lineardiff(self): - try: - import cvxpy - except: - __warning__ = '\nCannot import cvxpy, skipping lineardiff test.' - _logging.info(__warning__) - return +def test_lineardiff(): + try: + import cvxpy + except: + __warning__ = '\nCannot import cvxpy, skipping lineardiff test.' + _logging.info(__warning__) + return - params = [3, 5, 10] - x_hat, dxdt_hat = lineardiff(x, dt, params, options={'solver': 'CVXOPT'}) - x_smooth = np.array([3.070975, 3.435072, 6.363585, 10.276584, 12.033974, 10.594136, - 9.608228, 9.731326, 10.333255, 10.806791, 9.710448, 7.456045, - 5.70695, 4.856271, 5.685251]) - dxdt = np.array([36.409751, 164.630545, 342.075623, 283.519415, 15.877598, - -121.287252, -43.140514, 36.251305, 53.773231, -31.140351, - -167.537258, -200.174883, -129.988725, -1.084955, 82.897991]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) + params = [3, 5, 10] + x_hat, dxdt_hat = lineardiff(x, dt, params, options={'solver': 'CVXOPT'}) + x_smooth = np.array([3.070975, 3.435072, 6.363585, 10.276584, 12.033974, 10.594136, + 9.608228, 9.731326, 10.333255, 10.806791, 9.710448, 7.456045, + 5.70695, 4.856271, 5.685251]) + dxdt = np.array([36.409751, 164.630545, 342.075623, 283.519415, 15.877598, + -121.287252, -43.140514, 36.251305, 53.773231, -31.140351, + -167.537258, -200.174883, -129.988725, -1.084955, 82.897991]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) diff --git a/pynumdiff/tests/test_optimize.py b/pynumdiff/tests/test_optimize.py index 0256106..633153d 100644 --- a/pynumdiff/tests/test_optimize.py +++ b/pynumdiff/tests/test_optimize.py @@ -3,7 +3,6 @@ """ # pylint: skip-file import numpy as np -from unittest import TestCase import pytest from pynumdiff.optimize.finite_difference import first_order @@ -16,7 +15,6 @@ from pynumdiff.utils import simulate - # simulation noise_type = 'normal' noise_parameters = [0, 0.01] @@ -34,133 +32,129 @@ def get_err_msg(actual_params, desired_params): err_msg = 'Actual params were: ' + ', '.join(map(str, actual_params)) + ' instead of: ' + ', '.join(map(str, desired_params)) return err_msg -class TestOPT(TestCase): - def test_first_order(self): - params_1, val_1 = first_order(x, dt, params=None, options={'iterate': True}, - tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = first_order(x, dt, params=None, options={'iterate': True}, - tvgamma=0, dxdt_truth=None) - self.assertListEqual(params_1, [5]) - self.assertListEqual(params_2, [1]) - - def test_mediandiff(self): - params_1, val_1 = mediandiff(x, dt, params=None, options={'iterate': False}, - tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = mediandiff(x, dt, params=None, options={'iterate': False}, - tvgamma=0, dxdt_truth=None) - self.assertListEqual(params_1, [5]) - self.assertListEqual(params_2, [1]) - - def test_meandiff(self): - params_1, val_1 = meandiff(x, dt, params=None, options={'iterate': False}, + +def test_first_order(): + params_1, val_1 = first_order(x, dt, params=None, options={'iterate': True}, + tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = first_order(x, dt, params=None, options={'iterate': True}, + tvgamma=0, dxdt_truth=None) + assert params_1 == [5] + assert params_2 == [1] + +def test_mediandiff(): + params_1, val_1 = mediandiff(x, dt, params=None, options={'iterate': False}, + tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = mediandiff(x, dt, params=None, options={'iterate': False}, + tvgamma=0, dxdt_truth=None) + assert params_1 == [5] + assert params_2 == [1] + +def test_meandiff(): + params_1, val_1 = meandiff(x, dt, params=None, options={'iterate': False}, + tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = meandiff(x, dt, params=None, options={'iterate': False}, + tvgamma=0, dxdt_truth=None) + assert params_1 == [5] + assert params_2 == [1] + +def test_gaussiandiff(): + params_1, val_1 = gaussiandiff(x, dt, params=None, options={'iterate': False}, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = meandiff(x, dt, params=None, options={'iterate': False}, + params_2, val_2 = gaussiandiff(x, dt, params=None, options={'iterate': False}, tvgamma=0, dxdt_truth=None) - self.assertListEqual(params_1, [5]) - self.assertListEqual(params_2, [1]) - - def test_gaussiandiff(self): - params_1, val_1 = gaussiandiff(x, dt, params=None, options={'iterate': False}, - tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = gaussiandiff(x, dt, params=None, options={'iterate': False}, - tvgamma=0, dxdt_truth=None) - self.assertListEqual(params_1, [9]) - self.assertListEqual(params_2, [1]) - - def test_friedrichsdiff(self): - params_1, val_1 = friedrichsdiff(x, dt, params=None, options={'iterate': False}, - tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = friedrichsdiff(x, dt, params=None, options={'iterate': False}, - tvgamma=0, dxdt_truth=None) - self.assertListEqual(params_1, [9]) - self.assertListEqual(params_2, [1]) - - def test_butterdiff(self): - params_1, val_1 = butterdiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = butterdiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) - - np.testing.assert_array_less( np.abs(params_1[0] - 9), 1.001, err_msg=get_err_msg(params_1, [9, 0.157])) - np.testing.assert_array_less( np.abs(params_1[1] - 0.157), 0.01, err_msg=get_err_msg(params_1, [9, 0.157])) - #np.testing.assert_almost_equal(params_1, [9, 0.157], decimal=3, err_msg=get_err_msg(params_1, [9, 0.157])) - np.testing.assert_almost_equal(params_2, [2, 0.99], decimal=3, err_msg=get_err_msg(params_2, [2, 0.99])) - - def test_splinediff(self): - params_1, val_1 = splinediff(x, dt, params=None, options={'iterate': True}, + assert params_1 == [9] + assert params_2 == [1] + +def test_friedrichsdiff(): + params_1, val_1 = friedrichsdiff(x, dt, params=None, options={'iterate': False}, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = splinediff(x, dt, params=None, options={'iterate': True}, + params_2, val_2 = friedrichsdiff(x, dt, params=None, options={'iterate': False}, tvgamma=0, dxdt_truth=None) - np.testing.assert_almost_equal(params_1, [5, 0.0147, 1], decimal=2) - np.testing.assert_almost_equal(params_2, [5, 0.0147, 1], decimal=2) + assert params_1 == [9] + assert params_2 == [1] + +def test_butterdiff(): + params_1, val_1 = butterdiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = butterdiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) + + np.testing.assert_array_less( np.abs(params_1[0] - 9), 1.001, err_msg=get_err_msg(params_1, [9, 0.157])) + np.testing.assert_array_less( np.abs(params_1[1] - 0.157), 0.01, err_msg=get_err_msg(params_1, [9, 0.157])) + #np.testing.assert_almost_equal(params_1, [9, 0.157], decimal=3, err_msg=get_err_msg(params_1, [9, 0.157])) + np.testing.assert_almost_equal(params_2, [2, 0.99], decimal=3, err_msg=get_err_msg(params_2, [2, 0.99])) + +def test_splinediff(): + params_1, val_1 = splinediff(x, dt, params=None, options={'iterate': True}, + tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = splinediff(x, dt, params=None, options={'iterate': True}, + tvgamma=0, dxdt_truth=None) + np.testing.assert_almost_equal(params_1, [5, 0.0147, 1], decimal=2) + np.testing.assert_almost_equal(params_2, [5, 0.0147, 1], decimal=2) + +def test_iterative_velocity(): + params_1, val_1 = iterative_velocity(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = iterative_velocity(x, dt, params=None, tvgamma=0, dxdt_truth=None) + np.testing.assert_array_less( np.abs(params_1[0] - 2), 1.001) + np.testing.assert_array_less( np.abs(params_2[0] - 2), 1.001) - def test_iterative_velocity(self): - params_1, val_1 = iterative_velocity(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = iterative_velocity(x, dt, params=None, tvgamma=0, dxdt_truth=None) - np.testing.assert_array_less( np.abs(params_1[0] - 2), 1.001) - np.testing.assert_array_less( np.abs(params_2[0] - 2), 1.001) - - np.testing.assert_almost_equal(params_1[1], 0.0001, decimal=4) - np.testing.assert_almost_equal(params_2[1], 0.0001, decimal=4) - - #self.assertListEqual(params_1, [2, 0.0001]) - #self.assertListEqual(params_2, [2, 0.0001]) + np.testing.assert_almost_equal(params_1[1], 0.0001, decimal=4) + np.testing.assert_almost_equal(params_2[1], 0.0001, decimal=4) - def test_velocity(self): - try: - import cvxpy - except: - pytest.skip("could not import cvxpy, skipping test_velocity", allow_module_level=True) + #self.assertListEqual(params_1, [2, 0.0001]) + #self.assertListEqual(params_2, [2, 0.0001]) - params_1, val_1 = velocity(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = velocity(x, dt, params=None, tvgamma=0, dxdt_truth=None) - param_1_error = np.abs(params_1[0] - 0.07218) - param_2_error = np.abs(params_2[0] - 0.0001) +def test_velocity(): + try: + import cvxpy + except: + pytest.skip("could not import cvxpy, skipping test_velocity", allow_module_level=True) - np.testing.assert_array_less(param_1_error, 2) - np.testing.assert_array_less(param_2_error, 2) + params_1, val_1 = velocity(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = velocity(x, dt, params=None, tvgamma=0, dxdt_truth=None) + param_1_error = np.abs(params_1[0] - 0.07218) + param_2_error = np.abs(params_2[0] - 0.0001) - return - - def test_acceleration(self): - try: - import cvxpy - except: - pytest.skip("could not import cvxpy, skipping test_acceleration", allow_module_level=True) + np.testing.assert_array_less(param_1_error, 2) + np.testing.assert_array_less(param_2_error, 2) - params_1, val_1 = acceleration(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = acceleration(x, dt, params=None, tvgamma=0, dxdt_truth=None) - param_1_error = np.abs(params_1[0] - 0.1447) - param_2_error = np.abs(params_2[0] - 0.0001) +def test_acceleration(): + try: + import cvxpy + except: + pytest.skip("could not import cvxpy, skipping test_acceleration", allow_module_level=True) - np.testing.assert_array_less(param_1_error, 2) - np.testing.assert_array_less(param_2_error, 2) + params_1, val_1 = acceleration(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = acceleration(x, dt, params=None, tvgamma=0, dxdt_truth=None) + param_1_error = np.abs(params_1[0] - 0.1447) + param_2_error = np.abs(params_2[0] - 0.0001) - return - - def test_savgoldiff(self): - params_1, val_1 = savgoldiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = savgoldiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) - self.assertListEqual(params_1, [10, 59, 3]) - self.assertListEqual(params_2, [9, 3, 3]) - - def test_spectraldiff(self): - params_1, val_1 = spectraldiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = spectraldiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) - np.testing.assert_almost_equal(params_1, [0.0912], decimal=3) - np.testing.assert_almost_equal(params_2, [0.575], decimal=3) - - def test_polydiff(self): - params_1, val_1 = polydiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - params_2, val_2 = polydiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) - self.assertListEqual(params_1, [6, 50]) - self.assertListEqual(params_2, [4, 10]) - - # def test_chebydiff(self): - # try: - # import pychebfun - # except: - # pytest.skip("could not import pychebfun, skipping test_chebydiff", allow_module_level=True) - - # params_1, val_1 = chebydiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) - # params_2, val_2 = chebydiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) - # self.assertListEqual(params_1, [9, 108]) - # self.assertListEqual(params_2, [9, 94]) + np.testing.assert_array_less(param_1_error, 2) + np.testing.assert_array_less(param_2_error, 2) + +def test_savgoldiff(): + params_1, val_1 = savgoldiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = savgoldiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) + assert params_1 == [10, 59, 3] + assert params_2 == [9, 3, 3] + +def test_spectraldiff(): + params_1, val_1 = spectraldiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = spectraldiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) + np.testing.assert_almost_equal(params_1, [0.0912], decimal=3) + np.testing.assert_almost_equal(params_2, [0.575], decimal=3) + +def test_polydiff(): + params_1, val_1 = polydiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) + params_2, val_2 = polydiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) + assert params_1 == [6, 50] + assert params_2 == [4, 10] + +# def test_chebydiff(self): +# try: +# import pychebfun +# except: +# pytest.skip("could not import pychebfun, skipping test_chebydiff", allow_module_level=True) + +# params_1, val_1 = chebydiff(x, dt, params=None, tvgamma=tvgamma, dxdt_truth=dxdt_truth) +# params_2, val_2 = chebydiff(x, dt, params=None, tvgamma=0, dxdt_truth=None) +# self.assertListEqual(params_1, [9, 108]) +# self.assertListEqual(params_2, [9, 94]) diff --git a/pynumdiff/tests/test_smooth_finite_difference.py b/pynumdiff/tests/test_smooth_finite_difference.py index 403d79f..5429990 100644 --- a/pynumdiff/tests/test_smooth_finite_difference.py +++ b/pynumdiff/tests/test_smooth_finite_difference.py @@ -3,7 +3,6 @@ """ # pylint: skip-file import numpy as np -from unittest import TestCase from pynumdiff.smooth_finite_difference import mediandiff, meandiff, gaussiandiff, \ friedrichsdiff, butterdiff, splinediff @@ -14,69 +13,66 @@ dt = 0.01 -class TestSFD(TestCase): - def test_mediandiff(self): - 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_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(self): - 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_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(self): - 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_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(self): - 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(self): - 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(self): - 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) +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) diff --git a/pynumdiff/tests/test_total_variation_regularization.py b/pynumdiff/tests/test_total_variation_regularization.py index e38728f..61f4188 100644 --- a/pynumdiff/tests/test_total_variation_regularization.py +++ b/pynumdiff/tests/test_total_variation_regularization.py @@ -4,7 +4,6 @@ # pylint: skip-file import numpy as np -from unittest import TestCase import pytest from pynumdiff.total_variation_regularization import * @@ -14,104 +13,103 @@ dt = 0.01 -class TestTVR(TestCase): - def test_velocity(self): - try: - import cvxpy - except: - pytest.skip("could not import cvxpy, skipping test_velocity", allow_module_level=True) +def test_velocity(): + try: + import cvxpy + except: + pytest.skip("could not import cvxpy, skipping test_velocity", allow_module_level=True) - params = [0.5] - x_hat, dxdt_hat = velocity(x, dt, params, options={'solver': 'CVXOPT'}) - x_smooth = np.array([1.60206974, 3.84254116, 6.08301239, 8.32348272, 14.76608638, - 12.76589239, 10.76569864, 7.70248886, 11.43239643, 11.4325017, - 11.43260691, 6.42354819, 5.78305309, 5.14255819, 4.50206322]) - dxdt = np.array([2.24047187e+02, 2.24047133e+02, 2.24047078e+02, 4.34153700e+02, - 2.22120483e+02, -2.00019387e+02, -2.53170177e+02, 3.33348898e+01, - 1.86500642e+02, 1.05238579e-02, -2.50447675e+02, -2.82477691e+02, - -6.40494998e+01, -6.40494935e+01, -6.40494871e+01]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) + params = [0.5] + x_hat, dxdt_hat = velocity(x, dt, params, options={'solver': 'CVXOPT'}) + x_smooth = np.array([1.60206974, 3.84254116, 6.08301239, 8.32348272, 14.76608638, + 12.76589239, 10.76569864, 7.70248886, 11.43239643, 11.4325017, + 11.43260691, 6.42354819, 5.78305309, 5.14255819, 4.50206322]) + dxdt = np.array([2.24047187e+02, 2.24047133e+02, 2.24047078e+02, 4.34153700e+02, + 2.22120483e+02, -2.00019387e+02, -2.53170177e+02, 3.33348898e+01, + 1.86500642e+02, 1.05238579e-02, -2.50447675e+02, -2.82477691e+02, + -6.40494998e+01, -6.40494935e+01, -6.40494871e+01]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - def test_iterative_velocity(self): - params = [1, 0.05] - x_hat, dxdt_hat = iterative_velocity(x, dt, params) - x_smooth = np.array([1.256, 3.254, 5.249, 7.197, 8.96, 10.287, 11.08, 11.407, - 11.305, 10.875, 10.235, 9.371, 8.305, 7.174, 6.042]) - dxdt = np.array([199.802, 199.742, 199.222, 190.43, 162.105, 103.282, - 55.311, 10.12, -30.571, -55.409, -72.603, -100.119, - -113.097, -113.097, -113.464]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) +def test_iterative_velocity(): + params = [1, 0.05] + x_hat, dxdt_hat = iterative_velocity(x, dt, params) + x_smooth = np.array([1.256, 3.254, 5.249, 7.197, 8.96, 10.287, 11.08, 11.407, + 11.305, 10.875, 10.235, 9.371, 8.305, 7.174, 6.042]) + dxdt = np.array([199.802, 199.742, 199.222, 190.43, 162.105, 103.282, + 55.311, 10.12, -30.571, -55.409, -72.603, -100.119, + -113.097, -113.097, -113.464]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - def test_acceleration(self): - try: - import cvxpy - except: - pytest.skip("could not import cvxpy, skipping test_acceleration", allow_module_level=True) +def test_acceleration(): + try: + import cvxpy + except: + pytest.skip("could not import cvxpy, skipping test_acceleration", allow_module_level=True) - params = [1] - x_hat, dxdt_hat = acceleration(x, dt, params, options={'solver': 'CVXOPT'}) - x_smooth = np.array([0.87728375, 4.44441238, 7.31141687, 9.47829719, 10.94505335, - 11.7116852, 11.78131319, 11.5560333, 11.03584752, 10.2207553, - 9.11075633, 7.7058506, 6.41426253, 5.23599238, 4.17104012]) - dxdt = np.array([391.71907211, 321.70665613, 251.69424015, 181.6818242, - 111.66940057, 41.81299196, -7.78259499, -37.27328368, - -66.76389967, -96.25455924, -125.74523529, -134.82469003, - -123.49291116, -112.16112081, -100.82933046]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) + params = [1] + x_hat, dxdt_hat = acceleration(x, dt, params, options={'solver': 'CVXOPT'}) + x_smooth = np.array([0.87728375, 4.44441238, 7.31141687, 9.47829719, 10.94505335, + 11.7116852, 11.78131319, 11.5560333, 11.03584752, 10.2207553, + 9.11075633, 7.7058506, 6.41426253, 5.23599238, 4.17104012]) + dxdt = np.array([391.71907211, 321.70665613, 251.69424015, 181.6818242, + 111.66940057, 41.81299196, -7.78259499, -37.27328368, + -66.76389967, -96.25455924, -125.74523529, -134.82469003, + -123.49291116, -112.16112081, -100.82933046]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - def test_smooth_acceleration(self): - try: - import cvxpy - except: - pytest.skip("could not import cvxpy, skipping test_smooth_acceleration", allow_module_level=True) +def test_smooth_acceleration(): + try: + import cvxpy + except: + pytest.skip("could not import cvxpy, skipping test_smooth_acceleration", allow_module_level=True) - params = [5, 30] - x_hat, dxdt_hat = smooth_acceleration(x, dt, params, options={'solver': 'CVXOPT'}) - x_smooth = np.array([4.16480747, 5.52913444, 6.77037146, 7.87267273, 8.79483088, - 9.5044844, 9.97926076, 10.20730827, 10.18728338, 9.92792114, - 9.44728533, 8.77174094, 7.93472066, 6.97538656, 5.93725369]) - dxdt = np.array([136.43269721, 129.9388182, 118.30858578, 102.15166804, - 82.27996127, 59.65074227, 35.30453082, 10.30497111, - -14.30994982, -37.56249817, -58.56466324, -76.54421499, - -90.85984169, -101.00697716, -106.61959829]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) + params = [5, 30] + x_hat, dxdt_hat = smooth_acceleration(x, dt, params, options={'solver': 'CVXOPT'}) + x_smooth = np.array([4.16480747, 5.52913444, 6.77037146, 7.87267273, 8.79483088, + 9.5044844, 9.97926076, 10.20730827, 10.18728338, 9.92792114, + 9.44728533, 8.77174094, 7.93472066, 6.97538656, 5.93725369]) + dxdt = np.array([136.43269721, 129.9388182, 118.30858578, 102.15166804, + 82.27996127, 59.65074227, 35.30453082, 10.30497111, + -14.30994982, -37.56249817, -58.56466324, -76.54421499, + -90.85984169, -101.00697716, -106.61959829]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - def test_jerk(self): - try: - import cvxpy - except: - pytest.skip("could not import cvxpy, skipping test_jerk", allow_module_level=True) +def test_jerk(): + try: + import cvxpy + except: + pytest.skip("could not import cvxpy, skipping test_jerk", allow_module_level=True) - params = [10] - x_hat, dxdt_hat = jerk(x, dt, params, options={'solver': 'CVXOPT'}) - x_smooth = np.array([0.71013072, 4.51405229, 7.42619407, 9.53278029, 10.92003519, - 11.674183, 11.88144796, 11.6280543, 11.00022625, 10.08418804, - 8.9661639, 7.73237808, 6.4690548, 5.2624183, 4.19869281]) - dxdt = np.array([420.66993476, 335.80316742, 250.93640008, 174.69205619, - 107.07013572, 48.07063861, -2.30643522, -44.06108581, - -77.19331317, -101.70311726, -117.59049798, -124.85545525, - -123.49798898, -113.51809914, -103.5382093]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) + params = [10] + x_hat, dxdt_hat = jerk(x, dt, params, options={'solver': 'CVXOPT'}) + x_smooth = np.array([0.71013072, 4.51405229, 7.42619407, 9.53278029, 10.92003519, + 11.674183, 11.88144796, 11.6280543, 11.00022625, 10.08418804, + 8.9661639, 7.73237808, 6.4690548, 5.2624183, 4.19869281]) + dxdt = np.array([420.66993476, 335.80316742, 250.93640008, 174.69205619, + 107.07013572, 48.07063861, -2.30643522, -44.06108581, + -77.19331317, -101.70311726, -117.59049798, -124.85545525, + -123.49798898, -113.51809914, -103.5382093]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) - def test_jerk_sliding(self): - try: - import cvxpy - except: - pytest.skip("could not import cvxpy, skipping test_jerk_sliding", allow_module_level=True) +def test_jerk_sliding(): + try: + import cvxpy + except: + pytest.skip("could not import cvxpy, skipping test_jerk_sliding", allow_module_level=True) - params = [10] - x_hat, dxdt_hat = jerk_sliding(x, dt, params, options={'solver': 'CVXOPT'}) - x_smooth = np.array([0.71013072, 4.51405229, 7.42619407, 9.53278029, 10.92003519, - 11.674183, 11.88144796, 11.6280543, 11.00022625, 10.08418804, - 8.9661639, 7.73237808, 6.4690548, 5.2624183, 4.19869281]) - dxdt = np.array([420.66993476, 335.80316742, 250.93640008, 174.69205619, - 107.07013572, 48.07063861, -2.30643522, -44.06108581, - -77.19331317, -101.70311726, -117.59049798, -124.85545525, - -123.49798898, -113.51809914, -103.5382093]) - np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) - np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) + params = [10] + x_hat, dxdt_hat = jerk_sliding(x, dt, params, options={'solver': 'CVXOPT'}) + x_smooth = np.array([0.71013072, 4.51405229, 7.42619407, 9.53278029, 10.92003519, + 11.674183, 11.88144796, 11.6280543, 11.00022625, 10.08418804, + 8.9661639, 7.73237808, 6.4690548, 5.2624183, 4.19869281]) + dxdt = np.array([420.66993476, 335.80316742, 250.93640008, 174.69205619, + 107.07013572, 48.07063861, -2.30643522, -44.06108581, + -77.19331317, -101.70311726, -117.59049798, -124.85545525, + -123.49798898, -113.51809914, -103.5382093]) + np.testing.assert_almost_equal(x_smooth, x_hat, decimal=2) + np.testing.assert_almost_equal(dxdt, dxdt_hat, decimal=2) diff --git a/pynumdiff/tests/test_utils.py b/pynumdiff/tests/test_utils.py index 6f6df14..c5ad026 100644 --- a/pynumdiff/tests/test_utils.py +++ b/pynumdiff/tests/test_utils.py @@ -4,17 +4,14 @@ # pylint: skip-file import numpy as np -from unittest import TestCase from pynumdiff.utils import utility, simulate, evaluate -class TestKS(TestCase): - def test_utility(self): - return +def test_utility(): + return# TODO. There are a lot of basic integration functionalities, etc. that deserve tests. - def test_simulate(self): - return - - def test_evaluate(self): - return +def test_simulate(): + return +def test_evaluate(): + return diff --git a/pynumdiff/total_variation_regularization/__init__.py b/pynumdiff/total_variation_regularization/__init__.py index 93c4b77..4729519 100644 --- a/pynumdiff/total_variation_regularization/__init__.py +++ b/pynumdiff/total_variation_regularization/__init__.py @@ -1,36 +1,16 @@ +"""This module implements some common total variation regularization methods """ -Import useful functions from _total_variation_regularization module -""" -import logging as _logging -_logging.basicConfig( - level=_logging.INFO, - format="%(asctime)s [%(levelname)s] %(message)s", - handlers=[ - _logging.FileHandler("debug.log"), - _logging.StreamHandler() - ] -) - try: import cvxpy - __cvxpy_installed__ = True + from ._total_variation_regularization import velocity, acceleration, jerk, jerk_sliding, smooth_acceleration except: - __cvxpy_installed__ = False + from warnings import warn + warn("""Limited Total Variation Regularization Support Detected! CVXPY is not installed. + Many Total Variation Methods require CVXPY including: velocity, acceleration, jerk, jerk_sliding, smooth_acceleration. + Please install CVXPY to use these methods. Recommended to also install MOSEK and obtain a MOSEK license. + You can still use: total_variation_regularization.iterative_velocity.""") + +from ._total_variation_regularization import iterative_velocity -from pynumdiff.total_variation_regularization._total_variation_regularization import iterative_velocity +__all__ = ['velocity', 'acceleration', 'jerk', 'jerk_sliding', 'smooth_acceleration', 'iterative_velocity'] # So these get treated as direct members of the module by sphinx -if __cvxpy_installed__: - from pynumdiff.total_variation_regularization._total_variation_regularization import velocity - from pynumdiff.total_variation_regularization._total_variation_regularization import acceleration - from pynumdiff.total_variation_regularization._total_variation_regularization import jerk - from pynumdiff.total_variation_regularization._total_variation_regularization import jerk_sliding - from pynumdiff.total_variation_regularization._total_variation_regularization import smooth_acceleration -else: - __warning__ = '\nLimited Total Variation Regularization Support Detected! \n'\ - '---> CVXPY is not installed. \n'\ - '---> Many Total Variation Methods require CVXPY including: \n'\ - '---> velocity, acceleration, jerk, jerk_sliding, smooth_acceleration\n'\ - '---> Please install CVXPY to use these methods.\n'\ - '---> Recommended to also install MOSEK and obtain a MOSEK license.\n'\ - 'You can still use: total_variation_regularization.iterative_velocity\n' - _logging.info(__warning__) diff --git a/pynumdiff/total_variation_regularization/_total_variation_regularization.py b/pynumdiff/total_variation_regularization/_total_variation_regularization.py index 789f89d..04bff1a 100644 --- a/pynumdiff/total_variation_regularization/_total_variation_regularization.py +++ b/pynumdiff/total_variation_regularization/_total_variation_regularization.py @@ -1,13 +1,9 @@ -""" -This module implements some common total variation regularization methods -""" import logging import numpy as np from pynumdiff.total_variation_regularization import __chartrand_tvregdiff__ import pynumdiff.smooth_finite_difference from pynumdiff.utils import utility -__gaussian_kernel__ = utility.__gaussian_kernel__ try: import cvxpy @@ -283,8 +279,8 @@ def smooth_acceleration(x, dt, params, options=None): gamma, window_size = params x_hat, dxdt_hat = acceleration(x, dt, [gamma], options=options) - kernel = __gaussian_kernel__(window_size) - dxdt_hat = pynumdiff.smooth_finite_difference.__convolutional_smoother__(dxdt_hat, kernel, 1) + 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_initial_condition(x, x_hat) diff --git a/pynumdiff/utils/__init__.py b/pynumdiff/utils/__init__.py index c1e37aa..e69de29 100644 --- a/pynumdiff/utils/__init__.py +++ b/pynumdiff/utils/__init__.py @@ -1,6 +0,0 @@ -""" -Import useful functions from utils -""" -from pynumdiff.utils import utility -from pynumdiff.utils import simulate -from pynumdiff.utils import evaluate diff --git a/pynumdiff/utils/evaluate.py b/pynumdiff/utils/evaluate.py index a6a1f32..9ee3ee8 100644 --- a/pynumdiff/utils/evaluate.py +++ b/pynumdiff/utils/evaluate.py @@ -7,7 +7,7 @@ # local imports from pynumdiff.utils import utility as _utility -_finite_difference = _utility.finite_difference +from pynumdiff.finite_difference import first_order as _finite_difference # pylint: disable-msg=too-many-locals, too-many-arguments diff --git a/pynumdiff/utils/simulate.py b/pynumdiff/utils/simulate.py index b8f3f17..5a389d8 100644 --- a/pynumdiff/utils/simulate.py +++ b/pynumdiff/utils/simulate.py @@ -5,9 +5,9 @@ from scipy.integrate import odeint # local imports -from pynumdiff.utils import utility as _utility +from pynumdiff.utils.utility import peakdet from pynumdiff.utils import _pi_cruise_control -_finite_difference = _utility.finite_difference +from pynumdiff.finite_difference import first_order as _finite_difference # pylint: disable-msg=too-many-locals, too-many-arguments, no-member @@ -136,7 +136,7 @@ def triangle(timeseries_length=4, noise_type='normal', noise_parameters=(0, 0.5) continuous_x = _np.sin(t*t) # find peaks and valleys - peaks, valleys = _utility.peakdet(continuous_x, 0.1) + peaks, valleys = peakdet(continuous_x, 0.1) # organize peaks and valleys if len(peaks) > 0: diff --git a/pynumdiff/utils/utility.py b/pynumdiff/utils/utility.py index c561e71..e43c260 100644 --- a/pynumdiff/utils/utility.py +++ b/pynumdiff/utils/utility.py @@ -1,69 +1,5 @@ -""" -All kinds of utilities -""" -import os -import sys -import copy +import os, sys, copy, scipy import numpy as np -import scipy - - -def get_filenames(path, contains, does_not_contain=('~', '.pyc')): - """ - Create list of files found in given path that contain or do not contain certain strings. - - :param path: path in which to look for files - :type path: string (directory path) - - :param contains: string that filenames must contain - :type contains: string - - :param does_not_contain: list of strings that filenames must not contain, optional - :type does_not_contain: list of strings - - :return: list of filenames - :rtype: list of strings - """ - cmd = 'ls ' + '"' + path + '"' - ls = os.popen(cmd).read() - all_filelist = ls.split('\n') - all_filelist.remove('') - filelist = [] - for _, filename in enumerate(all_filelist): - if contains in filename: - fileok = True - for nc in does_not_contain: - if nc in filename: - fileok = False - if fileok: - filelist.append(os.path.join(path, filename)) - return filelist - - -def is_odd(num): - """ - tell if it is an odd number - - :param num: (integer) the number we need to tell - :return: (boolean) true or false - """ - return num & 0x1 - - -def isnotebook(): - """ - Checks to see if the environment is an interactive notebook or not. - :return: True or False - """ - try: - shell = get_ipython().__class__.__name__ - if shell == 'ZMQInteractiveShell': - return True # Jupyter notebook or qtconsole - # elif shell == 'TerminalInteractiveShell': - # return False # Terminal running IPython - return False # Other type (?) - except NameError: - return False # Probably standard Python interpreter def hankel_matrix(x, num_delays, pad=False): # fixed delay step of 1 @@ -84,7 +20,6 @@ def hankel_matrix(x, num_delays, pad=False): # fixed delay step of 1 ['b', 'c', 'd', 'e', 0], ['c', 'd', 'e', 0, 0]] """ - m = copy.copy(x) for d in range(1, num_delays): xi = x[:, d:] @@ -191,43 +126,18 @@ def peakdet(v, delta, x=None): return np.array(maxtab), np.array(mintab) -# simple finite difference -def finite_difference(x, dt): - """ - :param x: (np.array of floats, 1xN) time series to differentiate - :param dt: (float) time step - :param params: (list): (int, optional) number of iterations ignored if 'iterate' not in options - :param options: (bool) iterate the finite difference method (smooths the estimates) - :return: x_hat: smoothed x; dxdt_hat: derivative of x - """ - # Calculate the finite difference - dxdt_hat = np.diff(x)/dt - # Pad the data - dxdt_hat = np.hstack((dxdt_hat[0], dxdt_hat, dxdt_hat[-1])) - # Re-finite dxdt_hat using linear interpolation - dxdt_hat = np.mean((dxdt_hat[0:-1], dxdt_hat[1:]), axis=0) - - return x, dxdt_hat - - # Trapazoidal integration, with interpolated final point so that the lengths match. def integrate_dxdt_hat(dxdt_hat, dt): - """ - Wrapper for scipy.integrate.cumulative_trapezoid to integrate dxdt_hat that ensures the integral has the same length + """Wrapper for scipy.integrate.cumulative_trapezoid to integrate dxdt_hat that ensures the integral has the same length - :param dxdt_hat: estimate derivative of timeseries - :type dxdt_hat: np.array + :param np.array[float] dxdt_hat: estimate derivative of timeseries + :param float dt: time step in seconds - :param dt: time step in seconds - :type dt: float - - :return: integral of dxdt_hat - :rtype: np.array + :return: **x_hat** (np.array[float]) -- integral of dxdt_hat """ x = scipy.integrate.cumulative_trapezoid(dxdt_hat) - first_value = x[0] - np.mean(dxdt_hat[0:1]) - x = np.hstack((first_value, x))*dt - return x + first_value = x[0] - dxdt_hat[0] + return np.hstack((first_value, x))*dt # Optimization routine to estimate the integration constant. @@ -253,31 +163,44 @@ def f(x0, *args): # kernels -def __mean_kernel__(window_size): - """ - :param window_size: - :return: +def _mean_kernel(window_size): + """A uniform boxcar of total integral 1 """ return np.ones(window_size)/window_size -def __gaussian_kernel__(window_size): - """ - :param window_size: - :return: +def _gaussian_kernel(window_size): + """A truncated gaussian """ sigma = window_size / 6. t = np.linspace(-2.7*sigma, 2.7*sigma, window_size) - gaussian_func = lambda t, sigma: 1/np.sqrt(2*np.pi*sigma**2) * np.exp(-(t**2)/(2*sigma**2)) - ker = gaussian_func(t, sigma) + ker = 1/np.sqrt(2*np.pi*sigma**2) * np.exp(-(t**2)/(2*sigma**2)) # gaussian function itself return ker / np.sum(ker) -def __friedrichs_kernel__(window_size): - """ - :param window_size: - :return: +def _friedrichs_kernel(window_size): + """A bump function """ x = np.linspace(-0.999, 0.999, window_size) ker = np.exp(-1/(1-x**2)) return ker / np.sum(ker) + + +def convolutional_smoother(x, kernel, iterations): + """Perform smoothing by convolving x with a kernel. + + :param np.array[float] x: 1D data + :param np.array[float] kernel: kernel to use in convolution + :param int iterations: number of iterations, >=1 + :return: **x_hat** (np.array[float]) -- smoothed x + """ + x_hat = np.hstack((x[::-1], x, x[::-1])) # pad + w = np.arange(len(x_hat)) / (len(x_hat) - 1) # weights + + for _ in range(iterations): + x_hat_f = np.convolve(x_hat, kernel, 'same') + x_hat_b = np.convolve(x_hat[::-1], kernel, 'same')[::-1] + + x_hat = x_hat_f*w + x_hat_b*(1-w) + + return x_hat[len(x):len(x)*2] diff --git a/pyproject.toml b/pyproject.toml index 11dc537..a94d10c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,6 +11,7 @@ license = {text = "MIT"} maintainers = [ {name = "Floris van Breugel", email = "fvanbreugel@unr.edu"}, {name = "Yuying Liu", email = "yliu814@uw.edu"}, + {name = "Pavel Komarov", email = "pvlkmrv@uw.edu"} ] keywords = ["derivative", "smoothing", "curve fitting", "optimization", "total variation"] classifiers = [ @@ -25,7 +26,7 @@ classifiers = [ dependencies = [ "numpy", "scipy", - "matplotlib", + "matplotlib" ] [project.urls] @@ -39,81 +40,11 @@ advanced = [ "MOSEK", ] dev = [ -"appnope==0.1.2", -"argon2-cffi==20.1.0", -"astroid==2.6.2", -"async-generator==1.10", -"attrs==21.2.0", -"backcall==0.2.0", -"bleach==3.3.0", -"cffi~=1.14", -"cycler==0.10.0", -"debugpy==1.3.0", -"decorator==5.0.9", -"defusedxml==0.7.1", -"entrypoints==0.3", -"iniconfig==1.1.1", -"ipykernel==6.0.1", -"ipython==8.10.0", -"ipython-genutils==0.2.0", -"isort==5.9.1", -"jedi==0.18.0", -"Jinja2==3.0.1", -"jsonschema==3.2.0", -"jupyter-client==6.1.12", -"jupyter-core==4.11.2", -"jupyterlab-pygments==0.1.2", -"kiwisolver==1.3.1", -"lazy-object-proxy==1.6.0", -"MarkupSafe==2.0.1", -"matplotlib==3.4.2", -"matplotlib-inline==0.1.2", -"mccabe==0.6.1", -"mistune==0.8.4", -"Mosek>=9.2.47", -"nbclient==0.5.3", -"nbconvert==6.5.1", -"nbformat==5.1.3", -"nest-asyncio==1.5.1", -"notebook==6.4.12", -"numpy==1.24.0", -"osqp~=0.6.2.post0", -"packaging==21.0", -"pandocfilters==1.4.3", -"parso==0.8.2", -"pexpect==4.8.0", -"pickleshare==0.7.5", -"Pillow==9.3.0", -"pluggy==0.13.1", -"prometheus-client==0.11.0", -"prompt-toolkit==3.0.38", -"ptyprocess==0.7.0", -"py==1.10.0", -"pycparser==2.20", -"Pygments==2.9.0", -"pylint==2.9.3", -"pyparsing==2.4.7", -"pyrsistent==0.18.0", -"pytest==6.2.4", -"python-dateutil==2.8.1", -"pyzmq==22.1.0", -"qdldl==0.1.5.post0", -"scipy==1.10.0", -"Send2Trash>=1.7.1", -"six==1.16.0", -"terminado==0.10.1", -"testpath==0.5.0", -"toml==0.10.2", -"tornado==6.1", -"traitlets==5.0.5", -"wcwidth==0.2.5", -"webencodings==0.5.1", -"ecos==2.0.12", -"scs==3.2.3", -"wrapt==1.12.1", -"cvxopt==1.3.0", -"cvxpy==1.3.1", -"setuptools_scm" + "pylint", + "pytest", + "cvxopt", + "cvxpy", + "Mosek" ] [tool.setuptools_scm]