pbrod/numdifftools

View on GitHub
src/numdifftools/finite_difference.py

Summary

Maintainability
D
2 days
Test Coverage
"""
Finite difference methods module.
"""
from __future__ import absolute_import, division, print_function
import warnings
import numpy as np
from numpy import linalg
from scipy import special
from numdifftools.extrapolation import convolve
from numdifftools.multicomplex import Bicomplex

_SQRT_J = (1j + 1.0) / np.sqrt(2.0)  # = 1j**0.5

# step_ratio, parity, nterms
FD_RULES = {}
# (2.0, 1, 1): array([[1.]]),
# (2.0, 1, 2): array([[-0.333333333333333333333, 2.666666666666666666666666666],
#                     [8., -16.]]),
# (2.0, 1, 3): array([[2.22222222222222222e-02, -8.8888888888889e-01, 5.6888888888888889e+00],
#                     [-2.666666666666667e+00, 9.0666666666666667e+01, -1.7066666666666667e+02],
#                     [1.7066666666666667e+02, -1.7066666666666667e+03, 2.7306666666666667e+03]]),
# (2.0, 0, 2): array([[-1., 4.],
#                     [4., -8.]]),
# (2.0, 0, 4): array([[-4.76190476e-02, 1.33333333e+00, -1.06666667e+01, 2.43809524e+01],
#                     [1.33333333e+00, -3.46666667e+01, 2.34666667e+02, -3.41333333e+02],
#                     [-1.60000000e+01, 3.52000000e+02, -1.66400000e+03, 2.04800000e+03],
#                     [7.31428571e+01, -1.02400000e+03, 4.09600000e+03, -4.68114286e+03]])}

# @PydevIgnore


