diff --git a/AUTHORS.rst b/AUTHORS.rst index 180ebd3..f81966b 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -3,4 +3,4 @@ Developers ========== * Per A. Brodtkorb -* John D'Errico \ No newline at end of file +* John D'Errico diff --git a/README.rst b/README.rst index d50e582..4ceade2 100644 --- a/README.rst +++ b/README.rst @@ -17,12 +17,12 @@ numerical differentiation problems in one or more variables. Finite differences are used in an adaptive manner, coupled with a Romberg extrapolation methodology to provide a maximally accurate result. The user can configure many options like; changing the order of the method or -the extrapolation, even allowing the user to specify whether central, forward or +the extrapolation, even allowing the user to specify whether complex-step, central, forward or backward differences are used. The methods provided are: -- **Derivative**: Compute the derivatives of order 1 through 4 on any scalar function. +- **Derivative**: Compute the derivatives of order 1 through 10 on any scalar function. - **Gradient**: Compute the gradient vector of a scalar function of one or more variables. @@ -36,12 +36,19 @@ All of these methods also produce error estimates on the result. The documentation for numdifftools is available here http://numdifftools.readthedocs.org/ + Download the toolbox here: http://pypi.python.org/pypi/Numdifftools ---- News """" +2015 +---- +August 20 +^^^^^^^^^ +`New release of Numdifftools 0.9.0. `_ + 2014 ---- December 18 diff --git a/docs/index.rst b/docs/index.rst index 63b759f..54b9e6e 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -14,7 +14,7 @@ Contents Overview - Theory + Numerical differentiation License Module Reference <_rst/modules> diff --git a/docs/src/derivest.rst b/docs/src/derivest.rst index 65e9122..2a0fa01 100644 --- a/docs/src/derivest.rst +++ b/docs/src/derivest.rst @@ -1,5 +1,5 @@ -Introduction - Derivative Estimation -#################################### +Introduction +############ The general problem of differentiation of a function typically pops up in three ways in Python. @@ -26,13 +26,11 @@ Surely you recall the traditional definition of a derivative, in terms of a limi For small :math:`\delta`, the limit approaches :math:`f'(x)`. This is a one-sided approximation for the derivative. For a fixed value of :math:`\delta`, this is also known as a finite difference approximation (a forward difference.) Other approximations for the derivative are also available. We will see the origin of these approximations in the Taylor series expansion of a function :math:`f(x)` around some point :math:`x_0`. .. math:: - f(x_0+\delta) &= f(x_0) + \delta f'(x_0) + \frac{\delta^2}{2} f''(x_0) + + f(x_0+\delta) &= f(x_0) + \delta f'(x_0) + \frac{\delta^2}{2} f''(x_0) + \frac{\delta^3}{6} f^{(3)}(x_0) + \\ :label: 2 - & \frac{\delta^3}{6} f^{(3)}(x_0) + \frac{\delta^4}{24} f^{(4)}(x_0) + \\ - + & \frac{\delta^4}{24} f^{(4)}(x_0) + \frac{\delta^5}{120} f^{(5)}(x_0) + \frac{\delta^6}{720} f^{(6)}(x_0) +...\\ - & \frac{\delta^5}{120} f^{(5)}(x_0) + \frac{\delta^6}{720} f^{(6)}(x_0) +...\\ Truncate the series in :eq:`2` to the first three terms, divide by :math:`\delta` and rearrange yields the forward difference approximation :eq:`1`: @@ -69,24 +67,20 @@ Returning to the Taylor series expansion of :math:`f(x)` around some point :math f_{even}(x) = \frac{f(x_0 + x) - 2f(x_0) + f(x_0 - x)}{2} :label: 6 -The Taylor series expansion of :math:`f_{odd}(x-x_0)` has the useful property that we have killed off any even order terms, but the odd order terms are identical to :math:`f(x)`, as expanded around :math:`x_0`. +The Taylor series expansion of :math:`f_{odd}(x)` around zero has the useful property that we have killed off any even order terms, but the odd order terms are identical to :math:`f(x)`, as expanded around :math:`x_0`. .. math:: - f_{odd}(x-x_0) &= (x - x_0)f'(x_0) + \frac{(x - x_0)^3}{6} f^{(3)}(x_0) + + f_{odd}(\delta) = \delta f'(x_0) + \frac{\delta^3}{6} f^{(3)}(x_0) + \frac{\delta^5}{120} f^{(5)}(x_0) + \frac{\delta^7}{5040} f^{(7)}(x_0) +... :label: 7 - & \frac{(x - x_0)^5}{120} f^{(5)}(x_0) + \frac{(x - x_0)^7}{5040} f^{(7)}(x_0) +...\\ - Likewise, the Taylor series expansion of :math:`f_{even}(x)` has no odd order terms or a constant term, but other even order terms that are identical to :math:`f(x)`. .. index:: even transformation .. math:: - f_{even}(x-x_0) &= \frac{(x - x_0)^2}{2} f^{(2)}(x_0) + \frac{(x - x_0)^4}{24} f^{(4)}(x_0) + + f_{even}(\delta) = \frac{\delta^2}{2} f^{(2)}(x_0) + \frac{\delta^4}{24} f^{(4)}(x_0) + \frac{\delta^6}{720} f^{(6)}(x_0) + \frac{\delta^8}{40320} f^{(8)}(x_0) + ... :label: 8 - & \frac{(x - x_0)^6}{720} f^{(6)}(x_0) + \frac{(x - x_0)^8}{40320} f^{(8)}(x_0) + ...\\ - The point of these transformations is we can rather simply generate a higher order approximation for any odd order derivatives of :math:`f(x)` by working with :math:`f_{odd}(x)`. Even order derivatives of :math:`f(x)` are similarly generated from :math:`f_{even}(x)`. For example, a second order approximation for :math:`f'(x_0)` is trivially written in :eq:`9` as a function of :math:`\delta`. @@ -117,55 +111,64 @@ with a complex step :math:`i h`: f(x_0+ i h) &= f(x_0) + i h f'(x_0) - \frac{h^2}{2} f''(x_0) - \frac{i h^3}{6} f^{(3)}(x_0) + \frac{h^4}{24} f^{(4)}(x_0) + \\ :label: 12a - & \frac{i h^5}{120} f^{(5)}(x_0) - \frac{h^6}{720} f^{(6)}(x_0) +...\\ + & \frac{i h^5}{120} f^{(5)}(x_0) - \frac{h^6}{720} f^{(6)}(x_0) -...\\ Taking only the imaginary parts of both sides gives .. math:: - \Im(f(x_0+ i h)) &= h f'(x_0) - \frac{i h^3}{6} f^{(3)}(x_0) - ... + \, \Im \,(f(x_0+ i h)) &= h f'(x_0) - \frac{h^3}{6} f^{(3)}(x_0) + \frac{h^5}{120} f^{(5)}(x_0) -... :label: 12b Dividing with :math:`h` and rearranging yields: .. math:: - f'(x_0) = \Im(f(x_0+ i h))/ h + \frac{h^2}{6} f^{(3)}(x_0) - ... + f'(x_0) = \Im(f(x_0+ i h))/ h + \frac{h^2}{6} f^{(3)}(x_0) - \frac{h^4}{120} f^{(5)}(x_0) +... :label: 12c Terms with order :math:`h^2` or higher can safely be ignored since the interval :math:`h` can be chosen up to machine precision -without fear of rounding errors stemming from subtraction (since there are not any). Thus to within first-order the complex-step derivative approximation is given by: +without fear of rounding errors stemming from subtraction (since there are not any). Thus to within second-order the complex-step derivative approximation is given by: .. math:: f'(x_0) = \Im(f(x_0 + i h))/ h :label: 12d -Next, consider replacing the step :math:`\delta` in :Eq:`8` with a complex step :math:`i^\frac{1}{2} h` +Next, consider replacing the step :math:`\delta` in :Eq:`8` with the complex step :math:`i^\frac{1}{2} h`: .. math:: - f_{even}(i^\frac{1}{2} h) &= \frac{i h^2}{2} f^{(2)}(x_0) - \frac{h^4}{24} f^{(4)}(x_0) - + \quad f_{even}(i^\frac{1}{2} h) &= \frac{i h^2}{2} f^{(2)}(x_0) - \frac{h^4}{24} f^{(4)}(x_0) - \frac{i h^6}{720} f^{(6)}(x_0) + \\ :label: 12e - & \frac{i h^6}{720} f^{(6)}(x_0) + \frac{h^8}{40320} f^{(8)}(x_0) + ...\\ + & \frac{h^8}{40320} f^{(8)}(x_0) + \frac{i h^{10}}{3628800} f^{(10)}(x_0) -...\\ -Similarly dividing with :math:`h` and taking only the imaginary components yields: +Similarly dividing with :math:`h^2/2` and taking only the imaginary components yields: .. math:: - f^{(2)}(x_0) = 2 \Im(f_{even}(i^\frac{1}{2} h)) / h^2 + \frac{h^4}{360} f^{(6)}(x_0) + ... + \quad f^{(2)}(x_0) = \Im\,(2\,f_{even}(i^\frac{1}{2} h)) / h^2 + \frac{h^4}{360} f^{(6)}(x_0) - \frac{h^8}{1814400} f^{(10)}(x_0)... :label: 12f This approximation is still subject to difference errors, but the error associated with this approximation is proportional to :math:`h^4`. Neglecting these higher order terms yields: .. math:: - f^{(2)}(x_0) = 2 \Im(f_{even}(i^\frac{1}{2} h)) / h^2 = Im(f(x_0 + i^\frac{1}{2} h) + f(x_0-i^\frac{1}{2} h)) / h^2 + \quad f^{(2)}(x_0) = 2 \Im\,(f_{even}(i^\frac{1}{2} h)) / h^2 = \Im(f(x_0 + i^\frac{1}{2} h) + f(x_0-i^\frac{1}{2} h)) / h^2 :label: 12g See [LaiCrassidisCheng2005]_ and [Ridout2009]_ for more details. +The complex-step derivative in numdifftools.Derivative has truncation error +:math:`O(\delta^4)` for both odd and even order derivatives for :math:`n>1`. For :math:`n=1` +the truncation error is on the order of :math:`O(\delta^4)`, so +truncation error can be eliminated by choosing steps to be very small. The first order complex-step derivative avoids the problem of +round-off error with small steps because there is no subtraction. However, +the function to differentiate needs to be analytic. This method does not work if it does +not support complex numbers or involves non-analytic functions such as +e.g.: abs, max, min. For this reason the `central` method is the default method. + High order derivative ##################### -So how do we construct these higher order formulas? In order to compute the 6'th order approximation we simply set :math:`f_{odd}(\delta)` equal to its 3-term Taylor expansion: +So how do we construct these higher order approximation formulas? Here we will deomonstrate the principle by computing the 6'th order central approximation for the first-order derivative. In order to do so we simply set :math:`f_{odd}(\delta)` equal to its 3-term Taylor expansion: .. math:: f_{odd}(\delta) = \sum_{i=0}^{2} \frac{\delta^{2i+1}}{(2i+1)!} f^{(2i+1)}(x_0) @@ -209,9 +212,10 @@ The solution of these equations are simply: f_{odd}(\delta/2) \\ f_{odd}(\delta/4) \end{bmatrix} - :label: 14 + :label: 14a -This solution gives also the 3'rd and 5'th order derivatives as bonus. As previously noted these formulas have the additional benefit of beeing applicable to any scale, with only a scale factor applied. +The first row of :eq:`14a` gives the coefficients for 6'th order approximation. Looking at at row two and three, we see also that +this gives the 6'th order approximation for the 3'rd and 5'th order derivatives as bonus. Thus this is also a general method for obtaining high order differentiation rules. As previously noted these formulas have the additional benefit of beeing applicable to any scale, with only a scale factor applied. Richardson extrapolation methodology applied to derivative estimation @@ -279,7 +283,7 @@ A quick test in Python yields much better results yet. >>> allclose(1./3*df1 - 2*df2 + 8./3*df3, 1.00000539448361) True -Again, Derivative uses the appropriate multiple term Richardson extrapolants for all derivatives of :math:`f` up to the fourth order. This, combined with the use of high order approximations for the derivatives, allows the use of quite large step sizes. See [LynessMoler1966]_ and [LynessMoler1969]_. +Again, Derivative uses the appropriate multiple term Richardson extrapolants for all derivatives of :math:`f` up to any order [3]_. This, combined with the use of high order approximations for the derivatives, allows the use of quite large step sizes. See [LynessMoler1966]_ and [LynessMoler1969]_. How to compute the multiple term Richardson extrapolants will be elaborated further in the next section. Multiple term Richardson extrapolants @@ -360,32 +364,36 @@ These error estimates are also of value in a different sense. Since they are eff >>> import numdifftools as nd >>> from numpy import exp, allclose - >>> f = nd.Derivative(exp) - >>> allclose(f(1), 2.71828183) + >>> f = nd.Derivative(exp, full_output=True) + >>> val, info = f(1) + >>> allclose(val, 2.71828183) True - >>> allclose(f.error_estimate, 5.21804822e-14) + >>> allclose(info.error_estimate, 6.927791673660977e-14) True - >>> allclose(f.final_delta, 0.02166085) + >>> allclose(info.final_step, 0.0078125) True However, if we force the step size to be artificially large, then approximation error takes over. - >>> f = nd.Derivative(exp, delta=1, step_nom=1) - >>> allclose(f(1), 3.19452805) + >>> f = nd.Derivative(exp, step=1, full_output=True) + >>> val, info = f(1) + >>> allclose(val, 3.19452805) + True + >>> allclose(val-exp(1), 0.47624622) True - >>> allclose(f(1)-exp(1), 0.47624622) + >>> allclose(info.final_step, 1) True - >>> f.final_delta array([ 1.]) And if the step size is forced to be too small, then we see noise dominate the problem. - >>> f = nd.Derivative(exp, delta=.0000000001, step_nom=1) - >>> allclose(f(1), 2.71828093) + >>> f = nd.Derivative(exp, step=1e-10, full_output=True) + >>> val, info = f(1) + >>> allclose(val, 2.71828093) True - >>> allclose(f(1) - exp(1), -8.97648138e-07) + >>> allclose(val - exp(1), -8.97648138e-07) True - >>> allclose(f.final_delta, 1.00000008e-10) + >>> allclose(info.final_step, 1.0000000e-10) True @@ -397,30 +405,31 @@ Derivative in action How does numdifftools.Derivative work in action? A simple nonlinear function with a well known derivative is :math:`e^x`. At :math:`x = 0`, the derivative should be 1. - >>> f = nd.Derivative(exp) - >>> f(0) - array([ 1.]) + >>> f = nd.Derivative(exp, full_output=True) + >>> val, info = f(0) + >>> allclose(val, 1) + True - >>> allclose(f.error_estimate, 5.28466160e-14) + >>> allclose(info.error_estimate, 5.28466160e-14) True A second simple example comes from trig functions. The first four derivatives of the sine function, evaluated at :math:`x = 0`, should be respectively :math:`[cos(0), -sin(0), -cos(0), sin(0)]`, or :math:`[1,0,-1,0]`. >>> from numpy import sin, allclose >>> import numdifftools as nd - >>> df = nd.Derivative(sin, 1) + >>> df = nd.Derivative(sin, n=1) >>> allclose(df(0), 1.) True - >>> ddf = nd.Derivative(sin, 2) + >>> ddf = nd.Derivative(sin, n=2) >>> allclose(ddf(0), 0.) True - >>> dddf = nd.Derivative(sin, 3) + >>> dddf = nd.Derivative(sin, n=3) >>> allclose(dddf(0), -1.) True - >>> ddddf = nd.Derivative(sin, 4) + >>> ddddf = nd.Derivative(sin, n=4) >>> allclose(ddddf(0), 0.) True @@ -428,7 +437,7 @@ A second simple example comes from trig functions. The first four derivatives of Gradient and Hessian estimation ################################ -Estimation of the gradient vector (numdifftools.Gradient) of a function of multiple variables is a simple task, requiring merely repeated calls to numdifftools.Derivative. Likewise, the diagonal elements of the hessian matrix are merely pure second partial derivatives of a function. numdifftools.Hessdiag accomplishes this task, again calling numdifftools.Derivative multiple times. Efficient computation of the off-diagonal (mixed partial derivative) elements of the Hessian matrix uses a scheme much like that of numdifftools.Derivative, wherein numdifftools.Derivative is called to determine an initial step size, then Richardson extrapolation is used to improve a set of second order finite difference estimates of those mixed partials. +Estimation of the gradient vector (numdifftools.Gradient) of a function of multiple variables is a simple task, requiring merely repeated calls to numdifftools.Derivative. Likewise, the diagonal elements of the hessian matrix are merely pure second partial derivatives of a function. numdifftools.Hessdiag accomplishes this task, again calling numdifftools.Derivative multiple times. Efficient computation of the off-diagonal (mixed partial derivative) elements of the Hessian matrix uses a scheme much like that of numdifftools.Derivative, then Richardson extrapolation is used to improve a set of second order finite difference estimates of those mixed partials. Conclusion ########## @@ -439,7 +448,7 @@ numdifftools.Derivative is an a adaptive scheme that can compute the derivative Acknowledgments ############### -My thanks are due to Shaun Simmons for convincing me to learn enough LaTeX to write this document. +Thanks are due to John R. D’Errico for writing the first version of this document and `derivest` in which numdifftools is based on. References ########## @@ -473,3 +482,7 @@ References given point. An even symmetry has the property that :math:`f(x) = f(-x)`. Likewise, an odd function expresses an odd symmetry, wherein :math:`f(x) = -f(-x)`. + +.. [3] For practical purposes the maximum order of the derivative is between 4 and 10 + depending on the function to differentiate and also the method used + in the approximation. diff --git a/docs/src/theory.rst b/docs/src/theory.rst index 56f3aa0..59237e3 100644 --- a/docs/src/theory.rst +++ b/docs/src/theory.rst @@ -1,5 +1,5 @@ -Theory -====== +Numerical differentiation +========================= .. toctree:: :maxdepth: 2 diff --git a/numdifftools/core.py b/numdifftools/core.py index 69318bc..6ddafa4 100644 --- a/numdifftools/core.py +++ b/numdifftools/core.py @@ -1,40 +1,57 @@ -""" -Numdifftools implementation +# !/usr/bin/env python +"""numerical differentiation functions: +Derivative, Gradient, Jacobian, and Hessian -""" +Author: Per A. Brodtkorb +Created: 01.08.2008 +Copyright: (c) pab 2008 +Licence: New BSD -# Author: Per A. Brodtkorb -# -# Created: 01.08.2008 -# Copyright: (c) pab 2008 -# Licence: New BSD -# -# Based on matlab functions derivest.m gradest.m hessdiag.m, hessian.m -# and jacobianest.m by: -# -# Author: John D'Errico -# e-mail: woodchips@rochester.rr.com -# Release: 1.0 -# Release date: 12/27/2006 +Based on matlab functions derivest.m gradest.m hessdiag.m, hessian.m + and jacobianest.m version 1.0 released 12/27/2006 by John D'Errico +(e-mail: woodchips@rochester.rr.com) -# !/usr/bin/env python +Also based on the python functions approx_fprime, approx_fprime_cs, + approx_hess_cs, approx_hess1, approx_hess2 and approx_hess3 in the + statsmodels.tools.numdiff module released in 2014 written by Josef Perktold. + +""" from __future__ import division, print_function import numpy as np -import scipy.linalg as linalg -import scipy.misc as misc +from collections import namedtuple +from matplotlib import pyplot as plt +from numdifftools.multicomplex import bicomplex +from numdifftools.test_functions import get_test_function, function_names +from numpy import linalg +from scipy import misc from scipy.ndimage.filters import convolve1d -# import scipy.interpolate as si import warnings -import matplotlib.pyplot as plt __all__ = [ - 'dea3', 'Derivative', 'Jacobian', 'Gradient', 'Hessian', 'Hessdiag' + 'dea3', 'Derivative', 'Jacobian', 'Gradient', 'Hessian', 'Hessdiag', + 'MinStepGenerator', 'MaxStepGenerator' ] - +# NOTE: we only do double precision internally so far _TINY = np.finfo(float).tiny _EPS = np.finfo(float).eps +EPS = np.MachAr().eps +_SQRT_J = (1j + 1.0) / np.sqrt(2.0) # = 1j**0.5 + +_CENTRAL_WEIGHTS_AND_POINTS = { + (1, 3): (np.array([-1, 0, 1]) / 2.0, np.arange(-1, 2)), + (1, 5): (np.array([1, -8, 0, 8, -1]) / 12.0, np.arange(-2, 3)), + (1, 7): (np.array([-1, 9, -45, 0, 45, -9, 1]) / 60.0, np.arange(-3, 4)), + (1, 9): (np.array([3, -32, 168, -672, 0, 672, -168, 32, -3]) / 840.0, + np.arange(-4, 5)), + (2, 3): (np.array([1, -2.0, 1]), np.arange(-1, 2)), + (2, 5): (np.array([-1, 16, -30, 16, -1]) / 12.0, np.arange(-2, 3)), + (2, 7): (np.array([2, -27, 270, -490, 270, -27, 2]) / 180.0, + np.arange(-3, 4)), + (2, 9): (np.array([-9, 128, -1008, 8064, -14350, + 8064, -1008, 128, -9]) / 5040.0, + np.arange(-4, 5))} class Dea(object): @@ -299,400 +316,808 @@ def dea3(v0, v1, v2, symmetric=False): return result, abserr -def vec2mat(vec, n, m): - ''' forms the matrix M, such that M[i,j] = vec[i+j] +def fornberg_weights_all(x, x0, M=1): + ''' + Return finite difference weights_and_points for derivatives + of all orders 0, 1, ..., m + + Parameters + ---------- + x : vector, length n + x-coordinates for grid points + x0 : scalar + location where approximations are to be accurate + m : scalar integer + highest derivative that we want to find weights_and_points for + + Returns + ------- + C : array, shape n x m+1 + contains coefficients for the j'th derivative in column j (0 <= j <= m) + + See also: + --------- + fornberg_weights + + References + ---------- + B. Fornberg (1998) + "Calculation of weights_and_points in finite difference formulas", + SIAM Review 40, pp. 685-691. + + http://www.scholarpedia.org/article/Finite_difference_method + ''' + N = len(x) + if M >= N: + raise ValueError('length(x) must be larger than m') + + c1, c4 = 1, x[0] - x0 + C = np.zeros((N, M + 1)) + C[0, 0] = 1 + for n in range(1, N): + m = np.arange(0, min(n, M) + 1) + c2, c5, c4 = 1, c4, x[n] - x0 + for v in range(n): + c3 = x[n] - x[v] + c2, c6, c7 = c2 * c3, m * C[v, m-1], C[v, m] + C[v, m] = (c4 * c7 - c6) / c3 + else: + C[n, m] = c1 * (c6 - c5 * c7) / c2 + c1 = c2 + return C + + +def fornberg_weights(x, x0, m=1): ''' - [i, j] = np.ogrid[0:n, 0:m] - return np.matrix(vec[i + j]) + Return weights for finite difference approximation of the m'th derivative + U^m(x0), evaluated at x0, based on n values of U at x[0], x[1],... x[n-1]: + + U^m(x0) = sum weights[i] * U(x[i]) + + Parameters + ---------- + x : vector + abscissas used for the evaluation for the derivative at x0. + x0 : scalar + location where approximations are to be accurate + m : integer + order of derivative. Note for m=0 this can be used to evaluate the + interpolating polynomial itself. + + Notes + ----- + The x values can be arbitrarily spaced but must be distinct and len(x) > m. + + The Fornberg algorithm is much more stable numerically than regular + vandermonde systems for large values of n. + + See also + -------- + fornberg_weights_all + ''' + return fornberg_weights_all(x, x0, m)[:, -1] def _make_exact(h): - '''Make sure h is an exact representable number - This is important when calculating numerical derivatives and is - accomplished by adding 1 and then subtracting 1.. - ''' - return (h + 1.0) - 1.0 + '''Make sure h is an exact representable number + This is important when calculating numerical derivatives and is + accomplished by adding 1 and then subtracting 1.. + ''' + return (h + 1.0) - 1.0 + + +def default_scale(method='forward', n=1, order=2): + # is_odd = (n % 2) == 1 + high_order = int(n > 1 or order >= 4) + order2 = max(order // 2-1, 0) + n4 = n // 4 + return (dict(multicomplex=1.35, complex=1.35).get(method, 2.5) + + int((n - 1)) * dict(multicomplex=0, complex=0.0).get(method, 1.3) + + order2 * dict(central=3, forward=2, backward=2).get(method, 0) + + # is_odd * dict(complex=2.65*int(n//2)).get(method, 0) + + (n % 4 == 1) * high_order * dict(complex=3.65 + n4 * (5 + 1.5**n4) + ).get(method, 0) + + (n % 4 == 3) * dict(complex=3.65*2 + n4 * (5 + 2.1**n4) + ).get(method, 0) + + (n % 4 == 2) * dict(complex=3.65 + n4 * (5 + 1.7**n4) + ).get(method, 0) + + (n % 4 == 0) * dict(complex=(n//4) * (10 + 1.5*int(n > 10)) + ).get(method, 0)) + + +def valarray(shape, value=np.NaN, typecode=None): + """Return an array of all value. + """ + if typecode is None: + typecode = bool + out = np.ones(shape, dtype=typecode) * value + + if not isinstance(out, np.ndarray): + out = np.asarray(out) + return out + + +def nom_step(x=None): + '''Return nominal step''' + if x is None: + return 1.0 + return np.maximum(np.log1p(np.abs(x)), 1.0) -def _get_epsilon(x, s, epsilon, n): +def _default_base_step(x, scale, epsilon=None): if epsilon is None: - h = (10*_EPS) ** (1. / s) * np.maximum(np.log1p(np.abs(x)), 0.1) + h = EPS ** (1. / scale) * nom_step(x) else: - if np.isscalar(epsilon): - h = np.empty(n) - h.fill(epsilon) - else: # pragma : no cover - h = np.asarray(epsilon) - if h.shape != x.shape: - raise ValueError("If h is not a scalar it must have the same" - " shape as x.") - return _make_exact(h) - - -class RombergExtrapolation(object): + h = valarray(x.shape, value=epsilon) + return h + + +class MinStepGenerator(object): ''' + Generates a sequence of steps + + where + steps = base_step * step_ratio ** (np.arange(num_steps) + offset) + Parameters + ---------- + base_step : float, array-like, optional + Defines the base step, if None, then base_step is set to + EPS**(1/scale)*max(log(1+|x|), 1) + where x is supplied at runtime through the __call__ method. + step_ratio : real scalar, optional, default 2 + Ratio between sequential steps generated. + Note: Ratio > 1 + If None then step_ratio is 2 for n=1 otherwise step_ratio is 1.6 + num_steps : scalar integer, optional, default n + order - 1 + num_extrap + defines number of steps generated. It should be larger than + n + order - 1 + offset : real scalar, optional, default 0 + offset to the base step + scale : real scalar, optional + scale used in base step. If not None it will override the default + computed with the default_scale function. ''' - def __init__(self, step_ratio=2.0, method='central', order=2, num_terms=2): - self.num_terms = num_terms - self.order = order - self.method = method + + def __init__(self, base_step=None, step_ratio=2, num_steps=None, + offset=0, scale=None, num_extrap=0, use_exact_steps=True, + check_num_steps=True): + self.base_step = base_step + self.num_steps = num_steps self.step_ratio = step_ratio + self.offset = offset + self.scale = scale + self.check_num_steps = check_num_steps + self.use_exact_steps = use_exact_steps + self.num_extrap = num_extrap + + def __repr__(self): + class_name = self.__class__.__name__ + kwds = ['%s=%s' % (name, str(getattr(self, name))) + for name in self.__dict__.keys()] + return """%s(%s)""" % (class_name, ','.join(kwds)) + + def _default_scale(self, method, n, order): + scale = self.scale + if scale is None: + scale = default_scale(method, n, order) + return scale + + def _default_base_step(self, xi, method, n, order=2): + scale = self._default_scale(method, n, order) + base_step = _default_base_step(xi, scale, self.base_step) + if self.use_exact_steps: + base_step = _make_exact(base_step) + return base_step + + def _min_num_steps(self, method, n, order): + num_steps = n + order - 1 + + if method in ['central', 'central2', 'complex', 'multicomplex']: + step = 2 + if method == 'complex': + step = 4 if n > 2 or order >= 4 else 2 + num_steps = (n + order-1) // step + return max(int(num_steps), 1) + + def _default_num_steps(self, method, n, order): + min_num_steps = self._min_num_steps(method, n, order) + if self.num_steps is not None: + num_steps = int(self.num_steps) + if self.check_num_steps: + num_steps = max(num_steps, min_num_steps) + return num_steps + return min_num_steps + int(self.num_extrap) + + def _default_step_ratio(self, n): + if self.step_ratio is None: + step_ratio = {1: 2.0}.get(n, 1.6) + else: + step_ratio = float(self.step_ratio) + if self.use_exact_steps: + step_ratio = _make_exact(step_ratio) + return step_ratio - def _get_romberg_rule(self): - num_terms = self.num_terms - rmat = np.ones((num_terms + 2, num_terms + 1)) - if num_terms > 0: - fact = dict(central=2).get(self.method, 1) - i, j = np.ogrid[0:num_terms + 2, 0:num_terms] - srinv = _make_exact(1.0 / self.step_ratio) - rmat[:, 1:] = srinv ** (i*(fact*j + self.order)) - romberg_rule = linalg.pinv(rmat)[0] - return romberg_rule - - def _estimate_error(self, der_romb): - m, n = der_romb.shape - z = np.zeros((n, )) - tmperr = np.abs(np.diff(np.vstack((z, der_romb, z)), axis=0)) * (m > 1) - abserr = tmperr[:-1] + tmperr[1:] + np.abs(der_romb) * _EPS * 10.0 - return abserr + def __call__(self, x, method='central', n=1, order=2): + xi = np.asarray(x) + base_step = self._default_base_step(xi, method, n, order) + step_ratio = self._default_step_ratio(n) - def __call__(self, der, steps): - ne = der.shape[0] - if ne < self.num_terms + 2: - der_romb = der - nr = 0 - else: - rr_rule = self._get_romberg_rule() - nr = rr_rule.size - 1 - m = ne - nr - der_romb = convolve1d(der, rr_rule[::-1], axis=0, - origin=(nr//2))[:m] - abserr = self._estimate_error(der_romb) - return der_romb, abserr, steps[nr:] - - -def romberg_rule(step_ratio, order, method, num_terms=2): - rmat = np.ones((num_terms + 2, num_terms + 1)) - if num_terms > 0: - fact = dict(central=2).get(method, 1) - i, j = np.ogrid[0:num_terms + 2, 0:num_terms] - srinv = _make_exact(1.0 / step_ratio) - rmat[:, 1:] = srinv ** (i*(fact*j + order)) - romberg_rule = linalg.pinv(rmat)[0] - return romberg_rule + num_steps = self._default_num_steps(method, n, order) + offset = self.offset + for i in range(num_steps-1, -1, -1): + h = (base_step * step_ratio**(i + offset)) + if (np.abs(h) > 0).all(): + yield h -class _Derivative(object): +class MinMaxStepGenerator(object): + ''' + Generates a sequence of steps - ''' Object holding common variables and methods for the numdifftools + where + steps = logspace(log10(step_min), log10(step_max), num_steps) Parameters ---------- - fun : callable - function to differentiate. - n : Integer from 1 to 4 (Default 1) - defining derivative order. - order : Integer from 1 to 4 (Default 2) - defining order of basic method used. - For 'central' methods, it must be from the set [2,4]. - method : Method of estimation. Valid options are: - 'central', 'forward' or 'backward'. (Default 'central') - romberg_terms : integer from 0 to 3 (Default 2) - Number of Romberg terms used in the extrapolation. - Note: 0 disables the Romberg step completely. - step_nom : vector default maximum(log1p(abs(x0)), 0.1) + step_min : float, array-like, optional + Defines the minimim step. Default value is: + EPS**(1/scale)*max(log(1+|x|), 1) + where x and scale are supplied at runtime through the __call__ method. + step_max : real scalar, optional + maximum step generated. Default value is: + exp(log(step_min) * scale / (scale + 1.5)) + num_steps : scalar integer, optional + defines number of steps generated. + scale : real scalar, optional + scale used in base step. If set to a value it will override the scale + supplied at runtime. + ''' + + def __init__(self, step_min=None, step_max=None, num_steps=10, scale=None, + num_extrap=0): + self.step_min = step_min + self.num_steps = num_steps + self.step_max = step_max + self.scale = scale + self.num_extrap = num_extrap + + def __repr__(self): + class_name = self.__class__.__name__ + kwds = ['%s=%s' % (name, str(getattr(self, name))) + for name in self.__dict__.keys()] + return """%s(%s)""" % (class_name, ','.join(kwds)) + + def __call__(self, x, method='forward', n=1, order=None): + if self.scale is not None: + scale = self.scale + xi = np.asarray(x) + step_min, step_max = self.step_min, self.step_max + delta = _default_base_step(xi, scale, step_min) + if step_min is None: + step_min = (10 * EPS)**(1. / scale) + if step_max is None: + step_max = np.exp(np.log(step_min) * scale / (scale + 1.5)) + steps = np.logspace(0, np.log10(step_max) - np.log10(step_min), + self.num_steps)[::-1] + + for step in steps: + h = _make_exact(delta * step) + if (np.abs(h) > 0).all(): + yield h + +''' + step_nom : vector default maximum(log1p(abs(x0)), 1) Nominal step. (The steps: h_i = step_nom[i] * delta) step_max : real scalar (Default 2.0) Maximum allowed excursion from step_nom as a multiple of it. step_ratio: real scalar (Default 2.0) Ratio used between sequential steps in the estimation of the derivative step_num : integer (Default 26) - The minimum step_num for making romberg extrapolation work is - 7 + np.ceil(self.n/2.) + self.order + self.romberg_terms + The minimum step_num for making richardson extrapolation work is + 7 + np.ceil(self.n/2.) + self.order + self.richardson_terms delta : vector default step_max*step_ratio**(-arange(step_num)) Defines the steps sizes used in derivation: h_i = step_nom[i] * delta - vectorized : Bool - True - if your function is vectorized. - False - loop over the successive function calls (default). - - Uses a semi-adaptive scheme to provide the best estimate of the - derivative by its automatic choice of a differencing interval. It uses - finite difference approximations of various orders, coupled with a - generalized (multiple term) Romberg extrapolation. This also yields the - error estimate provided. See the document DERIVEST.pdf for more explanation - of the algorithms behind the parameters. - - Note on order: higher order methods will generally be more accurate, - but may also suffer more from numerical problems. First order - methods would usually not be recommended. - Note on method: Central difference methods are usually the most accurate, - but sometimes one can only allow evaluation in forward or backward - direction. +''' + + +class MaxStepGenerator(MinStepGenerator): ''' + Generates a sequence of steps - def __init__(self, fun, n=1, order=2, method='central', romberg_terms=2, - step_max=2.0, step_nom=None, step_ratio=2.0, step_num=26, - offset=-2, - delta=None, vectorized=False, verbose=False, - use_dea=True, transform=None): - self.fun = fun - self.n = n - self.order = order - self.method = method - self.romberg_terms = romberg_terms + where + steps = base_step * step_ratio ** (-np.arange(num_steps) + offset) + base_step = step_max * step_nom + + Parameters + ---------- + max_step : float, array-like, optional default 2 + Defines the maximum step + step_ratio : real scalar, optional, default 2 + Ratio between sequential steps generated. + Note: Ratio > 1 + num_steps : scalar integer, optional, default n + order - 1 + num_extrap + defines number of steps generated. It should be larger than + n + order - 1 + step_nom : default maximum(log1p(abs(x)), 1) + Nominal step. + offset : real scalar, optional, default 0 + offset to the base step: max_step * nom_step + + ''' + def __init__(self, step_max=2.0, step_ratio=2.0, num_steps=15, + step_nom=None, offset=0, num_extrap=0, + use_exact_steps=False, check_num_steps=True): + self.base_step = None self.step_max = step_max self.step_ratio = step_ratio - self.offset = offset + self.num_steps = num_steps self.step_nom = step_nom - self.step_num = step_num - self.delta = delta - self.vectorized = vectorized - self.verbose = verbose - self.use_dea = use_dea - self.transform = transform - - self._check_params() - - self.error_estimate = None - self.final_delta = None - - # The remaining member variables are set by _initialize - self._fd_rule = None - # self._rmat = None - self._qromb = None - self._rromb = None - self._diff_fun = None - - def _set_delta(self, delta=None): - ''' Set the steps to use in derivation. - - Member variables used: - step_ratio, step_num, step_max - ''' - if delta is None: - step_ratio = self._make_exact(self.step_ratio) - if self.step_num is None: - num_steps = self._get_min_num_steps() - else: - num_steps = max(self.step_num, 1) - step1 = self._make_exact(self.step_max) - self._delta = step1 * step_ratio ** (-np.arange(num_steps)) - else: - self._delta = np.atleast_1d(delta) + self.offset = offset + self.num_extrap = num_extrap + self.check_num_steps = check_num_steps + self.use_exact_steps = use_exact_steps - delta = property(lambda cls: cls._delta, fset=_set_delta) - finaldelta = property(lambda cls: cls.final_delta) + def _default_step_nom(self, x): + if self.step_nom is None: + return nom_step(x) + return valarray(x.shape, value=self.step_nom) - def _check_params(self): - ''' check the parameters for acceptability - ''' - atleast_1d = np.atleast_1d - kwds = self.__dict__ - for name in ['n', 'order']: - val = np.atleast_1d(kwds[name]) - if ((len(val) != 1) or (val != int(val)) or val < 0): - raise ValueError('%s must be positive scalar integer.' % name) - name = 'romberg_terms' - val = atleast_1d(kwds[name]) - if not ((len(val) == 1) and (val in (0, 1, 2, 3))): - raise ValueError('%s must be scalar, one of [0 1 2 3].' % name) - - for name in ('step_max', 'step_num'): - val = kwds[name] - if val is not None and (len(atleast_1d(val)) > 1 or val <= 0): - raise ValueError('%s must be None or a scalar, >0.' % name) - - valid_methods = dict(c='central', f='forward', b='backward') - method = valid_methods.get(kwds['method'][0]) - if method is None: - t = 'Invalid method: Must start with one of c, f, b characters!' - raise ValueError(t) - order = kwds['order'] - if method[0] == 'c' and (order % 2) == 1: - raise ValueError('order %d is not possible for central ' - 'difference methods' % order) - - def _initialize(self): - '''Set derivative parameters: - differention rule and romberg extrapolation matrices - ''' - self._set_fd_rule() - self._set_romb_qr() - self._set_difference_function() + def _default_base_step(self, xi, method, n): + base_step = self.base_step + if base_step is None: + base_step = self.step_max * self._default_step_nom(xi) + if self.use_exact_steps: + base_step = _make_exact(base_step) + return base_step - def _fder(self, fun, f_x0i, x0i, h): - ''' - Return derivative estimates of f at x0 for a sequence of stepsizes h + def __call__(self, x, method='forward', n=1, order=None): + xi = np.asarray(x) - Member variables used - --------------------- - n - _fd_rule - romberg_terms - ''' - fd_rule = self._fd_rule - n_fdr = fd_rule.size - n_h = h.size + offset = self.offset + + base_step = self._default_base_step(xi, method, n) + step_ratio = self._default_step_ratio(n) + + num_steps = self._default_num_steps(method, n, order) + for i in range(num_steps): + h = base_step * step_ratio**(-i + offset) + if (np.abs(h) > 0).all(): + yield h + + +class Richardson(object): + ''' + Extrapolates as sequence with Richardsons method + + Notes + ----- + Suppose you have series expansion that goes like this + + L = f(h) + a0 * h^p_0 + a1 * h^p_1+ a2 * h^p_2 + ... + + where p_i = order + step * i and f(h) -> L as h -> 0, but f(0) != L. + + If we evaluate the right hand side for different stepsizes h + we can fit a polynomial to that sequence of approximations. + This is exactly what this class does. + + Example + ------- + >>> import numpy as np + >>> import numdifftools.nd_cstep as nd + >>> n = 3 + >>> Ei = np.zeros((n,1)) + >>> h = np.zeros((n,1)) + >>> linfun = lambda i : np.linspace(0, np.pi/2., 2**(i+5)+1) + >>> for k in np.arange(n): + ... x = linfun(k) + ... h[k] = x[1] + ... Ei[k] = np.trapz(np.sin(x),x) + >>> En, err, step = nd.Richardson(step=1, order=1)(Ei, h) + >>> truErr = Ei-1. + >>> (truErr, err, En) + (array([[ -2.00805680e-04], + [ -5.01999079e-05], + [ -1.25498825e-05]]), array([[ 0.00111851]]), array([[ 1.]])) + + ''' + def __init__(self, step_ratio=2.0, step=1, order=1, num_terms=2): + self.num_terms = num_terms + self.order = order + self.step = step + self.step_ratio = step_ratio + + def _r_matrix(self, num_terms): + step = self.step + i, j = np.ogrid[0:num_terms+1, 0:num_terms] + r_mat = np.ones((num_terms + 1, num_terms + 1)) + r_mat[:, 1:] = (1.0 / self.step_ratio) ** (i*(step*j + self.order)) + return r_mat + + def _get_richardson_rule(self, sequence_length=None): + if sequence_length is None: + sequence_length = self.num_terms + 1 + num_terms = min(self.num_terms, sequence_length - 1) + if num_terms > 0: + r_mat = self._r_matrix(num_terms) + return linalg.pinv(r_mat)[0] + return np.ones((1,)) + + def _estimate_error(self, new_sequence, old_sequence, steps, rule): + m, _n = new_sequence.shape + + if m < 2: + return (np.abs(new_sequence) * EPS + steps) * 10.0 + cov1 = np.sum(rule**2) # 1 spare dof + fact = np.maximum(12.7062047361747 * np.sqrt(cov1), EPS * 10.) + err = np.abs(np.diff(new_sequence, axis=0)) * fact + tol = np.maximum(np.abs(new_sequence[1:]), + np.abs(new_sequence[:-1])) * EPS * fact + converged = err <= tol + abserr = err + np.where(converged, tol * 10, + abs(new_sequence[:-1]-old_sequence[1:])*fact) + # abserr = err1 + err2 + np.where(converged, tol2 * 10, abs(result-E2)) + # abserr = s * fact + np.abs(new_sequence) * EPS * 10.0 + return abserr + + def extrapolate(self, sequence, steps): + return self.__call__(sequence, steps) + + def __call__(self, sequence, steps): + ne = sequence.shape[0] + rule = self._get_richardson_rule(ne) + nr = rule.size - 1 + m = ne - nr + new_sequence = convolve1d(sequence, rule[::-1], axis=0, origin=(nr//2)) + abserr = self._estimate_error(new_sequence, sequence, steps, rule) + return new_sequence[:m], abserr[:m], steps[:m] - f_del = self._diff_fun(fun, f_x0i, x0i, h) +_cmn_doc = """ + Calculate %(derivative)s with finite difference approximation + + Parameters + ---------- + f : function + function of one array f(x, `*args`, `**kwargs`) + step : float, array-like or StepGenerator object, optional + Spacing used, if None, then the spacing is automatically chosen + according to (10*EPS)**(1/scale)*max(log(1+|x|), 1) where scale is + depending on method and derivative-order (see default_scale). + A StepGenerator can be used to extrapolate the results. However, + the generator must generate minimum 3 steps in order to extrapolate + the values. + method : string, optional + defines method used in the approximation + 'central': central difference derivative + 'complex': complex-step derivative + 'multicomplex': multicomplex derivative + 'backward': backward difference derivative + 'forward': forward difference derivative%(extra_parameter)s + full_output : bool, optional + If `full_output` is False, only the derivative is returned. + If `full_output` is True, then (der, r) is returned `der` is the + derivative, and `r` is a Results object. + + Call Parameters + --------------- + x : array_like + value at which function derivative is evaluated + args : tuple + Arguments for function `f`. + kwds : dict + Keyword arguments for function `f`. + %(returns)s + Notes + ----- + The complex-step derivative has truncation error O(steps**2) and + O(steps**4) for odd and even order derivatives respectively, so + truncation error can be eliminated by choosing steps to be very small. + Especially the first order complex-step derivative avoids the problem of + round-off error with small steps because there is no subtraction. However, + the function needs to be analytic. This method does not work if f(x) does + not support complex numbers or involves non-analytic functions such as + e.g.: abs, max, min. + For this reason the 'central' method is the default method. + This method is usually very accurate, but sometimes one can only allow + evaluation in forward or backward direction. + + Be careful in decreasing the step size too much due to round-off errors. + %(extra_note)s + References + ---------- + Ridout, M.S. (2009) Statistical applications of the complex-step method + of numerical differentiation. The American Statistician, 63, 66-74 + + K.-L. Lai, J.L. Crassidis, Y. Cheng, J. Kim (2005), New complex step + derivative approximations with application to second-order + kalman filtering, AIAA Guidance, Navigation and Control Conference, + San Francisco, California, August 2005, AIAA-2005-5944. + + Lyness, J. M., Moler, C. B. (1966). Vandermonde Systems and Numerical + Differentiation. *Numerische Mathematik*. + + Lyness, J. M., Moler, C. B. (1969). Generalized Romberg Methods for + Integrals of Derivatives. *Numerische Mathematik*. + %(example)s + %(see_also)s + """ - if f_del.size != n_h: + +class _Derivative(object): + + info = namedtuple('info', ['error_estimate', 'final_step', 'index']) + + def __init__(self, f, step=None, method='central', order=2, n=1, + full_output=False): + self.fun = f + self.n = n + self.order = order + self.method = method + self.full_output = full_output + self.richardson_terms = 2 + self.step = self._make_generator(step) + + def _make_generator(self, step): + if hasattr(step, '__call__'): + return step + if step is None and self.method not in ['complex', 'multicomplex']: + return MaxStepGenerator(step_ratio=None, num_extrap=14) + return MinStepGenerator(base_step=step, step_ratio=None, num_extrap=0) + + def _get_arg_min(self, errors): + shape = errors.shape + try: + arg_mins = np.nanargmin(errors, axis=0) + min_errors = np.nanmin(errors, axis=0) + except ValueError as msg: + warnings.warn(str(msg)) + ix = np.arange(shape[1]) + return ix + + for i, min_error in enumerate(min_errors): + idx = np.flatnonzero(errors[:, i] == min_error) + arg_mins[i] = idx[idx.size // 2] + ix = np.ravel_multi_index((arg_mins, np.arange(shape[1])), shape) + return ix + + def _get_best_estimate(self, der, errors, steps, shape): + ix = self._get_arg_min(errors) + final_step = steps.flat[ix].reshape(shape) + err = errors.flat[ix].reshape(shape) + return der.flat[ix].reshape(shape), self.info(err, final_step, ix) + + @property + def _method_order(self): + step = self._richardson_step() + # Make sure it is even and at least 2 or 4 + order = max((self.order // step) * step, step) + return order + + def _complex_high_order(self): + return self.method == 'complex' and (self.n > 1 or self.order >= 4) + + def _richardson_step(self): + # complex_step = 4 if self.n % 2 == 0 else 2 + complex_step = 4 if self._complex_high_order() else 2 + + return dict(central=2, central2=2, complex=complex_step, + multicomplex=2).get(self.method, 1) + + def _set_richardson_rule(self, step_ratio, num_terms=2): + order = self._method_order + step = self._richardson_step() + self._richardson_extrapolate = Richardson(step_ratio=step_ratio, + step=step, order=order, + num_terms=num_terms) + + def _wynn_extrapolate(self, der, steps): + der, errors = dea3(der[0:-2], der[1:-1], der[2:], symmetric=False) + return der, errors, steps[2:] + + def _extrapolate(self, results, steps, shape): + der, errors, steps = self._richardson_extrapolate(results, steps) + if len(der) > 2: + # der, errors, steps = self._richardson_extrapolate(results, steps) + der, errors, steps = self._wynn_extrapolate(der, steps) + der, info = self._get_best_estimate(der, errors, steps, shape) + return der, info + + def _get_function_name(self): + name = '_%s' % self.method + even_derivative_order = self._is_even_derivative() + if even_derivative_order and self.method in ('central', 'complex'): + name = name + '_even' + if self.method in ('complex') and self._is_fourth_derivative(): + name = name + '_higher' + else: + if self._complex_high_order() and self._is_odd_derivative(): + name = name + '_odd' + if self._is_third_derivative(): + name = name + '_higher' + elif self.method == 'multicomplex' and self.n > 1: + if self.n == 2: + name = name + '2' + else: + raise ValueError('Multicomplex method only support first ' + 'and second order derivatives.') + return name + + def _get_functions(self): + name = self._get_function_name() + return getattr(self, name), self.fun + + def _get_steps(self, xi): + method, n, order = self.method, self.n, self._method_order + return [step for step in self.step(xi, method, n, order)] + + def _is_odd_derivative(self): + return self.n % 2 == 1 + + def _is_even_derivative(self): + return self.n % 2 == 0 + + def _is_third_derivative(self): + return self.n % 4 == 3 + + def _is_fourth_derivative(self): + return self.n % 4 == 0 + + def _eval_first_condition(self): + even_derivative = self._is_even_derivative() + return ((even_derivative and self.method in ('central', 'central2')) or + self.method in ['forward', 'backward'] or + self.method == 'complex' and self._is_fourth_derivative()) + + def _eval_first(self, f, x, *args, **kwds): + if self._eval_first_condition(): + return f(x, *args, **kwds) + return 0.0 + + def _vstack(self, sequence, steps): + # sequence = np.atleast_2d(sequence) + original_shape = np.shape(sequence[0]) + f_del = np.vstack(list(np.ravel(r)) for r in sequence) + h = np.vstack(list(np.ravel(np.ones(original_shape)*step)) + for step in steps) + if f_del.size != h.size: raise ValueError('fun did not return data of correct size ' + '(it must be vectorized)') + return f_del, h, original_shape - # ne = max(n_h + 1 - n_fdr - self.romberg_terms, 1) - ne = max(n_h + 1 - n_fdr, 1) - der_init = np.asarray(vec2mat(f_del, ne, n_fdr) * fd_rule.T).ravel() - der_init = der_init / (h[:ne]) ** self.n - - return der_init, h[:ne] - -# def _trim_estimates(self, der_romb, errors, h): -# ''' -# trim off the estimates at each end of the scale -# ''' -# trimdelta = h.copy() -# der_romb = np.atleast_1d(der_romb) -# num_vals = len(der_romb) -# nr_rem_min = int((num_vals - 1) / 2) -# nr_rem = min(2 * max((self.n - 1), 1), nr_rem_min) -# if nr_rem > 0: -# tags = der_romb.argsort() -# tags = tags[nr_rem:-nr_rem] -# der_romb = der_romb[tags] -# errors = errors[tags] -# trimdelta = trimdelta[tags] -# return der_romb, errors, trimdelta - - def _plot_errors(self, h2, errors, step_nom_i, der_romb): - i = np.argsort(h2) - ii = np.arange(len(h2)) - ii[i] = np.arange(len(h2)) - print(' Stepsize Value Errors') - fmt = ''.join(('%10.2g %20.15g %10.2g\n',)*len(h2)) - tmp = (np.vstack((h2[i], der_romb[i], errors[i])).T).ravel().tolist() - print(fmt % tuple(tmp)) - # ind = self._get_arg_min(errors) - # plt.ioff() - try: - # diff_df = np.diff(np.hstack(((der_romb[i], 0))))[ii] - diff_df = np.diff(np.hstack(((0, der_romb[i], 0))), n=1) # [ii] - diff_df = np.diff(np.log(abs(diff_df)))[ii] - ip = diff_df >= 0 - plt.subplot(2, 1, 1) - plt.semilogx(h2[i], np.abs(diff_df[i]) + _EPS, 'k', - h2[ip], diff_df[ip] + _EPS, 'r.', - h2[~ip], -diff_df[~ip] + _EPS, 'b.') - small = 2 * np.sqrt(_EPS) ** (1. / np.sqrt(self.n)) - plt.vlines(small, 1e-15, 1) - plt.title('abs(df_i(x)-df_best(x)) as function of stepsize. ' + - '(nom=%g)' % step_nom_i) - plt.subplot(2, 1, 2) - plt.loglog(h2[i], errors[i], 'r--', h2, errors, 'r.') - plt.vlines(small, 1e-15, 1) - plt.title('Estimated error as function of stepsize. (nom=%g)' % - step_nom_i) - - # plt.show() - except: - pass - - def _get_arg_min(self, errest): - arg_mins = np.flatnonzero(errest == np.min(errest)) - n = arg_mins.size - return arg_mins[n // 2] - - def _get_step_nom(self, step_nom, x0): - # s = np.sqrt(4*self.n) - # return _get_epsilon(x0, s, step_nom, x0.shape) - if step_nom is None: - step_nom = np.maximum(np.log1p(np.abs(x0)), 0.1) - return self._make_exact(np.atleast_1d(step_nom)) * np.ones(x0.shape) - - def _eval_first(self, fun, x0): - f_x0 = np.zeros(x0.shape) - # will we need fun(x0)? - even_order = np.remainder(self.n, 2) == 0 - if even_order or not self.method[0] == 'c': - if self.vectorized: - f_x0 = fun(x0) - else: - f_x0 = np.asfarray([fun(x0j) for x0j in x0]) - return f_x0 - - def _remove_non_positive(self, h): - threshold = (self.n > 1) * 10.0 ** (-15 + self.n) - if (h <= threshold).any(): - warnings.warn('Some of the steps are too small, either because ' + - 'step_max*step_nom is too small or ' + - 'step_ratio**step_num too large!') - return h[h > threshold] - return h - - def _make_exact(self, h): - '''Make sure h is an exact representable number - This is important when calculating numerical derivatives and is - accomplished by adding 1 and then subtracting 1.. - ''' - return (h + 1.0) - 1.0 + def _compute_step_ratio(self, steps): + if len(steps) < 2: + return 1 + return np.unique(steps[0]/steps[1]).mean() - def _get_steps(self, step_nom): - h = self._make_exact((1.0 * step_nom) * self._delta) - return self._remove_non_positive(h) + def __call__(self, x, *args, **kwds): + xi = np.asarray(x) + results = self._derivative(xi, args, kwds) + derivative, info = self._extrapolate(*results) + if self.full_output: + return derivative, info + return derivative - def _best_der(self, der_romb, errors, h2): - # der_romb, errors, h2 = self._trim_estimates(der_romb, errors, h2) - i = self._get_arg_min(errors) - return der_romb[i], errors[i], h2[i] - def _derivative(self, fun, x00, step_nom=None): - x0 = np.atleast_1d(x00) - step_nom = self._get_step_nom(step_nom, x0) +class Derivative(_Derivative): + __doc__ = _cmn_doc % dict( + derivative='n-th derivative', + extra_parameter=""" + order : int, optional + defines the order of the error term in the Taylor approximation used. + For 'central' and 'complex' methods, it must be an even number. + n : int, optional + Order of the derivative.""", + extra_note=""" + Higher order approximation methods will generally be more accurate, but may + also suffer more from numerical problems. First order methods is usually + not recommended. + """, returns=""" + Returns + ------- + der : ndarray + array of derivatives + """, example=""" + Examples + -------- + >>> import numpy as np + >>> import numdifftools.nd_cstep as ndc - f_x0 = self._eval_first(fun, x0) - n, nx0 = x0.size, x0.shape - der, err, delta = np.zeros(nx0), np.zeros(nx0), np.zeros(nx0) - for i in range(n): - x0i, f_x0i = float(x0[i]), float(f_x0[i]) - h = self._get_steps(step_nom[i]) - der_init, h1 = self._fder(fun, f_x0i, x0i, h) - der_romb, errors, h2 = self._romb_extrap(der_init, h1) - if self.verbose: - self._plot_errors(h2, errors, step_nom[i], der_romb) - der[i], err[i], delta[i] = self._best_der(der_romb, errors, h2) - return der, err, delta - - def _fd_mat(self, parity, nterms): - ''' Return matrix for finite difference derivation. + # 1'st derivative of exp(x), at x == 1 + + >>> fd = ndc.Derivative(np.exp) + >>> np.allclose(fd(1), 2.71828183) + True + + >>> d2 = fd([1, 2]) + >>> np.allclose(d2, [ 2.71828183, 7.3890561 ]) + True + + >>> def f(x): + ... return x**3 + x**2 + + >>> df = ndc.Derivative(f) + >>> np.allclose(df(1), 5) + True + >>> ddf = ndc.Derivative(f, n=2) + >>> np.allclose(ddf(1), 8) + True + """, see_also=""" + See also + -------- + Gradient, + Hessian + """) + """ + Find the n-th derivative of a function at a point. + + Given a function, use a difference formula with spacing `dx` to + compute the `n`-th derivative at `x0`. + + Parameters + ---------- + f : function + Input function. + x0 : float + The point at which `n`-th derivative is found. + dx : float, optional + Spacing. + method : Method of estimation. Valid options are: + 'central', 'forward' or 'backward'. (Default 'central') + n : int, optional (Default 1) + Order of the derivative. + order : int, optional (Default 2) + defining order of basic method used. + For 'central' methods, it must be an even number eg. [2,4]. + + Notes + ----- + Decreasing the step size too small can result in round-off error. + + Note on order: higher order methods will generally be more accurate, + but may also suffer more from numerical problems. First order + methods would usually not be recommended. + Note on method: Central difference methods are usually the most accurate, + but sometimes one can only allow evaluation in forward or backward + direction. + + + """ + @staticmethod + def _fd_matrix(step_ratio, parity, nterms): + ''' Return matrix for finite difference and complex step derivation. Parameters ---------- + step_ratio : real scalar + ratio between steps in unequally spaced difference rule. parity : scalar, integer 0 (one sided, all terms included but zeroth order) 1 (only odd terms included) 2 (only even terms included) + 3 (only every 4'th order terms included starting from order 2) + 4 (only every 4'th order terms included starting from order 4) nterms : scalar, integer number of terms - - Member variables used - --------------------- - step_ratio ''' - srinv = 1.0 / self.step_ratio - [i, j] = np.ogrid[0:nterms, 0:nterms] - try: - fact = [1, 2, 2][parity] + step = [1, 2, 2, 4, 4, 4, 4][parity] except Exception as msg: - raise ValueError('%s. Parity must be 0, 1 or 2! (%d)' % (str(msg), - parity)) - offset = max(1, parity) - c = 1.0 / misc.factorial(np.arange(offset, fact * nterms + 1, fact)) - mat = c[j] * srinv ** (i * (fact * j + offset)) - return np.matrix(mat) - - def _set_fd_rule(self): + raise ValueError('%s. Parity must be 0, 1, 2, 3, 4, 5 or 6! ' + + '(%d)' % (str(msg), parity)) + inv_sr = 1.0 / step_ratio + offset = [1, 1, 2, 2, 4, 1, 3][parity] + c0 = [1.0, 1.0, 1.0, 2.0, 24.0, 1.0, 6.0][parity] + c = c0/misc.factorial(np.arange(offset, step * nterms + offset, step)) + [i, j] = np.ogrid[0:nterms, 0:nterms] + return np.atleast_2d(c[j] * inv_sr ** (i * (step * j + offset))) + + def _flip_fd_rule(self): + n = self.n + return ((self._is_even_derivative() and (self.method == 'backward')) or + (self.method == 'complex' and (n % 8 in [3, 4, 5, 6]))) + + def _get_finite_difference_rule(self, step_ratio): ''' Generate finite differencing rule in advance. @@ -701,7 +1126,7 @@ def _set_fd_rule(self): Member methods used ------------------- - _fd_mat + _fd_matrix Member variables used --------------------- @@ -709,445 +1134,149 @@ def _set_fd_rule(self): order method ''' + method = self.method + if method in ('multicomplex', ): + return np.ones((1,)) - order, method_order = self.n - 1, self.order + order, method_order = self.n - 1, self._method_order parity = 0 - num_terms = order + method_order - if self.method[0] == 'c': # 'central' + if (method.startswith('central') or + (method.startswith('complex') and self.n == 1 and + method_order < 4)): parity = (order % 2) + 1 - num_terms, order = num_terms // 2, order // 2 + elif self.method == 'complex': + if self._is_odd_derivative(): + parity = 6 if self._is_third_derivative() else 5 + else: + parity = 4 if self._is_fourth_derivative() else 3 - fd_rule = linalg.pinv(self._fd_mat(parity, num_terms))[order] - self._fd_rule = np.matrix(fd_rule) + step = self._richardson_step() + num_terms, ix = (order + method_order) // step, order // step + fd_mat = self._fd_matrix(step_ratio, parity, num_terms) + fd_rule = linalg.pinv(fd_mat)[ix] - def _get_min_num_steps(self): - n0 = 5 if self.method[0] == 'c' else 7 - return int(n0 + np.ceil(self.n / 2.) + self.order + self.romberg_terms) + if self._flip_fd_rule(): + fd_rule *= -1 + return fd_rule - def _set_romb_qr(self): + def _apply_fd_rule(self, fd_rule, sequence, steps): ''' + Return derivative estimates of f at x0 for a sequence of stepsizes h + Member variables used - order - method - romberg_terms - ''' - num_terms = self.romberg_terms - rmat = np.ones((num_terms + 2, num_terms + 1)) - if num_terms > 0: - i, j = np.ogrid[0:num_terms + 2, 0:num_terms] - fact = dict(c=2).get(self.method[0], 1) - srinv = self._make_exact(1.0 / self.step_ratio) - rmat[:, 1:] = srinv ** (i * (fact*j + self.order)) -# add1 = self.method[0] == 'c' -# rombexpon = (1 + add1) * np.arange(num_terms) + self.order -# for n in range(1, num_terms + 2): -# rmat[n, 1:] = srinv ** (n * rombexpon) - rmat = np.matrix(rmat) - self._qromb, self._rromb = linalg.qr(rmat) - self._rmat = rmat - self._romberg_rule = linalg.pinv(rmat)[0] - - def _set_difference_function(self): - ''' Set _diff_fun function according to method + --------------------- + n ''' - vectorized = ' v'[int(self.vectorized)].rstrip() - even_order = 'oe'[(np.remainder(self.n, 2) == 0)] - code = (self.method[0] + even_order + vectorized) - self._diff_fun = dict(cov=self._central_odd_vectorized, - cev=self._central_even_vectorized, - bov=self._backward_vectorized, - bev=self._backward_vectorized, - fov=self._forward_vectorized, - fev=self._forward_vectorized, - co=self._central_odd, - ce=self._central_even, - bo=self._backward, - be=self._backward, - fo=self._forward, - fe=self._forward)[code] + f_del, h, original_shape = self._vstack(sequence, steps) + + ne = h.shape[0] + if ne < fd_rule.size: + raise ValueError('num_steps (%d) must be larger than ' + '(%d) n + order - 1 = %d + %d -1' + ' (%s)' % (ne, fd_rule.size, self.n, self.order, + self.method) + ) + nr = (fd_rule.size-1) + f_diff = convolve1d(f_del, fd_rule[::-1], axis=0, origin=nr//2) + + der_init = f_diff / (h ** self.n) + ne = max(ne - nr, 1) + return der_init[:ne], h[:ne], original_shape + + def _derivative(self, xi, args, kwds): + diff, f = self._get_functions() + steps = self._get_steps(xi) + fxi = self._eval_first(f, xi, *args, **kwds) + results = [diff(f, fxi, xi, h, *args, **kwds) for h in steps] + step_ratio = self._compute_step_ratio(steps) + + self._set_richardson_rule(step_ratio, self.richardson_terms) + fd_rule = self._get_finite_difference_rule(step_ratio) + return self._apply_fd_rule(fd_rule, results, steps) @staticmethod - def _central_even_vectorized(fun, f_x0i, x0i, h): - return (fun(x0i + h) + fun(x0i - h)).ravel() / 2.0 - f_x0i + def _central_even(fun, f_x0i, x0i, h, *args, **kwds): + return (fun(x0i + h, *args, **kwds) + + fun(x0i - h, *args, **kwds)) / 2.0 - f_x0i @staticmethod - def _central_odd_vectorized(fun, f_x0i, x0i, h): - return (fun(x0i + h) - fun(x0i - h)).ravel() / 2.0 + def _central(fun, f_x0i, x0i, h, *args, **kwds): + return (fun(x0i + h, *args, **kwds) - + fun(x0i - h, *args, **kwds)) / 2.0 @staticmethod - def _central_even(fun, f_x0i, x0i, h): - return np.asfarray([fun(x0i + h_j) + fun(x0i - h_j) - for h_j in h]).ravel() / 2.0 - f_x0i + def _forward(fun, f_x0i, x0i, h, *args, **kwds): + return (fun(x0i + h, *args, **kwds) - f_x0i) @staticmethod - def _central_odd(fun, f_x0i, x0i, h): - return np.asfarray([fun(x0i + h_j) - fun(x0i - h_j) - for h_j in h]).ravel() / 2.0 + def _backward(fun, f_x0i, x0i, h, *args, **kwds): + return (f_x0i - fun(x0i - h)) @staticmethod - def _forward_vectorized(fun, f_x0i, x0i, h): - return (fun(x0i + h) - f_x0i).ravel() + def _complex(f, fx, x, h, *args, **kwds): + return f(x + 1j * h, *args, **kwds).imag @staticmethod - def _forward(fun, f_x0i, x0i, h): - return np.asfarray([fun(x0i + h_j) - f_x0i for h_j in h]).ravel() + def _complex_odd(f, fx, x, h, *args, **kwds): + ih = h * _SQRT_J + return ((_SQRT_J/2.) * (f(x + ih, *args, **kwds) - + f(x - ih, *args, **kwds))).imag @staticmethod - def _backward_vectorized(fun, f_x0i, x0i, h): - return (f_x0i - fun(x0i - h)).ravel() + def _complex_odd_higher(f, fx, x, h, *args, **kwds): + ih = h * _SQRT_J + return ((3 * _SQRT_J) * (f(x + ih, *args, **kwds) - + f(x - ih, *args, **kwds))).real @staticmethod - def _backward(fun, f_x0i, x0i, h): - return np.asfarray([f_x0i - fun(x0i - h_j) for h_j in h]).ravel() - - def _remove_non_finite(self, der_init, h1): - isnonfinite = 1 - np.isfinite(der_init) - i_nonfinite, = isnonfinite.ravel().nonzero() - if i_nonfinite.size > 0: - i = np.max(i_nonfinite) - im = np.min(i_nonfinite) - if i < der_init.size - 1: - msg = 'The stepsize (%g) is possibly too large!' % h1[i] - der_init = der_init[i+1:] - h1 = h1[i+1:] - elif 0 < im: - msg = 'The stepsize (%g) is possibly too small!' % h1[im] - der_init = der_init[:im] - h1 = h1[:im] - else: - msg = '' - k, = (1-isnonfinite).nonzero() - der_init = der_init[k] - h1 = h1[k] - warnings.warn(msg) - return der_init, h1 - - def _predict_uncertainty(self, rombcoefs): - '''uncertainty estimate of derivative prediction''' - coefs = rombcoefs[0][0] - s = np.sqrt(rombcoefs[1]) - rinv = np.asarray(linalg.pinv(self._rromb)) - cov1 = np.sum(rinv ** 2, axis=1) # 1 spare dof - errest = np.maximum(s * 12.7062047361747 * np.sqrt(cov1[0]), - s * _EPS * 10.) - - tmp_err = np.abs(np.diff(np.hstack((0, coefs, 0)))) * (len(coefs) > 1) - abserr = tmp_err[:-1] + tmp_err[1:] + np.abs(coefs) * _EPS * 10.0 - return np.maximum(errest, abserr) - - def _romb_extrap(self, der_init, h1): - ''' Return Romberg extrapolated derivatives and error estimates - based on the initial derivative estimates - - Parameter - --------- - der_init - initial derivative estimates - h1 - stepsizes used in the derivative estimates - - Returns - ------- - der_romb - derivative estimates returned - errest - error estimates - hout - stepsizes returned - - Member variables used - --------------------- - step_ratio - Ratio decrease in step - romberg_terms - higher order terms to cancel using the romberg step - ''' - # amp = np.linalg.cond(self._rromb) - # amp - noise amplification factor due to the romberg step - # the noise amplification is further amplified by the Romberg step - der_romb, hout = self._remove_non_finite(der_init, h1) - # this does the extrapolation to a zero step size. - num_terms = self.romberg_terms - ne = der_romb.size - - if ne < num_terms + 2: - errest = np.ones(der_init.shape) * hout - else: - rr_rule = self._romberg_rule - der_romb1 = convolve1d(der_romb, rr_rule[::-1], - origin=(rr_rule.size-1)//2) - rhs = vec2mat(der_romb, num_terms + 2, max(1, ne - num_terms - 2)) - - rombcoefs = linalg.lstsq(self._rromb, (self._qromb.T * rhs)) - der_romb = rombcoefs[0][0, :] - hout = hout[:der_romb.size] - errest = self._predict_uncertainty(rombcoefs) - - if self.use_dea and der_romb.size > 2: - der_romb, errest = dea3(der_romb[0:-2], der_romb[1:-1], der_romb[2:], symmetric=False) - hout = hout[2:] - - return der_romb, errest + _EPS, hout - - -class _PartialDerivative(_Derivative): - - def _partial_der(self, x00): - ''' Return partial derivatives - ''' - x0 = np.atleast_1d(x00) - nx = len(x0) - df, err, delta = np.zeros(nx), np.zeros(nx), np.zeros(nx) - - step_nom = [None, ] * nx if self.step_nom is None else self.step_nom - - fun = self._fun - self._x = np.asarray(x0, dtype=float) - for i in range(nx): - self._ix = i - df[i], err[i], delta[i] = self._derivative(fun, x0[i], step_nom[i]) - - return df, err, delta - - def _fun(self, xi): - x = self._x.copy() - x[self._ix] = xi - return self.fun(x) - - -class Derivative(_Derivative): - __doc__ = '''Estimate n'th derivative of fun at x0, with error estimate - - %s - - Examples - -------- - >>> import numpy as np - >>> import numdifftools as nd - - # 1'st and 2'nd derivative of exp(x), at x == 1 - - >>> fd = nd.Derivative(np.exp) # 1'st derivative - >>> fdd = nd.Derivative(np.exp,n=2) # 2'nd derivative - >>> fd(1) - array([ 2.71828183]) - - >>> d2 = fdd([1, 2]) - >>> d2 - array([ 2.71828183, 7.3890561 ]) - - >>> np.abs(d2-np.exp([1,2]))< fdd.error_estimate # Check error estimate - array([ True, True], dtype=bool) - - # 3'rd derivative of x.^3+x.^4, at x = [0,1] - - >>> fun = lambda x: x**3 + x**4 - >>> dfun = lambda x: 6 + 4*3*2*np.asarray(x) - >>> fd3 = nd.Derivative(fun,n=3) - >>> fd3([0,1]) # True derivatives: [6,30] - array([ 6., 30.]) + def _complex_even(f, fx, x, h, *args, **kwargs): + ih = h * _SQRT_J + return (f(x + ih, *args, **kwargs) + + f(x - ih, *args, **kwargs)).imag - >>> np.abs(fd3([0,1])-dfun([0,1])) <= fd3.error_estimate - array([ True, True], dtype=bool) - - See also - -------- - Gradient, - Hessdiag, - Hessian, - Jacobian - ''' % _Derivative.__doc__.partition('\n')[2] if _Derivative.__doc__ else '' - - def __call__(self, x): - return self.derivative(x) - - def _get_transformed_fun(self, x): - if self.transform is None: - fun = self.fun - f = 1 - elif self.transform == 'exp': - fun = lambda x: np.exp(self.fun(x)) - f = np.exp(-self.fun(x)) - elif self.transform == 'log': - fun = lambda x: np.log(np.abs(self.fun(x))) - f = self.fun(x) - return fun, f - - def derivative(self, x): - ''' Return estimate of n'th derivative of fun at x - using romberg extrapolation - ''' - self._initialize() - x0 = np.atleast_1d(x) - shape = x0.shape - fun, f0 = self._get_transformed_fun(x0) - der, err, delta = self._derivative(fun, x0.ravel(), self.step_nom) - self.error_estimate = err.reshape(shape) * f0 - self.final_delta = delta.reshape(shape) - return der.reshape(shape) * f0 - - -class Jacobian(_Derivative): - __doc__ = ('''Estimate Jacobian matrix, with error estimate - %s - - The Jacobian matrix is the matrix of all first-order partial derivatives - of a vector-valued function. - - Assumptions - ----------- - fun : (vector valued) - analytical function to differentiate. - fun must be a function of the vector or array x0. - - x0 : vector location at which to differentiate fun - If x0 is an N x M array, then fun is assumed to be - a function of N*M variables. - - Examples - -------- - >>> import numpy as np - >>> import numdifftools as nd - - #(nonlinear least squares) - - >>> xdata = np.reshape(np.arange(0,1,0.1),(-1,1)) - >>> ydata = 1+2*np.exp(0.75*xdata) - >>> fun = lambda c: (c[0]+c[1]*np.exp(c[2]*xdata) - ydata)**2 - >>> Jfun = nd.Jacobian(fun) - >>> h = Jfun([1., 2., 0.75]) # should be numerically zero - >>> np.abs(h) < 1e-14 - array([[ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True]], dtype=bool) - - >>> np.abs(h) <= 100 * Jfun.error_estimate - array([[ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True], - [ True, True, True]], dtype=bool) - - See also - -------- - Gradient, - Derivative, - Hessdiag, - Hessian - ''' % _Derivative.__doc__.partition('\n')[2].replace( - 'Integer from 1 to 4 (Default 1)', '1').replace( - 'defining derivative order.', - 'Derivative order is always 1.') if _Derivative.__doc__ else '') - - def __call__(self, x): - return self.jacobian(x) - - def jacobian(self, x): - ''' - Return Jacobian matrix of a vector valued function of n variables - - - Parameter - --------- - x : vector - location at which to differentiate fun. - If x is an nxm array, then fun is assumed to be - a function of n*m variables. - - Member variable used - -------------------- - fun : (vector valued) analytical function to differentiate. - fun must be a function of the vector or array x0. - - Returns - ------- - jac : array-like - first partial derivatives of fun. Assuming that x0 - is a vector of length p and fun returns a vector - of length n, then jac will be an array of size (n,p) - - err - vector - of error estimates corresponding to each partial - derivative in jac. - - See also - -------- - Derivative, - Gradient, - Hessian, - Hessdiag - ''' - self.n = 1 - fun = self.fun - self._initialize() - - zeros = np.zeros - newaxis = np.newaxis - x0 = np.atleast_1d(x) - nx = x0.size - - f0 = fun(x0) - f0 = f0.ravel() - n = f0.size - - jac = zeros((n, nx)) - if n == 0: - self.error_estimate = jac - return jac - - step_nom = self._get_step_nom(self.step_nom, x0) - - err, delta = jac.copy(), jac.copy() - for i in range(nx): - x0_i = x0[i] - h = self._get_steps(step_nom[i]) - nsteps = h.size - # evaluate at each step, centered around x0_i - # difference to give a second order estimate - fdel = zeros((n, nsteps)) - xp, xm = x0.copy(), x0.copy() - for j in range(nsteps): - xp[i], xm[i] = x0_i + h[j], x0_i - h[j] - fdif = fun(xp) - fun(xm) - fdel[:, j] = 0.5 * fdif.ravel() - derest = fdel / h[newaxis, :] - - for j in range(n): - der_romb, errors, h1 = self._romb_extrap(derest[j, :], h) - jac[j, i], err[j, i], delta[j, i] = self._best_der(der_romb, - errors, h1) - - self.final_delta = delta - self.error_estimate = err - return jac - - -class Gradient(_PartialDerivative): - __doc__ = ('''Estimate gradient of fun at x, with error estimate - %s - - Assumptions - ----------- - fun : SCALAR analytical function to differentiate. - fun must be a function of the vector or array x, - but it needs not to be vectorized. - - x : vector location at which to differentiate fun - If x is an N x M array, then fun is assumed to be - a function of N*M variables. + @staticmethod + def _complex_even_higher(f, fx, x, h, *args, **kwargs): + ih = h * _SQRT_J + return 12.0 * (f(x + ih, *args, **kwargs) + + f(x - ih, *args, **kwargs) - 2 * fx).real + @staticmethod + def _multicomplex(f, fx, x, h, *args, **kwds): + z = bicomplex(x + 1j * h, 0) + return f(z, *args, **kwds).imag + @staticmethod + def _multicomplex2(f, fx, x, h, *args, **kwds): + z = bicomplex(x + 1j * h, h) + return f(z, *args, **kwds).imag12 + + +class Gradient(Derivative): + def __init__(self, f, step=None, method='central', order=2, + full_output=False): + super(Gradient, self).__init__(f, step=step, method=method, n=1, + order=order, full_output=full_output) + __doc__ = _cmn_doc % dict( + derivative='Gradient', + extra_parameter="""order : int, optional + defines the order of the error term in the Taylor approximation used. + For 'central' and 'complex' methods, it must be an even number.""", + returns=""" + Returns + ------- + grad : array + gradient + """, extra_note=""" + Higher order approximation methods will generally be more accurate, but may + also suffer more from numerical problems. First order methods is usually + not recommended. + """, example=""" Examples -------- >>> import numpy as np - >>> import numdifftools as nd + >>> import numdifftools.nd_cstep as ndc >>> fun = lambda x: np.sum(x**2) - >>> dfun = nd.Gradient(fun) + >>> dfun = ndc.Gradient(fun) >>> dfun([1,2,3]) array([ 2., 4., 6.]) @@ -1156,7 +1285,7 @@ class Gradient(_PartialDerivative): >>> sin = np.sin; exp = np.exp >>> z = lambda xy: sin(xy[0]-xy[1]) + xy[1]*exp(xy[0]) - >>> dz = nd.Gradient(z) + >>> dz = ndc.Gradient(z) >>> grad2 = dz([1, 1]) >>> grad2 array([ 3.71828183, 1.71828183]) @@ -1165,285 +1294,578 @@ class Gradient(_PartialDerivative): # compute the gradient. It should be essentially zero. >>> rosen = lambda x : (1-x[0])**2 + 105.*(x[1]-x[0]**2)**2 - >>> rd = nd.Gradient(rosen) + >>> rd = ndc.Gradient(rosen) >>> grad3 = rd([1,1]) - >>> np.abs(grad3 - 0)<=rd.error_estimate - array([ True, True], dtype=bool) + >>> np.allclose(grad3,[0, 0]) + True""", see_also=""" + See also + -------- + Derivative, Hessian, Jacobian + """) + @staticmethod + def _central(f, fx, x, h, *args, **kwds): + n = len(x) + increments = np.identity(n) * h + partials = [(f(x + hi, *args, **kwds) - f(x - hi, *args, **kwds)) / 2.0 + for hi in increments] + return np.array(partials).T - See also + @staticmethod + def _backward(f, fx, x, h, *args, **kwds): + n = len(x) + increments = np.identity(n) * h + partials = [(fx - f(x - hi, *args, **kwds)) for hi in increments] + return np.array(partials).T + + @staticmethod + def _forward(f, fx, x, h, *args, **kwds): + n = len(x) + increments = np.identity(n) * h + partials = [(f(x + hi, *args, **kwds) - fx) for hi in increments] + return np.array(partials).T + + @staticmethod + def _complex(f, fx, x, h, *args, **kwds): + # From Guilherme P. de Freitas, numpy mailing list + # http://mail.scipy.org/pipermail/numpy-discussion/2010-May/050250.html + n = len(x) + increments = np.identity(n) * 1j * h + partials = [f(x + ih, *args, **kwds).imag for ih in increments] + return np.array(partials).T + + @staticmethod + def _complex_odd(f, fx, x, h, *args, **kwds): + n = len(x) + increments = np.identity(n) * _SQRT_J * h + + partials = [((_SQRT_J/2.) * (f(x + ih, *args, **kwds) - + f(x - ih, *args, **kwds))).imag + for ih in increments] + return np.array(partials).T + + @staticmethod + def _multicomplex(f, fx, x, h, *args, **kwds): + n = len(x) + increments = np.identity(n) * 1j * h + partials = [f(bicomplex(x + hi, 0), *args, **kwds).imag + for hi in increments] + return np.array(partials).T + + +class Jacobian(Gradient): + __doc__ = _cmn_doc % dict( + derivative='Jacobian', + extra_parameter="""order : int, optional + defines the order of the error term in the Taylor approximation used. + For 'central' and 'complex' methods, it must be an even number.""", + returns=""" + Returns + ------- + jacob : array + Jacobian + """, extra_note=""" + Higher order approximation methods will generally be more accurate, but may + also suffer more from numerical problems. First order methods is usually + not recommended. + + If f returns a 1d array, it returns a Jacobian. If a 2d array is returned + by f (e.g., with a value for each observation), it returns a 3d array + with the Jacobian of each observation with shape xk x nobs x xk. I.e., + the Jacobian of the first observation would be [:, 0, :] + """, example=''' + Examples -------- - Derivative, Hessdiag, Hessian, Jacobian - ''' % _Derivative.__doc__.partition('\n')[2].replace( - 'Integer from 1 to 4 (Default 1)', '1').replace( - 'defining derivative order.', - 'Derivative order is always 1.') if _Derivative.__doc__ else '') + >>> import numdifftools.nd_cstep as ndc - def __call__(self, x): - return self.gradient(x) + #(nonlinear least squares) - def gradient(self, x): - '''Returns gradient + >>> xdata = np.reshape(np.arange(0,1,0.1),(-1,1)) + >>> ydata = 1+2*np.exp(0.75*xdata) + >>> fun = lambda c: (c[0]+c[1]*np.exp(c[2]*xdata) - ydata)**2 - See also derivative, hessian, jacobian - ''' - self.n = 1 - self.vectorized = False - self._initialize() - pder, self.error_estimate, self.final_delta = self._partial_der(x) - return pder - - -class Hessdiag(_PartialDerivative): - __doc__ = (''' - Estimate diagonal elements of Hessian of fun at x, with error estimate - %s - - HESSDIAG return a vector of second order partial derivatives of fun. - These are the diagonal elements of the Hessian matrix, evaluated - at x0. When all that you want are the diagonal elements of the hessian - matrix, it will be more efficient to call HESSDIAG than HESSIAN. - HESSDIAG uses DERIVATIVE to provide both second derivative estimates - and error estimates. - - Assumptions - ------------ - fun : SCALAR analytical function to differentiate. - fun must be a function of the vector or array x0, - but it needs not to be vectorized. - - x : vector location at which to differentiate fun - If x is an N x M array, then fun is assumed to be - a function of N*M variables. + >>> Jfun = ndc.Jacobian(fun) + >>> val = Jfun([1,2,0.75]) + >>> np.allclose(val, np.zeros((10,3))) + True + >>> fun2 = lambda x : x[0]*x[1]*x[2] + np.exp(x[0])*x[1] + >>> Jfun3 = ndc.Jacobian(fun2) + >>> Jfun3([3.,5.,7.]) + array([ 135.42768462, 41.08553692, 15. ]) + ''', see_also=""" + See also + -------- + Derivative, Hessian, Gradient + """) + + +class Hessdiag(Derivative): + def __init__(self, f, step=None, method='central', order=2, + full_output=False): + super(Hessdiag, self).__init__(f, step=step, method=method, n=2, + order=order, full_output=full_output) + __doc__ = _cmn_doc % dict( + derivative='Hessian diagonal', + extra_parameter="""order : int, optional + defines the order of the error term in the Taylor approximation used. + For 'central' and 'complex' methods, it must be an even number.""", + returns=""" + Returns + ------- + hessdiag : array + hessian diagonal + """, extra_note=""" + Higher order approximation methods will generally be more accurate, but may + also suffer more from numerical problems. First order methods is usually + not recommended. + """, example=""" Examples -------- >>> import numpy as np - >>> import numdifftools as nd + >>> import numdifftools.nd_cstep as ndc >>> fun = lambda x : x[0] + x[1]**2 + x[2]**3 - >>> ddfun = lambda x : np.asarray((0, 2, 6*x[2])) - >>> Hfun = nd.Hessdiag(fun) - >>> hd = Hfun([1,2,3]) # HD = [ 0,2,18] - >>> hd - array([ 0., 2., 18.]) - >>> np.abs(ddfun([1,2,3])-hd) <= Hfun.error_estimate - array([ True, True, True], dtype=bool) - + >>> Hfun = ndc.Hessdiag(fun, full_output=True) + >>> hd, info = Hfun([1,2,3]) + >>> np.allclose(hd, [ 0., 2., 18.]) + True + >>> info.error_estimate < 1e-11 + array([ True, True, True], dtype=bool) + """, see_also=""" See also -------- - Gradient, Derivative, Hessian, Jacobian - ''' % _Derivative.__doc__.partition('\n')[2].replace( - 'Integer from 1 to 4 (Default 1)', '2').replace( - 'defining derivative order.', - 'Derivative order is always 2.') if _Derivative.__doc__ else '') - - def __call__(self, x): - return self.hessdiag(x) - - def hessdiag(self, x): - ''' Diagonal elements of Hessian matrix - - See also derivative, gradient, hessian, jacobian - ''' - self.n = 2 - self.vectorized = False - self._initialize() - dder, self.error_estimate, self.final_delta = self._partial_der(x) - return dder + Derivative, Hessian, Jacobian, Gradient + """) + @staticmethod + def _central2(f, fx, x, h, *args, **kwds): + '''Eq. 8''' + n = len(x) + increments = np.identity(n) * h + partials = [(f(x + 2*hi, *args, **kwds) + + f(x - 2*hi, *args, **kwds) + 2*fx - + 2*f(x + hi, *args, **kwds) - + 2*f(x - hi, *args, **kwds)) / 4.0 + for hi in increments] + return np.array(partials) -class Hessian(Hessdiag): - __doc__ = (''' Estimate Hessian matrix, with error estimate - %s + @staticmethod + def _central_even(f, fx, x, h, *args, **kwds): + '''Eq. 9''' + n = len(x) + increments = np.identity(n) * h + partials = [(f(x + hi, *args, **kwds) + + f(x - hi, *args, **kwds)) / 2.0 - fx + for hi in increments] + return np.array(partials) - HESSIAN estimate the matrix of 2nd order partial derivatives of a real - valued function FUN evaluated at X. HESSIAN is NOT a tool for frequent - use on an expensive to evaluate objective function, especially in a large - number of dimensions. Its computation will use roughly O(6*n^2) function - evaluations for n parameters. + @staticmethod + def _backward(f, fx, x, h, *args, **kwds): + n = len(x) + increments = np.identity(n) * h + partials = [(fx - f(x - hi, *args, **kwds)) for hi in increments] + return np.array(partials) - Assumptions - ----------- - fun : SCALAR analytical function - to differentiate. fun must be a function of the vector or array x0, - but it needs not to be vectorized. + @staticmethod + def _forward(f, fx, x, h, *args, **kwds): + n = len(x) + increments = np.identity(n) * h + partials = [(f(x + hi, *args, **kwds) - fx) for hi in increments] + return np.array(partials) - x : vector location - at which to differentiate fun - If x is an N x M array, then fun is assumed to be a function - of N*M variables. + @staticmethod + def _multicomplex2(f, fx, x, h, *args, **kwds): + n = len(x) + increments = np.identity(n) * h + partials = [f(bicomplex(x + 1j * hi, hi), *args, **kwds).imag12 + for hi in increments] + return np.array(partials) + @staticmethod + def _complex_even(f, fx, x, h, *args, **kwargs): + n = len(x) + increments = np.identity(n) * h * (1j+1) / np.sqrt(2) + partials = [(f(x + hi, *args, **kwargs) + + f(x - hi, *args, **kwargs)).imag + for hi in increments] + return np.array(partials) + + +class Hessian(_Derivative): + def __init__(self, f, step=None, method='central', full_output=False): + order = dict(backward=1, forward=1, complex=2).get(method, 2) + super(Hessian, self).__init__(f, n=2, step=step, method=method, + order=order, full_output=full_output) + + __doc__ = _cmn_doc % dict( + derivative='Hessian', + extra_parameter="", + returns=""" + Returns + ------- + hess : ndarray + array of partial second derivatives, Hessian + """, extra_note=""" + Computes the Hessian according to method as: + 'forward', Eq. (7): + 1/(d_j*d_k) * ((f(x + d[j]*e[j] + d[k]*e[k]) - f(x + d[j]*e[j]))) + 'central', Eq. (9): + 1/(4*d_j*d_k) * ((f(x + d[j]*e[j] + d[k]*e[k]) - + f(x + d[j]*e[j] - d[k]*e[k])) - + (f(x - d[j]*e[j] + d[k]*e[k]) - + f(x - d[j]*e[j] - d[k]*e[k])) + 'complex', Eq. (10): + 1/(2*d_j*d_k) * imag(f(x + i*d[j]*e[j] + d[k]*e[k]) - + f(x + i*d[j]*e[j] - d[k]*e[k])) + where e[j] is a vector with element j == 1 and the rest are zero and + d[i] is steps[i]. + """, example=""" Examples -------- >>> import numpy as np - >>> import numdifftools as nd + >>> import numdifftools.nd_cstep as ndc # Rosenbrock function, minimized at [1,1] >>> rosen = lambda x : (1.-x[0])**2 + 105*(x[1]-x[0]**2)**2 - >>> Hfun = nd.Hessian(rosen) + >>> Hfun = ndc.Hessian(rosen) >>> h = Hfun([1, 1]) >>> h array([[ 842., -420.], [-420., 210.]]) - >>> Hfun.error_estimate < 5.e-10 - array([[ True, True], - [ True, True]], dtype=bool) # cos(x-y), at (0,0) >>> cos = np.cos >>> fun = lambda xy : cos(xy[0]-xy[1]) - >>> Hfun2 = nd.Hessian(fun) + >>> Hfun2 = ndc.Hessian(fun) >>> h2 = Hfun2([0, 0]) >>> h2 array([[-1., 1.], - [ 1., -1.]]) - >>> np.abs(h2-np.array([[-1, 1],[ 1, -1]])) < 10*Hfun2.error_estimate - array([[ True, True], - [ True, True]], dtype=bool) - - >>> Hfun2.romberg_terms = 3 - >>> h3 = Hfun2([0,0]) - >>> h3 - array([[-1., 1.], - [ 1., -1.]]) - >>> np.abs(h3-np.array([[-1., 1.],[ 1., -1.]])) < 10*Hfun2.error_estimate - array([[ True, True], - [ True, True]], dtype=bool) - - + [ 1., -1.]])""", see_also=""" See also -------- - Gradient, - Derivative, - Hessdiag, - Jacobian - ''' % _Derivative.__doc__.partition('\n')[2].replace( - 'Integer from 1 to 4 (Default 1)', '2').replace( - 'defining derivative order.', - 'Derivative order is always 2.') if _Derivative.__doc__ else '') - - def __call__(self, x): - return self.hessian(x) - - def _get_step_max(self): - # Decide on intelligent step sizes for the mixed partials - best_step_size = self.final_delta - num_steps = self._get_min_num_steps() - if self.step_num is not None: - num_steps = min(num_steps, self.step_num) - deltas = (1.0 * self.step_ratio) ** (-np.arange(num_steps)) - deltas = self._make_exact(deltas) - stepmax = best_step_size / deltas[num_steps // 2] - return stepmax, deltas - - def hessian(self, x): - '''Hessian matrix i.e., array of 2nd order partial derivatives - - See also derivative, gradient, hessdiag, jacobian + Derivative, Hessian + """) + + def _complex_high_order(self): + return False + + def _derivative(self, xi, args, kwds): + diff, f = self._get_functions() + steps = self._get_steps(xi) + + fxi = self._eval_first(f, xi, *args, **kwds) + results = [diff(f, fxi, xi, h, *args, **kwds) for h in steps] + step_ratio = self._compute_step_ratio(steps) + self._set_richardson_rule(step_ratio, self.richardson_terms) + return self._vstack(results, steps) + + @staticmethod + def _complex_even(f, fx, x, h, *args, **kwargs): + '''Calculate Hessian with complex-step derivative approximation + The stepsize is the same for the complex and the finite difference part ''' - x0 = np.atleast_1d(x) - nx = len(x0) - self.method = 'central' - - hess = self.hessdiag(x0) - err = self.error_estimate - - hess, err = np.diag(hess), np.diag(err) - if nx < 2: - return hess # the hessian matrix is 1x1. all done - - stepmax, dfac = self._get_step_max() - ndel = dfac.size - fun = self.fun - zeros = np.zeros - for i in range(1, nx): - for j in range(i): - dij, step = zeros(ndel), zeros(nx) - step[[i, j]] = stepmax[[i, j]] - for k in range(ndel): - x1 = x0 + step * dfac[k] - x2 = x0 - step * dfac[k] - step[j] = -step[j] - x3 = x0 + step * dfac[k] - step = -step - x4 = x0 + step * dfac[k] - step[i] = -step[i] - dij[k] = fun(x1) + fun(x2) - fun(x3) - fun(x4) - h2 = stepmax[[i, j]].prod() * (dfac ** 2) - dij = dij / (4 * h2) - - hess_romb, errors, h = self._romb_extrap(dij, np.sqrt(h2)) - hess[i, j], err[i, j], h = self._best_der(hess_romb, errors, h) - hess[j, i], err[j, i] = hess[i, j], err[i, j] - - self.error_estimate = err + n = len(x) + # h = _default_base_step(x, 3, base_step, n) + ee = np.diag(h) + hes = 2. * np.outer(h, h) + + for i in range(n): + for j in range(i, n): + hes[i, j] = (f(x + 1j * ee[i] + ee[j], *args, **kwargs) - + f(x + 1j * ee[i] - ee[j], *args, **kwargs) + ).imag / hes[j, i] + hes[j, i] = hes[i, j] + return hes + + @staticmethod + def _multicomplex2(f, fx, x, h, *args, **kwargs): + '''Calculate Hessian with bicomplex-step derivative approximation + ''' + n = len(x) + ee = np.diag(h) + hess = np.outer(h, h) + for i in range(n): + for j in range(i, n): + zph = bicomplex(x + 1j * ee[i, :], ee[j, :]) + hess[i, j] = (f(zph, *args, **kwargs)).imag12 / hess[j, i] + hess[j, i] = hess[i, j] + return hess + + @staticmethod + def _central_even(f, fx, x, h, *args, **kwargs): + '''Eq 9.''' + n = len(x) + # h = _default_base_step(x, 4, base_step, n) + ee = np.diag(h) + hess = np.outer(h, h) + + for i in range(n): + hess[i, i] = (f(x + 2*ee[i, :], *args, **kwargs) - 2*fx + + f(x - 2*ee[i, :], *args, **kwargs) + ) / (4. * hess[i, i]) + for j in range(i+1, n): + hess[i, j] = (f(x + ee[i, :] + ee[j, :], *args, **kwargs) - + f(x + ee[i, :] - ee[j, :], *args, **kwargs) - + f(x - ee[i, :] + ee[j, :], *args, **kwargs) + + f(x - ee[i, :] - ee[j, :], *args, **kwargs) + ) / (4. * hess[j, i]) + hess[j, i] = hess[i, j] + return hess + + @staticmethod + def _central2(f, fx, x, h, *args, **kwargs): + '''Eq. 8''' + n = len(x) + # NOTE: ridout suggesting using eps**(1/4)*theta + # h = _default_base_step(x, 3, base_step, n) + ee = np.diag(h) + dtype = np.result_type(fx) + g = np.empty(n, dtype=dtype) + gg = np.empty(n, dtype=dtype) + for i in range(n): + g[i] = f(x + ee[i], *args, **kwargs) + gg[i] = f(x - ee[i], *args, **kwargs) + + hess = np.empty((n, n), dtype=dtype) + np.outer(h, h, out=hess) + for i in range(n): + for j in range(i, n): + hess[i, j] = (f(x + ee[i, :] + ee[j, :], *args, **kwargs) - + g[i] - g[j] + fx + + f(x - ee[i, :] - ee[j, :], *args, **kwargs) - + gg[i] - gg[j] + fx) / (2 * hess[j, i]) + hess[j, i] = hess[i, j] + + return hess + + @staticmethod + def _forward(f, fx, x, h, *args, **kwargs): + '''Eq. 7''' + n = len(x) + ee = np.diag(h) + + dtype = np.result_type(fx) + g = np.empty(n, dtype=dtype) + for i in range(n): + g[i] = f(x + ee[i, :], *args, **kwargs) + + hess = np.empty((n, n), dtype=dtype) + np.outer(h, h, out=hess) + for i in range(n): + for j in range(i, n): + hess[i, j] = (f(x + ee[i, :] + ee[j, :], *args, **kwargs) - + g[i] - g[j] + fx) / hess[j, i] + hess[j, i] = hess[i, j] return hess + def _backward(self, f, fx, x, h, *args, **kwargs): + return self._forward(f, fx, x, -h, *args, **kwargs) + + +def main(): + import statsmodels.api as sm + + data = sm.datasets.spector.load() + data.exog = sm.add_constant(data.exog, prepend=False) + mod = sm.Probit(data.endog, data.exog) + _res = mod.fit(method="newton") + _test_params = [1, 0.25, 1.4, -7] + _llf = mod.loglike + _score = mod.score + _hess = mod.hessian + + def fun(beta, x): + return np.dot(x, beta).sum(0) + + def fun1(beta, y, x): + # print(beta.shape, x.shape) + xb = np.dot(x, beta) + return (y - xb) ** 2 # (xb-xb.mean(0))**2 + + def fun2(beta, y, x): + # print(beta.shape, x.shape) + return fun1(beta, y, x).sum(0) + + nobs = 200 + x = np.random.randn(nobs, 3) + + # xk = np.array([1, 2, 3]) + xk = np.array([1., 1., 1.]) + # xk = np.zeros(3) + beta = xk + y = np.dot(x, beta) + 0.1 * np.random.randn(nobs) + xk = np.dot(np.linalg.pinv(x), y) + + epsilon = 1e-6 + args = (y, x) + from scipy import optimize + _xfmin = optimize.fmin(fun2, (0, 0, 0), args) # @UndefinedVariable + # print(approx_fprime((1, 2, 3), fun, steps, x)) + jac = Gradient(fun1, epsilon, method='forward')(xk, *args) + jacmin = Gradient(fun1, -epsilon, method='forward')(xk, *args) + # print(jac) + print(jac.sum(0)) + print('\nnp.dot(jac.T, jac)') + print(np.dot(jac.T, jac)) + print('\n2*np.dot(x.T, x)') + print(2 * np.dot(x.T, x)) + jac2 = (jac + jacmin) / 2. + print(np.dot(jac2.T, jac2)) + + # he = approx_hess(xk,fun2,steps,*args) + print(Hessian(fun2, 1e-3, method='central2')(xk, *args)) + he = Hessian(fun2, method='central2')(xk, *args) + print('hessfd') + print(he) + print('base_step =', None) + print(he - 2 * np.dot(x.T, x)) + + for eps in [1e-3, 1e-4, 1e-5, 1e-6]: + print('eps =', eps) + print(Hessian(fun2, eps, method='central2')(xk, *args) - + 2 * np.dot(x.T, x)) + + hcs2 = Hessian(fun2, method='hybrid')(xk, *args) + print('hcs2') + print(hcs2 - 2 * np.dot(x.T, x)) + + hfd3 = Hessian(fun2, method='central')(xk, *args) + print('hfd3') + print(hfd3 - 2 * np.dot(x.T, x)) + + hfi = [] + epsi = np.array([1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]) * 10. + for eps in epsi: + h = eps * np.maximum(np.log1p(np.abs(xk)), 0.1) + hfi.append(Hessian(fun2, h, method='hybrid')(xk, *args)) + print('hfi, eps =', eps) + print(hfi[-1] - 2 * np.dot(x.T, x)) + + import numdifftools as nd + print('Dea3') + err = 1000 * np.ones(hfi[0].shape) + val = np.zeros(err.shape) + errt = [] + for i in range(len(hfi) - 2): + tval, terr = nd.dea3(hfi[i], hfi[i + 1], hfi[i + 2]) + errt.append(terr) + k = np.flatnonzero(terr < err) + if k.size > 0: + np.put(val, k, tval.flat[k]) + np.put(err, k, terr.flat[k]) + print(val - 2 * np.dot(x.T, x)) + print(err) + erri = [v.max() for v in errt] + + plt.loglog(epsi[1:-1], erri) + plt.show('hold') + hnd = nd.Hessian(lambda a: fun2(a, y, x)) + hessnd = hnd(xk) + print('numdiff') + print(hessnd - 2 * np.dot(x.T, x)) + # assert_almost_equal(hessnd, he[0]) + gnd = nd.Gradient(lambda a: fun2(a, y, x)) + _gradnd = gnd(xk) + + print(Derivative(np.cosh)(0)) + print(nd.Derivative(np.cosh)(0)) + + +def _example3(x=0.0001, fun_name='cos', epsilon=None, method='central', + scale=None, n=1, order=2): + fun0, dfun = get_test_function(fun_name, n) + if dfun is None: + return dict(n=n, order=order, method=method, fun=fun_name, + error=np.nan, scale=np.nan) + fd = Derivative(fun0, step=epsilon, method=method, n=n, order=order) + t = [] + scales = np.arange(1.0, 45, 0.25) + for scale in scales: + fd.step.scale = scale + try: + val = fd(x) + except Exception: + val = np.nan + t.append(val) + t = np.array(t) + tt = dfun(x) + relativ_error = np.abs(t - tt) / (np.maximum(np.abs(tt), 1)) + 1e-16 + + weights = np.ones((3,))/3 + relativ_error = convolve1d(relativ_error, weights) # smooth curve + + if np.isnan(relativ_error).all(): + return dict(n=n, order=order, method=method, fun=fun_name, + error=np.nan, scale=np.nan) + if True: # False: # + plt.semilogy(scales, relativ_error) + plt.vlines(default_scale(fd.method, n, order), + np.nanmin(relativ_error), 1) + plt.xlabel('scales') + plt.ylabel('Relative error') + txt = ['', "1'st", "2'nd", "3'rd", "4'th", "5'th", "6'th", + "7th"] + ["%d'th" % i for i in range(8, 25)] + + plt.title("The %s derivative of %s using %s, order=%d" % (txt[n], + fun_name, + method, + order)) + + plt.axis([min(scales), max(scales), np.nanmin(relativ_error), 1]) + plt.figure() + # plt.show('hold') + i = np.nanargmin(relativ_error) + return dict(n=n, order=order, method=method, fun=fun_name, + error=relativ_error[i], scale=scales[i]) + + +def _example2(x=0.0001, fun_name='inv', epsilon=None, method='central', + scale=None, n=1): + fun0, dfun = get_test_function(fun_name, n) + + fd = Derivative(fun0, step=epsilon, method=method, n=n) + t = [] + orders = n + (n % 2) + np.arange(0, 12, 2) + + for order in orders: + fd.order = order + fd.step.num_steps = n + order - 1 + t.append(fd(x)) + t = np.array(t) + tt = dfun(x) + plt.semilogy(orders, np.abs(t - tt) / (np.abs(tt) + 1e-17) + 1e-17) + + plt.show('hold') + -def _example(x=0.0001, fun_name='inv', n=1, method='central', step_max=100, - step_ratio=2, step_num=30, romberg_terms=2, use_dea=True, - transform=None): +def _example(x=0.0001, fun_name='inv', epsilon=None, method='central', + scale=None): ''' ''' - sinh, cosh, tanh = np.sinh, np.cosh, np.tanh - f_dic = dict(cos=(np.cos, - lambda x: -np.sin(x), - lambda x: -np.cos(x), - lambda x: np.sin(x), - lambda x: np.cos(x),), - tanh=(tanh, - lambda x: 1. / cosh(x) ** 2, - lambda x: -2 * sinh(x) / cosh(x) ** 3, - lambda x: 2.*(3 * tanh(x)**2 - 1) / cosh(x)**2, - lambda x: (6 + 4*sinh(x)*(cosh(x) - - 3*tanh(x))) / cosh(x)**4), - log=(np.log, - lambda x: 1. / x, - lambda x: -1. / x ** 2, - lambda x: 2. / x ** 3, - lambda x: -6. / x ** 4), - inv=(lambda x: 1. / x, - lambda x: -1. / x ** 2, - lambda x: 2. / x ** 3, - lambda x: -6. / x ** 4, - lambda x: 24. / x ** 5)) - funs = f_dic.get(fun_name) - fun0 = funs[0] - dfun = funs[n] - - h = 1e-4 - fd = Derivative(fun0, n=n, method=method, step_max=step_max, - step_ratio=step_ratio, offset=1, - step_num=step_num, verbose=True, vectorized=True, - romberg_terms=romberg_terms, use_dea=use_dea, - transform=transform) - - t = fd(x) - txt = '' - if n == 1: - txt = (' (f(x+h)-f(x))/h = %g\n' % - ((fun0(x + h) - fun0(x)) / h)) + fun0, dfun = get_test_function(fun_name) + + h = _default_base_step(x, scale=2, epsilon=None) # 1e-4 + + fd = Derivative(fun0, step=epsilon, method=method, full_output=True) + + t, res = fd(x) + + txt = (' (f(x+h)-f(x))/h = %g\n' % + ((fun0(x + h) - fun0(x)) / h)) + deltas = np.array([h for h in epsilon(x, fd.scale)]) + print((txt + ' true df(x) = %20.15g\n' + ' estimated df(x) = %20.15g\n' + - ' true err = %g\n err estimate = %g\n relative err = %g\n' + - ' final_delta = %g\n') % (dfun(x), t, dfun(x) - t, - fd.error_estimate, - fd.error_estimate / t, - fd.final_delta)) - plt.show('hold') - - -def _test_rosen(): - def rosen(x): - return (1.-x[0])**2 + 105*(x[1]-x[0]**2)**2 - Hfun = Hessian(rosen, delta=[1e-6], step_num=1) - print(Hfun([1, 1])) - print(Hfun.final_delta) - print(Hfun.error_estimate) + ' true err = %g\n err estimate = %g\n relative err = %g\n' + ' delta = %g\n') % (dfun(x), t, dfun(x) - t, + res.error_estimate, + res.error_estimate / t, + deltas.flat[res.index])) + # plt.show('hold') def test_docstrings(): @@ -1451,12 +1873,45 @@ def test_docstrings(): doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE) -if __name__ == '__main__': - # _test_rosen() - # test_dea() - # test_epsal() - _example(x=0.01, fun_name='inv', n=1, method='c', step_max=1., - step_ratio=2., step_num=15, romberg_terms=2, use_dea=True, - transform=None) +def main2(): + import pandas as pd + num_extrap = 0 + method = 'complex' + data = [] + for name in ['exp', 'expm1', 'sin', 'cos', 'square']: + # function_names[:-3]: + for order in range(4, 5, 1): + # order = 1 + for n in range(1, 16, 1): + num_steps = n + order - 1 + num_extrap + if method in ['central', 'complex']: + step = 2 + if (n > 1 or order >= 4) and method == 'complex': + step = 4 + num_steps = (n + order-1) // step + num_extrap + + step_ratio = 1.6 # 4**(1./n) + epsilon = MinStepGenerator(num_steps=num_steps, + step_ratio=step_ratio, + offset=0, use_exact_steps=True) + data.append(pd.DataFrame(_example3(x=0.7, fun_name=name, + epsilon=epsilon, + method=method, + scale=None, n=n, + order=order), + index=np.arange(1))) + df = pd.concat(data) + # sprint(df) + print(df.groupby(['n']).mean()) + print(np.diff(df.groupby(['n']).mean(), axis=0)) + plt.show('hold') + +if __name__ == '__main__': # pragma : no cover + test_docstrings() + # main() + # main2() - #test_docstrings() +# r = _example3(x=1, fun_name='sin', epsilon=None, method='complex', +# scale=None, n=4, order=2) +# print(r) +# plt.show('hold') diff --git a/numdifftools/info.py b/numdifftools/info.py index fac7234..86196c5 100644 --- a/numdifftools/info.py +++ b/numdifftools/info.py @@ -7,10 +7,11 @@ manner, coupled with a Richardson extrapolation methodology to provide a maximally accurate result. The user can configure many options like; changing the order of the method or the extrapolation, even allowing the user to specify -whether central, forward or backward differences are used. The methods provided +whether complex-step, central, forward or backward differences are used. The +methods provided are: -*Derivative:* Computate derivatives of order 1 through 4 on any scalar function. +*Derivative:* Computate derivatives of order 1 through 10 on any scalar function. *Gradient:* Computes the gradient vector of a scalar function of one or more variables. @@ -23,6 +24,7 @@ All of these methods also produce error estimates on the result. Documentation is at: http://numdifftools.readthedocs.org/ + Code and issue tracker is at https://github.com/pbrod/numdifftools . To test if the toolbox is working paste the following in an interactive diff --git a/numdifftools/multicomplex.py b/numdifftools/multicomplex.py index a59f4d8..b767453 100644 --- a/numdifftools/multicomplex.py +++ b/numdifftools/multicomplex.py @@ -33,6 +33,24 @@ _TINY = np.finfo(float).machar.tiny +def c_atan2(x, y): + a, b = np.real(x), np.imag(x) + c, d = np.real(y), np.imag(y) + return np.arctan2(a, c) + 1j * (c * b - a * d) / (a**2 + c**2) + + +def c_max(x, y): + return np.where(x.real < y.real, y, x) + + +def c_min(x, y): + return np.where(x.real > y.real, y, x) + + +def c_abs(z): + return np.where(np.real(z) >= 0, z, -z) + + class bicomplex(object): '''BICOMPLEX(z1, z2) Creates an instance of a bicomplex object. diff --git a/numdifftools/tests/test_hessian.py b/numdifftools/tests/test_hessian.py index fd2bd73..a9bf3ee 100644 --- a/numdifftools/tests/test_hessian.py +++ b/numdifftools/tests/test_hessian.py @@ -9,7 +9,7 @@ # import numdifftools.nd_algopy as algopy # import numdifftools.nd_scientific as scientific from scipy import linalg, optimize, constants - +from numdifftools.multicomplex import c_abs as abs _TINY = np.finfo(float).machar.tiny @@ -79,24 +79,25 @@ def _run_hamiltonian(verbose=True): # Important to restrict the step in order to avoid the discontinutiy at # x=[0,0] # hessian = nd.Hessian(c.potential, step_max=1.0, step_nom=np.abs(xopt)) - hessian = nd.Hessian(c.potential) + step = nd.MaxStepGenerator(step_max=1e-1, num_steps=26) + hessian = nd.Hessian(c.potential, step=step, method='complex', full_output=True) # hessian = algopy.Hessian(c.potential) # Does not work # hessian = scientific.Hessian(c.potential) # does not work - H = hessian(xopt) + H, info = hessian(xopt) true_H = np.array([[5.23748385e-12, -2.61873829e-12], [-2.61873829e-12, 5.23748385e-12]]) if verbose: print(xopt) print('H', H) print('H-true_H', np.abs(H-true_H)) - print('error_estimate', hessian.error_estimate) + print('error_estimate', info.error_estimate) eigenvalues = linalg.eigvals(H) normal_modes = c.normal_modes(eigenvalues) print('eigenvalues', eigenvalues) print('normal_modes', normal_modes) - return H, hessian.error_estimate, true_H + return H, info.error_estimate, true_H class TestHessian(unittest.TestCase): diff --git a/numdifftools/tests/test_numdifftools.py b/numdifftools/tests/test_numdifftools.py index 0042fae..7ff0b51 100644 --- a/numdifftools/tests/test_numdifftools.py +++ b/numdifftools/tests/test_numdifftools.py @@ -2,35 +2,357 @@ """ import unittest -import numdifftools as nd +import numdifftools.core as nd import numpy as np from numpy.testing import assert_array_almost_equal +class TestGlobalFunctions(unittest.TestCase): + + def testdea3(self): + def linfun(k): + return np.linspace(0, np.pi / 2., 2. ** (k + 5) + 1) + Ei = np.zeros(3) + for k in np.arange(3): + x = linfun(k) + Ei[k] = np.trapz(np.sin(x), x) + [En, err] = nd.dea3(Ei[0], Ei[1], Ei[2]) + self.assertTrue(np.abs(En - 1) < err) + assert_array_almost_equal(En, 1.0, decimal=8) + + +class TestRichardson(unittest.TestCase): + + def test_order_step_combinations(self): + true_vals = { + (1, 1, 1): [-0.9999999999999998, 1.9999999999999998], + (1, 1, 2): [-0.33333333333333304, 1.333333333333333], + (1, 1, 3): [-0.14285714285714307, 1.142857142857143], + (1, 1, 4): [-0.06666666666666654, 1.0666666666666664], + (1, 1, 5): [-0.03225806451612906, 1.0322580645161292], + (1, 1, 6): [-0.015873015873015872, 1.0158730158730154], + (1, 2, 1): [-0.9999999999999998, 1.9999999999999998], + (1, 2, 2): [-0.33333333333333304, 1.333333333333333], + (1, 2, 3): [-0.14285714285714307, 1.142857142857143], + (1, 2, 4): [-0.06666666666666654, 1.0666666666666664], + (1, 2, 5): [-0.03225806451612906, 1.0322580645161292], + (1, 2, 6): [-0.015873015873015872, 1.0158730158730154], + (2, 1, 1): [0.33333333333333337, -2.0, 2.666666666666667], + (2, 1, 2): [0.04761904761904753, -0.5714285714285693, + 1.523809523809522], + (2, 1, 3): [0.009523809523810024, -0.2285714285714322, + 1.2190476190476225], + (2, 1, 4): [0.002150537634408055, -0.10322580645160284, + 1.1010752688171945], + (2, 1, 5): [0.0005120327700975248, -0.04915514592935677, + 1.0486431131592595], + (2, 1, 6): [0.0001249843769525012, -0.02399700037493191, + 1.0238720159979793], + (2, 2, 1): [0.1428571428571428, -1.428571428571427, + 2.2857142857142843], + (2, 2, 2): [0.022222222222222185, -0.444444444444444, + 1.4222222222222216], + (2, 2, 3): [0.004608294930875861, -0.1843317972350207, + 1.179723502304145], + (2, 2, 4): [0.0010582010582006751, -0.08465608465608221, + 1.0835978835978812], + (2, 2, 5): [0.0002540005080009476, -0.040640081280166496, + 1.0403860807721657], + (2, 2, 6): [6.224712107032182e-05, -0.01991907874258203, + 1.0198568316215115], + (3, 1, 1): [-0.04761904761904734, 0.6666666666666641, + -2.6666666666666594, 3.047619047619042], + (3, 1, 2): [-0.003174603174603108, 0.08888888888889318, + -0.7111111111111337, 1.6253968253968432], + (3, 1, 3): [-0.0003072196620577672, 0.01720430107525861, + -0.27526881720422713, 1.258371735791026], + (3, 1, 4): [-3.4135518007183396e-05, 0.003823178016754525, + -0.12234169653539884, 1.1185526540366513], + (3, 1, 5): [-4.031754094968587e-06, 0.0009031129172963892, + -0.0577992267083981, 1.056900145545197], + (3, 1, 6): [-4.901348115149418e-07, 0.00021958039560535103, + -0.02810629063481751, 1.0278872003740238], + (3, 2, 1): [-0.004608294930874168, 0.1935483870967698, + -1.5483870967741966, 2.359447004608302], + (3, 2, 2): [-0.00035273368606647537, 0.02962962962962734, + -0.47407407407406155, 1.444797178130501], + (3, 2, 3): [-3.628578685754835e-05, 0.006096012192020994, + -0.19507239014474764, 1.1890126637395837], + (3, 2, 4): [-4.149808071229888e-06, 0.0013943355119737377, + -0.08923747276669802, 1.0878472870627958], + (3, 2, 5): [-4.970655732572382e-07, 0.0003340280653114369, + -0.042755592360228634, 1.0424220613604906], + (3, 2, 6): [-6.08476257157875e-08, 8.177920896951241e-05, + -0.02093547748207586, 1.0208537591207332]} + for num_terms in [1, 2, 3]: + for step in [1, 2]: + for order in range(1, 7): + r_extrap = nd.Richardson(step_ratio=2.0, step=step, + num_terms=num_terms, order=order) + rule = r_extrap._get_richardson_rule() + # print((num_terms, step, order), rule.tolist()) + assert_array_almost_equal(rule, + true_vals[(num_terms, step, + order)]) + # self.assert_(False) + + def test_central(self): + method = 'central' + true_vals = {(1, 1): [-0.33333333, 1.33333333], + (1, 2): [-0.33333333, 1.33333333], + (1, 3): [-0.33333333, 1.33333333], + (1, 4): [-0.06666667, 1.06666667], + (1, 5): [-0.06666667, 1.06666667], + (1, 6): [-0.01587302, 1.01587302], + (2, 1): [0.02222222, -0.44444444, 1.42222222], + (2, 2): [0.02222222, -0.44444444, 1.42222222], + (2, 3): [0.02222222, -0.44444444, 1.42222222], + (2, 4): [1.05820106e-03, -8.46560847e-02, 1.08359788e+00], + (2, 5): [1.05820106e-03, -8.46560847e-02, 1.08359788e+00], + (2, 6): [6.22471211e-05, -1.99190787e-02, 1.01985683e+00]} + + for num_terms in [1, 2]: + for order in range(1, 7): + d = nd.Derivative(np.exp, method=method, order=order) + d._set_richardson_rule(step_ratio=2.0, num_terms=num_terms) + rule = d._richardson_extrapolate._get_richardson_rule() + assert_array_almost_equal(rule, + true_vals[(num_terms, order)]) + + def test_complex(self): + truth = { + (1, 2, 8): [9.576480164718605e-07, -0.004167684167715291, + 1.004166726519699], + (4, 2, 2): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (1, 2, 4): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (4, 1, 8): [-0.0039215686274510775, 1.0039215686274505], + (2, 2, 4): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (4, 2, 8): [9.576480164718605e-07, -0.004167684167715291, + 1.004166726519699], + (3, 1, 8): [-0.0039215686274510775, 1.0039215686274505], + (4, 1, 2): [-0.06666666666666654, 1.0666666666666664], + (3, 1, 6): [-0.06666666666666654, 1.0666666666666664], + (1, 1, 8): [-0.0039215686274510775, 1.0039215686274505], + (2, 1, 8): [-0.0039215686274510775, 1.0039215686274505], + (4, 1, 4): [-0.06666666666666654, 1.0666666666666664], + (3, 1, 4): [-0.06666666666666654, 1.0666666666666664], + (2, 1, 4): [-0.06666666666666654, 1.0666666666666664], + (3, 2, 2): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (2, 2, 8): [9.576480164718605e-07, -0.004167684167715291, + 1.004166726519699], + (2, 1, 6): [-0.06666666666666654, 1.0666666666666664], + (3, 1, 2): [-0.06666666666666654, 1.0666666666666664], + (4, 1, 6): [-0.06666666666666654, 1.0666666666666664], + (1, 1, 6): [-0.06666666666666654, 1.0666666666666664], + (1, 2, 2): [0.022222222222222185, -0.444444444444444, + 1.4222222222222216], + (3, 2, 6): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (1, 1, 4): [-0.06666666666666654, 1.0666666666666664], + (2, 1, 2): [-0.06666666666666654, 1.0666666666666664], + (4, 2, 4): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (3, 2, 4): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (2, 2, 2): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (1, 2, 6): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (4, 2, 6): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616], + (1, 1, 2): [-0.33333333333333304, 1.333333333333333], + (3, 2, 8): [9.576480164718605e-07, -0.004167684167715291, + 1.004166726519699], + (2, 2, 6): [0.0002614379084968331, -0.07111111111111235, + 1.070849673202616]} + # t = dict() + for n in [1, 2, 3, 4]: + for num_terms in [1, 2]: + for order in range(2, 9, 2): + d = nd.Derivative(np.exp, n=n, method='complex', + order=order) + d._set_richardson_rule(step_ratio=2.0, num_terms=num_terms) + rule = d._richardson_extrapolate._get_richardson_rule() + # t[(n, num_terms, order)] = rule.tolist() + assert_array_almost_equal(rule, + truth[(n, num_terms, order)]) + # print(t) + # self.assert_(False) + + def test_forward_backward(self): + truth = {(1, 1): [-1., 2.], + (1, 2): [-0.33333333, 1.33333333], + (1, 3): [-0.14285714, 1.14285714], + (1, 4): [-0.06666667, 1.06666667], + (1, 5): [-0.03225806, 1.03225806], + (1, 6): [-0.01587302, 1.01587302], + (2, 1): [0.33333333, -2., 2.66666667], + (2, 2): [0.04761905, -0.57142857, 1.52380952], + (2, 3): [0.00952381, -0.22857143, 1.21904762], + (2, 4): [0.00215054, -0.10322581, 1.10107527], + (2, 5): [5.12032770e-04, -4.91551459e-02, 1.04864311e+00], + (2, 6): [1.24984377e-04, -2.39970004e-02, 1.02387202e+00]} + for method in ['forward', 'backward']: + for num_terms in [1, 2]: + for order in range(1, 7): + d = nd.Derivative(np.exp, method=method, order=order) + d._set_richardson_rule(step_ratio=2.0, num_terms=num_terms) + rule = d._richardson_extrapolate._get_richardson_rule() + assert_array_almost_equal(rule, + truth[(num_terms, order)]) + + def _example_(self): + def f(x, h): + return (np.exp(x + h) - np.exp(x - h)) / (2.) + # f = lambda x, h: (np.exp(x+h)-np.exp(x)) + steps = [h for h in 2.0**-np.arange(10)] + df = [f(1, h) for h in steps] + print([dfi / hi for dfi, hi in zip(df, steps)]) + step = nd.MaxStepGenerator(step_ratio=2.0) + for method in ['central']: + d = nd.Derivative(np.exp, step=step, method=method) + for order in [2, 6]: + d.order = order + r_extrap = nd.Richardson(step_ratio=2.0, method=method, + num_terms=2, order=order) + + fd_rule = d._get_finite_difference_rule(step_ratio=2.0) + print(fd_rule) + df1, stepsi, _shape = d._apply_fd_rule(fd_rule, df, steps) + + rule = r_extrap._get_richardson_rule() + df2, error, hi = r_extrap(df1, stepsi) + + print(rule) + print(np.hstack((df2, error))) + + self.assert_(False) + + +class TestStepGenerator(unittest.TestCase): + + def test_default_generator(self): + step_gen = nd.MinStepGenerator(base_step=None, num_steps=10, + step_ratio=4, offset=-1) + h = np.array([h for h in step_gen(0)]) + desired = np.array([3.58968236e-02, 8.97420590e-03, 2.24355147e-03, + 5.60887869e-04, 1.40221967e-04, 3.50554918e-05, + 8.76387295e-06, 2.19096824e-06, 5.47742059e-07, + 1.36935515e-07]) + + assert_array_almost_equal((h - desired) / desired, 0) + + def test_default_base_step(self): + step_gen = nd.MinStepGenerator(num_steps=1, offset=0) + h = [h for h in step_gen(0)] + desired = nd.EPS ** (1. / 2.5) + assert_array_almost_equal((h[0] - desired) / desired, 0) + + def test_fixed_base_step(self): + desired = 0.1 + step_gen = nd.MinStepGenerator(base_step=desired, num_steps=1, scale=2, + offset=0) + h = [h for h in step_gen(0)] + assert_array_almost_equal((h[0] - desired) / desired, 0) + + +class TestFornbergWeights(unittest.TestCase): + + def test_weights(self): + x = np.r_[-1, 0, 1] + xbar = 0 + k = 1 + weights = nd.fornberg_weights(x, xbar, k) + np.testing.assert_allclose(weights, [-.5, 0, .5]) + + class TestDerivative(unittest.TestCase): + + def test_infinite_functions(self): + def finf(x): + return np.inf + df = nd.Derivative(finf) + val = df(0) + self.assert_(np.isnan(val)) + + def _example_fd_mat(self): + fdmat = nd.Derivative._fd_matrix(step_ratio=2.0, parity=1, nterms=3) + _fd_rules = np.linalg.pinv(fdmat) + self.assert_(False) + + def test_high_order_derivative_cos(self): + true_vals = (-1.0, 0.0, 1.0, 0.0, -1.0, 0.0) + methods = ['complex', 'multicomplex', 'central', + 'forward', 'backward'] + for method in methods: + n_max = dict(multicomplex=2, central=6).get(method, 5) + for n in range(1, n_max + 1): + true_val = true_vals[n - 1] + for order in range(2, 9, 2): + d3cos = nd.Derivative(np.cos, n=n, order=order, + method=method, full_output=True) + y, info = d3cos(np.pi / 2.0) + error = np.abs(y - true_val) + small = error <= info.error_estimate + if not small: + small = error < 10**(-12 + n) + if not small: + print('method=%s, n=%d, order=%d' % (method, n, order)) + print error, info.error_estimate + self.assertTrue(small) + # self.assert_(False) + + def test_derivative_of_cos_x(self): + x = np.r_[0, np.pi / 6.0, np.pi / 2.0] + true_vals = (-np.sin(x), -np.cos(x), np.sin(x), np.cos(x), -np.sin(x), + -np.cos(x)) + for method in ['complex', 'central', 'forward', 'backward']: + n_max = dict(complex=2, central=6).get(method, 5) + for n in range(1, n_max + 1): + true_val = true_vals[n - 1] + start, stop, step = dict(central=(2, 7, 2), + complex=(2, 3, 1)).get(method, + (1, 5, 1)) + for order in range(start, stop, step): + d3cos = nd.Derivative(np.cos, n=n, order=order, + method=method, full_output=True) + y, info = d3cos(x) + error = np.abs(y - true_val) + small = error <= info.error_estimate + if not small.all(): + small = np.where(small, small, error <= 10**(-11 + n)) + self.assertTrue(small.all()) + + def test_default_scale(self): + for method, scale in zip(['complex', 'central', 'forward', 'backward', + 'multicomplex'], + [1.35, 2.5, 2.5, 2.5, 1.35]): + np.testing.assert_allclose(scale, nd.default_scale(method, n=1)) + def test_derivative_cube(self): '''Test for Issue 7''' def cube(x): return x * x * x dcube = nd.Derivative(cube) shape = (3, 2) - x = np.ones(shape)*2 + x = np.ones(shape) * 2 dx = dcube(x) assert_array_almost_equal(list(dx.shape), list(shape), - decimal=12, + decimal=8, err_msg='Shape mismatch') txt = 'First differing element %d\n value = %g,\n true value = %g' - for i, (val, tval) in enumerate(zip(dx.ravel(), (3*x**2).ravel())): - assert_array_almost_equal(val, tval, decimal=12, + for i, (val, tval) in enumerate(zip(dx.ravel(), (3 * x**2).ravel())): + assert_array_almost_equal(val, tval, decimal=8, err_msg=txt % (i, val, tval)) def test_derivative_exp(self): # derivative of exp(x), at x == 0 dexp = nd.Derivative(np.exp) assert_array_almost_equal(dexp(0), np.exp(0), decimal=8) - dexp.n = 2 - t = dexp(0) - assert_array_almost_equal(t, np.exp(0)) def test_derivative_sin(self): # Evaluate the indicated (default = first) @@ -38,78 +360,25 @@ def test_derivative_sin(self): dsin = nd.Derivative(np.sin) x = np.linspace(0, 2. * np.pi, 13) y = dsin(x) - small = np.abs(y - np.cos(x)) <= dsin.error_estimate * 100 - self.assertTrue(np.all(small)) - - def test_second_and_fourth_derivative_of_sin(self): - # Higher order derivatives (second derivative) - # Truth: 0 - d2sin = nd.Derivative(np.sin, n=2, step_max=0.5) - assert_array_almost_equal(d2sin(np.pi), 0.0, decimal=8) - - # Higher order derivatives (up to the fourth derivative) - # Truth: sqrt(2)/2 = 0.707106781186548 - d2sin.n = 4 - y = d2sin(np.pi / 4) - small = np.abs(y - np.sqrt(2.) / 2.) < d2sin.error_estimate - self.assertTrue(small) - - def test_high_order_derivative_cos(self): - for method, orders in zip(['central', 'forward'], - [[2, 4], [1, 2, 3, 4]]): - for n, true_val in zip([1, 2, 3, 4, 5, 6], - (-1.0, 0.0, 1.0, 0.0, -1.0, 0.0)): - for order in orders: - d3cos = nd.Derivative(np.cos, n=n, order=order, - method=method) - y = d3cos(np.pi / 2.0) - small = np.abs(y - true_val) < 10 * d3cos.error_estimate - self.assertTrue(small) + np.testing.assert_almost_equal(y, np.cos(x), decimal=8) def test_backward_derivative_on_sinh(self): # Compute the derivative of a function using a backward difference # scheme. A backward scheme will only look below x0. dsinh = nd.Derivative(np.sinh, method='backward') - small = np.abs(dsinh(0.0) - np.cosh(0.0)) < dsinh.error_estimate - self.assertTrue(small) + self.assertAlmostEqual(dsinh(0.0), np.cosh(0.0)) def test_central_and_forward_derivative_on_log(self): # Although a central rule may put some samples in the wrong places, it # may still succeed - dlog = nd.Derivative(np.log, method='central') - x = 0.001 - small = np.abs(dlog(x) - 1.0 / x) < dlog.error_estimate - self.assertTrue(small) + epsilon = nd.MinStepGenerator(num_steps=15, offset=0, step_ratio=2) + dlog = nd.Derivative(np.log, method='central', step=epsilon) + x = 0.01 + self.assertAlmostEqual(dlog(x), 1.0 / x) # But forcing the use of a one-sided rule may be smart anyway - dlog.method = 'forward' - small = np.abs(dlog(x) - 1 / x) < dlog.error_estimate - self.assertTrue(small) - - def test_forward_derivative_on_tan(self): - # Control the behavior of Derivative - forward 2nd order method, with - # only 1 Romberg term. - dtan = nd.Derivative(np.tan, n=1, order=2, method='central', - romberg_terms=1, step_num=2+1+3, step_max=2**6*1e-15**(1./3)) - y = dtan(np.pi) - abserr = dtan.error_estimate - self.assertTrue(np.abs(y - 1.0) < abserr) - assert_array_almost_equal(y, 1.0, decimal=8) - - def test_derivative_poly1d(self): - p0 = np.poly1d(range(1, 6)) - fd = nd.Derivative(p0, n=4, romberg_terms=0) - p4 = p0.deriv(4) - assert_array_almost_equal(fd(1), p4(1), decimal=4) - - def test_vectorized_derivative_of_x2(self): - # Functions should be vectorized for speed, but its not - # always easy to do. - def fun(x): - return x**2 - df = nd.Derivative(fun, vectorized=True) - x = np.linspace(0, 5, 6) - assert_array_almost_equal(df(x), 2*x) + dlog = nd.Derivative(np.log, method='forward', step=epsilon) + self.assertAlmostEqual(dlog(x), 1 / x) class TestJacobian(unittest.TestCase): @@ -120,72 +389,97 @@ def testjacobian(self): def fun(c): return (c[0] + c[1] * np.exp(c[2] * xdata) - ydata) ** 2 - Jfun = nd.Jacobian(fun) - J = Jfun([1, 2, 0.75]) # should be numerically zero - for ji in J.ravel(): - assert_array_almost_equal(ji, 0.0) + + for method in ['complex', 'central', 'forward', 'backward']: + for order in [2, 4]: + Jfun = nd.Jacobian(fun, method=method, order=order) + J = Jfun([1, 2, 0.75]) # should be numerically zero + assert_array_almost_equal(J, np.zeros(J.shape)) class TestGradient(unittest.TestCase): + def testgradient(self): def fun(x): return np.sum(x ** 2) - dfun = nd.Gradient(fun) - d = dfun([1, 2, 3]) dtrue = [2., 4., 6.] - for (di, dit) in zip(d, dtrue): - assert_array_almost_equal(di, dit) - - -class TestHessian(unittest.TestCase): - - def testhessian(self): - # cos(x-y), at (0,0) - cos = np.cos - def fun(xy): - return cos(xy[0] - xy[1]) - Hfun2 = nd.Hessian(fun) - h2 = Hfun2([0, 0]) # h2 = [-1 1; 1 -1]; - htrue = [-1., 1., 1., -1.] - for (hi, hit) in zip(h2.ravel(), htrue): - assert_array_almost_equal(hi, hit) + for method in ['complex', 'central', 'backward', 'forward']: + for order in [2, 4]: + dfun = nd.Gradient(fun, method=method, order=order) + d = dfun([1, 2, 3]) + assert_array_almost_equal(d, dtrue) + # self.assert_(False) class TestHessdiag(unittest.TestCase): - def testhessdiag(self): + def test_complex(self): + def fun(x): + return x[0] + x[1] ** 2 + x[2] ** 3 + htrue = np.array([0., 2., 18.]) + method = 'complex' + for num_steps in range(3, 7, 1): + steps = nd.MinStepGenerator(num_steps=num_steps, + use_exact_steps=True, + step_ratio=2.0, offset=4) + Hfun = nd.Hessdiag(fun, step=steps, method=method, + full_output=True) + hd, _info = Hfun([1, 2, 3]) + _error = hd - htrue + assert_array_almost_equal(hd, htrue) + + def test_fixed_step(self): def fun(x): return x[0] + x[1] ** 2 + x[2] ** 3 - Hfun = nd.Hessdiag(fun) - hd = Hfun([1, 2, 3]) - htrue = [0., 2., 18.] - for (hi, hit) in zip(hd, htrue): - assert_array_almost_equal(hi, hit) + htrue = np.array([0., 2., 18.]) + + methods = ['multicomplex', 'complex', 'central', 'forward', 'backward'] + for order in range(2, 7, 2): + steps = nd.MinStepGenerator(num_steps=order + 1, + use_exact_steps=True, + step_ratio=3., offset=0) + for method in methods: + Hfun = nd.Hessdiag(fun, step=steps, method=method, order=order, + full_output=True) + hd, _info = Hfun([1, 2, 3]) + _error = hd - htrue + assert_array_almost_equal(hd, htrue) + + def test_default_step(self): + def fun(x): + return x[0] + x[1] ** 2 + x[2] ** 3 + htrue = np.array([0., 2., 18.]) + methods = ['central2', 'central', 'multicomplex', 'complex', 'forward', + 'backward'] + for order in range(2, 7, 2): + for method in methods: + Hfun = nd.Hessdiag(fun, method=method, order=order, + full_output=True) + hd, _info = Hfun([1, 2, 3]) + _error = hd - htrue + assert_array_almost_equal(hd, htrue) -class TestGlobalFunctions(unittest.TestCase): - def test_vec2mat(self): - mat = nd.core.vec2mat(np.arange(6), n=2, m=3) - assert_array_almost_equal(mat.tolist(), [[0, 1, 2], [1, 2, 3]], - decimal=12) +class TestHessian(unittest.TestCase): - mat = nd.core.vec2mat(np.arange(12), 3, 4) - assert_array_almost_equal(mat.tolist(), [[0, 1, 2, 3], - [1, 2, 3, 4], - [2, 3, 4, 5]], - decimal=12) + def test_hessian_cosIx_yI_at_I0_0I(self): + # cos(x-y), at (0,0) + + def fun(xy): + return np.cos(xy[0] - xy[1]) + htrue = [[-1., 1.], [1., -1.]] + methods = ['multicomplex', 'complex', 'central', 'central2', 'forward', + 'backward'] + for num_steps in [10, 1]: + step = nd.MinStepGenerator(num_steps=num_steps) + for method in methods: + Hfun2 = nd.Hessian(fun, method=method, step=step, + full_output=True) + h2, _info = Hfun2([0, 0]) + # print(method, (h2-np.array(htrue))) + assert_array_almost_equal(h2, htrue) - def testdea3(self): - def linfun(k): - return np.linspace(0, np.pi / 2., 2. ** (k + 5) + 1) - Ei = np.zeros(3) - for k in np.arange(3): - x = linfun(k) - Ei[k] = np.trapz(np.sin(x), x) - [En, err] = nd.dea3(Ei[0], Ei[1], Ei[2]) - self.assertTrue(np.abs(En - 1) < err) - assert_array_almost_equal(En, 1.0, decimal=8) if __name__ == '__main__': unittest.main() diff --git a/setup.py b/setup.py index a090b5c..f2e144d 100644 --- a/setup.py +++ b/setup.py @@ -9,17 +9,26 @@ Usage: Run all tests: -python setup.py test + python setup.py test -python setup.py install [, --prefix=$PREFIX] + python setup.py doctests -python setup.py bdist_wininst +Build documentation + + python setup.py docs + +Install + python setup.py install [, --prefix=$PREFIX] + +Build wininstaller + + python setup.py bdist_wininst PyPi upload: -python setup.py sdist bdist_wininst upload --show-response + python setup.py sdist bdist_wininst upload --show-response -python setup.py register sdist bdist_wininst upload --show-response + python setup.py register sdist bdist_wininst upload --show-response """ import os