def _assert(cond, msg):
    if not cond:
        raise ValueError(msg)


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.0 and then subtracting 1.0.
    """
    return (h + 1.0) - 1.0


class DifferenceFunctions(object):
    """
    Class defining difference functions

    Notes
    -----
    The d
    """
    # pylint: disable=unused-argument
    @staticmethod
    def _central_even(f, f_x0i, x0i, h):
        return (f(x0i + h) + f(x0i - h)) / 2.0 - f_x0i

    @staticmethod
    def _central(f, f_x0i, x0i, h):  # @UnusedVariable
        return (f(x0i + h) - f(x0i - h)) / 2.0

    @staticmethod
    def _forward(f, f_x0i, x0i, h):
        return f(x0i + h) - f_x0i

    @staticmethod
    def _backward(f, f_x0i, x0i, h):
        return f_x0i - f(x0i - h)

    @staticmethod
    def _complex(f, f_x, x, h):  # @UnusedVariable
        return f(x + 1j * h).imag

    @staticmethod
    def _complex_odd(f, f_x, x, h):
        i_h = h * _SQRT_J
        return ((_SQRT_J / 2.) * (f(x + i_h) - f(x - i_h))).imag

    @staticmethod
    def _complex_odd_higher(f, f_x, x, h):
        i_h = h * _SQRT_J
        return ((3 * _SQRT_J) * (f(x + i_h) - f(x - i_h))).real

    @staticmethod
    def _complex_even(f, f_x, x, h):
        i_h = h * _SQRT_J
        return (f(x + i_h) + f(x - i_h)).imag

    @staticmethod
    def _complex_even_higher(f, f_x, x, h):
        i_h = h * _SQRT_J
        return 12.0 * (f(x + i_h) + f(x - i_h) - 2 * f_x).real

    @staticmethod
    def _multicomplex(f, f_x, x, h):
        z = Bicomplex(x + 1j * h, 0)
        return Bicomplex.__array_wrap__(f(z)).imag

    @staticmethod
    def _multicomplex2(f, f_x, x, h):
        z = Bicomplex(x + 1j * h, h)
        return Bicomplex.__array_wrap__(f(z)).imag12


class JacobianDifferenceFunctions(object):
    """Class defining Jacobian difference functions"""
    # pylint: disable=unused-argument
    @staticmethod
    def increments(n, h):
        """Returns Jacobian steps"""
        e_i = np.zeros(np.shape(h), float)
        for k in range(n):
            e_i[k] = h[k]
            yield e_i
            e_i[k] = 0

    @staticmethod
    def _central(f, f_x, x, h):
        n = len(x)
        steps = JacobianDifferenceFunctions.increments(n, h)
        return np.array([(f(x + hi) - f(x - hi)) / 2.0 for hi in steps])

    @staticmethod
    def _central_even(f, f_x, x, h):
        n = len(x)
        steps = JacobianDifferenceFunctions.increments(n, h)
        return np.array([(f(x + hi) + f(x - hi)) / 2.0 - f_x for hi in steps])

    @staticmethod
    def _backward(f, f_x, x, h):
        n = len(x)
        steps = JacobianDifferenceFunctions.increments(n, h)
        return np.array([f_x - f(x - hi) for hi in steps])

    @staticmethod
    def _forward(f, f_x, x, h):
        n = len(x)
        steps = JacobianDifferenceFunctions.increments(n, h)
        return np.array([f(x + hi) - f_x for hi in steps])

    @staticmethod
    def _complex(f, f_x, x, h):
        n = len(x)
        steps = JacobianDifferenceFunctions.increments(n, h)
        return np.array([f(x + 1j * ih).imag for ih in steps])

    @staticmethod
    def _complex_even(f, f_x, x, h):
        n = len(x)
        j_1 = _SQRT_J
        steps = JacobianDifferenceFunctions.increments(n, h)
        return np.array([(f(x + j_1*ih) + f(x - j_1*ih)).imag for ih in steps])

    @staticmethod
    def _complex_odd(f, f_x, x, h):
        n = len(x)
        j_1 = _SQRT_J
        steps = JacobianDifferenceFunctions.increments(n, h)
        return np.array([((j_1 / 2.) * (f(x + j_1 * ih) - f(x - j_1 * ih))).imag for ih in steps])

    @staticmethod
    def _multicomplex(f, f_x, x, h):
        n = len(x)
        steps = JacobianDifferenceFunctions.increments(n, h)
        cmplx_wrap = Bicomplex.__array_wrap__
        partials = [cmplx_wrap(f(Bicomplex(x + 1j * hi, 0))).imag for hi in steps]
        return np.array(partials)


class HessdiagDifferenceFunctions(object):
    """Class defining Hessdiag difference functions

    References
    ----------
    Ridout, M.S. (2009) Statistical applications of the complex-step method
        of numerical differentiation. The American Statistician, 63, 66-74
    """
    # pylint: disable=unused-argument
    @staticmethod
    def _central2(f, f_x, x, h):
        """Eq. 8 in Ridout (2009)."""
        n = len(x)
        increments = np.identity(n) * h
        partials = [(f(x + 2 * hi) + f(x - 2 * hi)
                     + 2 * f_x - 2 * f(x + hi) - 2 * f(x - hi)) / 4.0
                    for hi in increments]
        return np.array(partials)

    @staticmethod
    def _central_even(f, f_x, x, h):
        """Eq. 9 in Ridout (2009)."""
        n = len(x)
        increments = np.identity(n) * h
        partials = [(f(x + hi) + f(x - hi)) / 2.0 - f_x for hi in increments]
        return np.array(partials)

    @staticmethod
    def _backward(f, f_x, x, h):
        n = len(x)
        increments = np.identity(n) * h
        partials = [f_x - f(x - hi) for hi in increments]
        return np.array(partials)

    @staticmethod
    def _forward(f, f_x, x, h):
        n = len(x)
        increments = np.identity(n) * h
        partials = [f(x + hi) - f_x for hi in increments]
        return np.array(partials)

    @staticmethod
    def _multicomplex2(f, f_x, x, h):
        n = len(x)
        increments = np.identity(n) * h
        cmplx_wrap = Bicomplex.__array_wrap__
        partials = [cmplx_wrap(f(Bicomplex(x + 1j * hi, hi))).imag12
                    for hi in increments]
        return np.array(partials)

    @staticmethod
    def _complex_even(f, f_x, x, h):
        n = len(x)
        increments = np.identity(n) * h * (1j + 1) / np.sqrt(2)
        partials = [(f(x + hi) + f(x - hi)).imag for hi in increments]
        return np.array(partials)


class HessianDifferenceFunctions(object):
    """Class defining Hessian difference functions

    References
    ----------
    Ridout, M.S. (2009)
    "Statistical applications of the complex-step method of numerical differentiation",
    The American Statistician, 63, 66-74
    """
    # pylint: disable=unused-argument
    @staticmethod
    def _complex_even(f, f_x, x, h):
        """
        Calculate Hessian with complex-step derivative approximation

        The stepsize is the same for the complex and the finite difference part

        Eq 10 in Ridout (2009).
        """
        n = len(x)
        eee = np.diag(h)
        hess = 2. * np.outer(h, h)
        for i in range(n):
            for j in range(i, n):
                hess[i, j] = (f(x + 1j * eee[i] + eee[j])
                              - f(x + 1j * eee[i] - eee[j])).imag / hess[j, i]
                hess[j, i] = hess[i, j]
        return hess

    @staticmethod
    def _multicomplex2(f, f_x, x, h):
        """Calculate Hessian with Bicomplex-step derivative approximation"""
        n = len(x)
        eee = np.diag(h)
        hess = np.outer(h, h)
        cmplx_wrap = Bicomplex.__array_wrap__
        for i in range(n):
            for j in range(i, n):
                zph = Bicomplex(x + 1j * eee[i, :], eee[j, :])
                hess[i, j] = cmplx_wrap(f(zph)).imag12 / hess[j, i]
                hess[j, i] = hess[i, j]
        return hess

    @staticmethod
    def _central_even(f, f_x, x, h):
        """Eq 9 in Ridout (2009)."""
        n = len(x)
        eee = np.diag(h)
        dtype = np.result_type(f_x, float)  # make sure it is at least float64
        hess = np.empty((n, n), dtype=dtype)
        np.outer(h, h, out=hess)
        for i in range(n):
            e_i = eee[i, :]
            hess[i, i] = (f(x + 2 * e_i) - 2 * f_x + f(x - 2 * e_i)) / (4. * hess[i, i])
            for j in range(i + 1, n):
                e_j = eee[j, :]
                hess[i, j] = (f(x + e_i + e_j) - f(x + e_i - e_j)
                              - f(x - e_i + e_j) + f(x - e_i - e_j)) / (4. * hess[j, i])
                hess[j, i] = hess[i, j]
        return hess

    @staticmethod
    def _central2(f, f_x, x, h):
        """Eq. 8 in Ridout (2009)"""
        n = len(x)
        eee = np.diag(h)
        dtype = np.result_type(f_x, float)
        f_xpe = np.empty(n, dtype=dtype)
        f_xme = np.empty(n, dtype=dtype)
        for i in range(n):
            f_xpe[i] = f(x + eee[i])
            f_xme[i] = f(x - eee[i])

        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 + eee[i, :] + eee[j, :])
                              + f(x - eee[i, :] - eee[j, :])
                              - f_xpe[i] - f_xpe[j] + f_x
                              - f_xme[i] - f_xme[j] + f_x) / (2 * hess[j, i])
                hess[j, i] = hess[i, j]
        return hess

    @staticmethod
    def _forward(f, f_x, x, h):
        """Eq. 7 in Ridout (2009)"""
        n = len(x)
        eee = np.diag(h)
        dtype = np.result_type(f_x, float)
        g = np.empty(n, dtype=dtype)
        for i in range(n):
            g[i] = f(x + eee[i, :])

        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 + eee[i, :] + eee[j, :]) - g[i] - g[j] + f_x) / hess[j, i]
                hess[j, i] = hess[i, j]
        return hess

    @staticmethod
    def _backward(f, f_x, x, h):
        return HessianDifferenceFunctions._forward(f, f_x, x, -h)


class LogRule(object):
    """Log spaced finite difference rule class

    Parameters
    ----------
    n : int, optional
        Order of the derivative.
    method : {'central', 'complex', 'multicomplex', 'forward', 'backward'}
        defines the method used in the approximation
    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.


    Examples
    --------
    >>> from numdifftools.finite_difference import LogRule
    >>> np.allclose(LogRule(n=1, method='central', order=2).rule(step_ratio=2.0), 1)
    True
    >>> np.allclose(LogRule(n=1, method='central', order=4).rule(step_ratio=2.),
    ...             [-0.33333333,  2.66666667])
    True
    >>> np.allclose(LogRule(n=1, method='central', order=6).rule(step_ratio=2.),
    ...             [ 0.02222222, -0.88888889,  5.68888889])
    True

    >>> np.allclose(LogRule(n=1, method='forward', order=2).rule(step_ratio=2.), [-1.,  4.])
    True

    >>> np.allclose(LogRule(n=1, method='forward', order=4).rule(step_ratio=2.),
    ...             [ -0.04761905,   1.33333333, -10.66666667,  24.38095238])
    True
    >>> np.allclose(LogRule(n=1, method='forward', order=6).rule(step_ratio=2.),
    ...    [ -1.02406554e-04,   1.26984127e-02,  -5.07936508e-01,
    ...       8.12698413e+00,  -5.20126984e+01,   1.07381055e+02])
    True
    >>> step_ratio=2.0
    >>> fd_rule = LogRule(n=2, method='forward', order=4)
    >>> h = 0.002*(1./step_ratio)**np.arange(6)
    >>> x0 = 1.
    >>> f = np.exp
    >>> f_x0 = f(x0)
    >>> f_del = f(x0+h) - f_x0  # forward difference
    >>> f_del = fd_rule.diff(f, f_x0, x0, h)  # or alternatively
    >>> fder, h, shape = fd_rule.apply(f_del, h, step_ratio)
    >>> np.allclose(fder, f(x0))
    True

    """
    _difference_functions = DifferenceFunctions()

    def __init__(self, n=1, method='central', order=2):
        self.n = n
        self.method = method
        self.order = order

# --- properties ---

    @property
    def _odd_derivative(self):
        return self.n % 2 == 1

    @property
    def _even_derivative(self):
        return self.n % 2 == 0

    @property
    def _derivative_mod_four_is_three(self):
        return self.n % 4 == 3

    @property
    def _derivative_mod_four_is_zero(self):
        return self.n % 4 == 0

    @property
    def eval_first_condition(self):
        """True if f(x0) needs to be evaluated given the differentiation method."""
        even_derivative = self._even_derivative
        return ((even_derivative and self.method in ('central', 'central2')) or
                self.method in ['forward', 'backward'] or
                self.method == 'complex' and self._derivative_mod_four_is_zero)

    @property
    def _complex_high_order(self):
        return self.method == 'complex' and (self.n > 1 or self.order >= 4)

    @property
    def richardson_step(self):
        """The step between exponents in the error polynomial of the Richardson extrapolation."""
        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)

    @property
    def method_order(self):
        """The leading order of the truncation error of the Richardson extrapolation."""
        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 _parity_complex(self, order, method_order):
        if self.n == 1 and method_order < 4:
            return (order % 2) + 1
        return (3
                + 2 * int(self._odd_derivative)
                + int(self._derivative_mod_four_is_three)
                + int(self._derivative_mod_four_is_zero))

    def _parity(self, method, order, method_order):
        if method.startswith('central'):
            return (order % 2) + 1
        if method == 'complex':
            return self._parity_complex(order, method_order)
        return 0

    @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)
            5 (only every 4'th order terms included starting from order 1)
            6 (only every 4'th order terms included starting from order 3)
        nterms : scalar, integer
            number of terms
        """
        _assert(0 <= parity <= 6,
                'Parity must be 0, 1, 2, 3, 4, 5 or 6! ({0:d})'.format(parity))
        step = [1, 2, 2, 4, 4, 4, 4][parity]
        inv_sr = 1.0 / step_ratio
        offset = [1, 1, 2, 2, 4, 1, 3][parity]
        c_0 = [1.0, 1.0, 1.0, 2.0, 24.0, 1.0, 6.0][parity]
        c = c_0 / special.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)))

    @property
    def _flip_fd_rule(self):
        return ((self._even_derivative and (self.method == 'backward'))
                or (self.method == 'complex' and (self.n % 8 in [3, 4, 5, 6])))

    @property
    def _multicomplex_middle_name(self):
        if self.method == 'multicomplex' and self.n > 1:
            _assert(self.n <= 2, 'Multicomplex method only support first '
                    'and second order derivatives.')
            return '2'
        return ''

    def _get_middle_name(self):
        if self._even_derivative and self.method in ('central', 'complex'):
            return '_even'
        if self._complex_high_order and self._odd_derivative:
            return '_odd'
        return self._multicomplex_middle_name

    def _get_last_name(self):
        last = ''
        if (self.method == 'complex' and self._derivative_mod_four_is_zero or
                self._complex_high_order and
                self._derivative_mod_four_is_three):
            last = '_higher'
        return last

    @property
    def diff(self):
        "The difference function"
        first = '_{0!s}'.format(self.method)
        middle = self._get_middle_name()
        last = self._get_last_name()
        name = first + middle + last
        return getattr(self._difference_functions, name)

    def rule(self, step_ratio=2.0):
        """
        Return finite differencing rule.

        Parameters
        ----------
        step_ratio : real scalar, optional, default 2.0
            Ratio between sequential steps generated.

        Notes
        -----
        The rule is for a nominal unit step size, and must be scaled later
        to reflect the local step size.

        Member method used:  _fd_matrix

        Member variables used:
        n
        order
        method
        """
        method = self.method
        if method in ('multicomplex',) or self.n == 0:
            return np.ones((1,))
        step_ratio = make_exact(step_ratio)
        order, method_order = self.n - 1, self.method_order
        parity = self._parity(method, order, method_order)
        step = self.richardson_step
        num_terms = (order + method_order) // step
        fd_rules = FD_RULES.get((step_ratio, parity, num_terms))
        if fd_rules is None:
            fd_mat = self._fd_matrix(step_ratio, parity, num_terms)
            fd_rules = linalg.pinv(fd_mat)
            FD_RULES[(step_ratio, parity, num_terms)] = fd_rules

        rule_index = order // step
        if self._flip_fd_rule:
            return -fd_rules[rule_index]
        return fd_rules[rule_index]

    @staticmethod
    def _vstack(sequence, steps):
        original_shape = np.shape(sequence[0])
        f_del = np.vstack([np.ravel(r) for r in sequence])
        one = np.ones(original_shape)
        h = np.vstack([np.ravel(one * step) for step in steps])
        _assert(f_del.size == h.size, 'fun did not return data of correct '
                'size (it must be vectorized)')
        return f_del, h, original_shape

    def apply(self, sequence, steps, step_ratio=2.0):
        """
        Apply finite difference rule along the first axis.

        Return derivative estimates of fun at x0 for a sequence of stepsizes h

        Parameters
        ----------
        sequence: finite differences
        steps: steps

        """
        f_del, h, original_shape = self._vstack(sequence, steps)
        der_init, h = self._apply(f_del, h, step_ratio)
        return der_init, h, original_shape

    def _apply(self, f_del, h, step_ratio=2.0):
        """
        Apply finite difference rule along the first axis.

        Return derivative estimates of fun at x0 for a sequence of stepsizes h

        Parameters
        ----------
        f_del: finite differences
        h: steps

        """
        fd_rule = self.rule(step_ratio)

        num_steps = h.shape[0]
        n_r = fd_rule.size - 1
        _assert(n_r < num_steps, 'num_steps ({0:d}) must  be larger than '
                '({1:d}) n + order - 1 = {2:d} + {3:d} -1'
                ' ({4:s})'.format(num_steps, n_r + 1, self.n, self.order, self.method))
        f_diff = convolve(f_del, fd_rule[::-1], axis=0, origin=n_r // 2)

        der_init = f_diff / (h ** self.n)
        num_steps = max(num_steps - n_r, 1)
        return der_init[:num_steps], h[:num_steps]


class LogJacobianRule(LogRule):
    """Log spaced finite difference Jacobian rule class

    Parameters
    ----------
    n : 1
        Order of the derivative.
    method : {'central', 'complex', 'multicomplex', 'forward', 'backward'}
        defines the method used in the approximation
    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.

    Examples
    --------
    >>> from numdifftools.finite_difference import LogJacobianRule as Rule
    >>> np.allclose(Rule(n=1, method='central', order=2).rule(step_ratio=2.0), 1)
    True
    >>> np.allclose(Rule(n=1, method='central', order=4).rule(step_ratio=2.),
    ...             [-0.33333333,  2.66666667])
    True
    >>> np.allclose(Rule(n=1, method='central', order=6).rule(step_ratio=2.),
    ...             [ 0.02222222, -0.88888889,  5.68888889])
    True

    >>> np.allclose(Rule(n=1, method='forward', order=2).rule(step_ratio=2.), [-1.,  4.])
    True

    >>> np.allclose(Rule(n=1, method='forward', order=4).rule(step_ratio=2.),
    ...             [ -0.04761905,   1.33333333, -10.66666667,  24.38095238])
    True
    >>> np.allclose(Rule(n=1, method='forward', order=6).rule(step_ratio=2.),
    ...    [ -1.02406554e-04,   1.26984127e-02,  -5.07936508e-01,
    ...       8.12698413e+00,  -5.20126984e+01,   1.07381055e+02])
    True
    >>> step_ratio=2.0
    >>> fd_rule = Rule(n=1, method='forward', order=4)
    >>> steps = 0.002*(1./step_ratio)**np.arange(6)

    >>> x0 = np.atleast_1d(1.)
    >>> f = np.exp
    >>> f_x0 = f(x0)
    >>> f_del = [f(x0+h) - f_x0  for h in steps] # forward difference
    >>> f_del = [fd_rule.diff(f, f_x0, x0, h) for h in steps[:, None]]  # or alternatively
    >>> fder, h, shape = fd_rule.apply(f_del, steps, step_ratio)
    >>> np.allclose(fder, f(x0))
    True

    """
    _difference_functions = JacobianDifferenceFunctions()

    # n = property(fget=lambda cls: 1, fset=lambda cls, n: None)

    @staticmethod
    def _atleast_2d(original_shape, ndim):
        if ndim == 1:
            original_shape = (1,) + tuple(original_shape)
        return tuple(original_shape)

    def _vstack(self, sequence, steps):
        original_shape = list(np.shape(np.atleast_1d(sequence[0])))
        ndim = len(original_shape)
        if sum(original_shape) == ndim:
            f_del = np.vstack(sequence)
            h = np.vstack(steps)
        else:
            axes = [0, 1, 2][:ndim]
            axes[:2] = axes[1::-1]
            original_shape[:2] = original_shape[1::-1]

            f_del = np.vstack([np.atleast_1d(r).transpose(axes).ravel() for r in sequence])
            h = np.vstack([np.atleast_1d(r).transpose(axes).ravel() for r in steps])
        _assert(f_del.size == h.size, 'fun did not return data of correct '
                'size (it must be vectorized)')
        return f_del, h, self._atleast_2d(original_shape, ndim)


class LogHessdiagRule(LogRule):
    """Log spaced finite difference Hessdiag rule class

    Parameters
    ----------
    n : 2
        Order of the derivative.
    method : {'central', 'complex', 'multicomplex', 'forward', 'backward'}
        defines the method used in the approximation
    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.

    Examples
    --------
    >>> import numpy as np
    >>> from numdifftools.finite_difference import LogHessdiagRule as Rule

    >>> np.allclose(Rule(method='central', order=2).rule(step_ratio=2.0), 2)
    True
    >>> np.allclose(Rule(method='central', order=4).rule(step_ratio=2.),
    ...             [-0.66666667, 10.66666667])
    True
    >>> np.allclose(Rule(method='central', order=6).rule(step_ratio=2.),
    ...             [ 4.44444444e-02, -3.55555556e+00,  4.55111111e+01])
    True
    >>> np.allclose(Rule(method='forward', order=2).rule(step_ratio=2.), [ -4.,  40., -64.])
    True
    >>> np.allclose(Rule(method='forward', order=4).rule(step_ratio=2.),
    ...          [-1.90476190e-01,  1.10476190e+01, -1.92000000e+02,
    ...          1.12152381e+03, -1.56038095e+03])
    True
    >>> np.allclose(Rule(method='forward', order=6).rule(step_ratio=2.),
    ...    [-4.09626216e-04,  1.02406554e-01, -8.33015873e+00,  2.76317460e+02,
    ...     -3.84893968e+03,  2.04024004e+04, -2.74895500e+04])
    True
    >>> step_ratio=2.0
    >>> fd_rule = Rule(method='forward', order=4)
    >>> steps = 0.002*(1./step_ratio)**np.arange(6)
    >>> x0 = np.array([0., 0.])
    >>> f = lambda xy : np.cos(xy[0]-xy[1])
    >>> f_x0 = f(x0)
    >>> f_del = [f(x0+h) - f_x0  for h in steps] # forward difference
    >>> f_del = [fd_rule.diff(f, f_x0, x0, h) for h in steps]  # or alternatively
    >>> fder, h, shape = fd_rule.apply(f_del, steps, step_ratio)

    >>> np.allclose(fder, [[-1., -1.], [-1., -1.]])
    True

    """
    _difference_functions = HessdiagDifferenceFunctions()

    n = property(fget=lambda cls: 2, fset=lambda cls, n: None)


class LogHessianRule(LogRule):
    """Log spaced finite difference Hessian rule class

    Parameters
    ----------
    n : 2
        Order of the derivative.
    method : {'central', 'complex', 'multicomplex', 'forward', 'backward'}
        defines the method used in the approximation
    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.
    """
    _difference_functions = HessianDifferenceFunctions()

    n = property(fget=lambda cls: 2, fset=lambda cls, unused_n: None)

    @property
    def order(self):
        """The order of the error term in the Taylor approximation used"""
        return dict(backward=1, forward=1).get(self.method, 2)

    @order.setter
    def order(self, order):
        valid_order = self.order

        if order is not None and order != valid_order:
            msg = 'Can not change order to {}! The only valid order is {} for method={}.'
            warnings.warn(msg.format(order, valid_order, self.method))

    @property
    def _complex_high_order(self):
        return False

    def apply(self, sequence, steps, step_ratio=2.0):
        """
        Apply finite difference rule along the first axis.

        Return derivative estimates of fun at x0 for a sequence of stepsizes h

        Parameters
        ----------
        sequence: finite differences
        steps: steps

        """
        f_del, h, original_shape = self._vstack(sequence, steps)
        return f_del, h, original_shape


if __name__ == '__main__':
    from numdifftools.testing import test_docstrings
    test_docstrings()