Skip to content

Commit

Permalink
Ongoing work to make nd_cstep the new standard tool
Browse files Browse the repository at this point in the history
  • Loading branch information
@pbrod
pbrod committed Aug 19, 2015
1 parent f55e7ed commit 2a88fb2
Show file tree
Hide file tree
Showing 3 changed files with 301 additions and 49 deletions.
322 changes: 280 additions & 42 deletions numdifftools/nd_cstep.py
View file
Original file line number Diff line number Diff line change
@@ -1,51 +1,24 @@
"""numerical differentiation functions:
Derivative, Gradient, Jacobian, and Hessian
Author : pbrod, josef-pkt
License : BSD
Notes
-----
These are simple forward differentiation, so that we have them available
without dependencies.
* Jacobian should be faster than numdifftools.core because it doesn't use loop
over observations.
* numerical precision will vary and depend on the choice of stepsizes
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 version 1.0 released 12/27/2006 by John D'Errico
(e-mail: [email protected])
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.
"""

# TODO:
# * some cleanup
# * check numerical accuracy (and bugs) with numdifftools and analytical
# derivatives
# - linear least squares case: (hess - 2*X'X) is 1e-8 or so
# - gradient and Hessian agree with numdifftools when evaluated away from
# minimum
# - forward gradient, Jacobian evaluated at minimum is inaccurate, centered
# (+/- base_step) is ok
# * dot product of Jacobian is different from Hessian, either wrong example or
# a bug (unlikely), or a real difference
#
#
# What are the conditions that Jacobian dotproduct and Hessian are the same?
#
# See also:
#
# BHHH: Greene p481 17.4.6, MLE Jacobian = d loglike / d beta , where loglike
# is vector for each observation
# see also example 17.4 when J'J is very different from Hessian
# also does it hold only at the minimum, what's relationship to covariance
# of Jacobian matrix
# http://projects.scipy.org/scipy/ticket/1157
# http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
# objective: sum((y-f(beta,x)**2), Jacobian = d f/d beta
# and not d objective/d beta as in MLE Greene similar:
# http://crsouza.blogspot.com/2009/11/neural-network-learning-by-levenberg_18.html#hessian
#
# in example: if J = d x*beta / d beta then J'J == X'X
# similar to
# http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
from __future__ import print_function
import numpy as np
from numdifftools.core import dea3
from collections import namedtuple
from matplotlib import pyplot as plt
from numdifftools.multicomplex import bicomplex
Expand All @@ -55,6 +28,8 @@
from scipy.ndimage.filters import convolve1d
import warnings
# 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

Expand All @@ -73,6 +48,268 @@
np.arange(-4, 5))}


class Dea(object):
'''
LIMEXP is the maximum number of elements the
epsilon table data can contain. The epsilon table
is stored in the first (LIMEXP+2) entries of EPSTAB.
LIST OF MAJOR VARIABLES
-----------------------
E0,E1,E2,E3 - DOUBLE PRECISION
The 4 elements on which the computation of
a new element in the epsilon table is based.
NRES - INTEGER
Number of extrapolation results actually
generated by the epsilon algorithm in prior
calls to the routine.
NEWELM - INTEGER
Number of elements to be computed in the
new diagonal of the epsilon table. The
condensed epsilon table is computed. Only
those elements needed for the computation of
the next diagonal are preserved.
RES - DOUBLE PREISION
New element in the new diagonal of the
epsilon table.
ERROR - DOUBLE PRECISION
An estimate of the absolute error of RES.
Routine decides whether RESULT=RES or
RESULT=SVALUE by comparing ERROR with
ABSERR from the previous call.
RES3LA - DOUBLE PREISION
Vector of DIMENSION 3 containing at most
the last 3 results.
'''
def __init__(self, limexp=3):
self.limexp = 2 * (limexp // 2) + 1
self.epstab = np.zeros(limexp+5)
self.ABSERR = 10.
self._n = 0
self._nres = 0
if (limexp < 3):
raise ValueError('LIMEXP IS LESS THAN 3')

def _compute_error(self, RES3LA, NRES, RES):
fact = [6.0, 2.0, 1.0][min(NRES-1, 2)]
error = fact * np.abs(RES - RES3LA[:NRES]).sum()
return error

def _shift_table(self, EPSTAB, N, NEWELM, NUM):
i_0 = 1 if ((NUM // 2) * 2 == NUM - 1) else 0
i_n = 2 * NEWELM + 2
EPSTAB[i_0:i_n:2] = EPSTAB[i_0 + 2:i_n + 2:2]

if (NUM != N):
i_n = NUM - N
EPSTAB[:N + 1] = EPSTAB[i_n:i_n + N + 1]
return EPSTAB

def _update_RES3LA(self, RES3LA, RESULT, NRES):
if NRES > 2:
RES3LA[:2] = RES3LA[1:]
RES3LA[2] = RESULT
else:
RES3LA[NRES] = RESULT

def __call__(self, SVALUE):

EPSTAB = self.epstab
RES3LA = EPSTAB[-3:]
RESULT = SVALUE
N = self._n
NRES = self._nres
EPSTAB[N] = SVALUE
if (N == 0):
ABSERR = abs(RESULT)
elif (N == 1):
ABSERR = 6.0 * abs(RESULT - EPSTAB[0])
else:
ABSERR = self.ABSERR
EPSTAB[N + 2] = EPSTAB[N]
NEWELM = N // 2
NUM = N
K1 = N
for I in range(NEWELM):
E0 = EPSTAB[K1 - 2]
E1 = EPSTAB[K1 - 1]
E2 = RES = EPSTAB[K1 + 2]
DELTA2, DELTA3 = E2 - E1, E1 - E0
ERR2, ERR3 = abs(DELTA2), abs(DELTA3)
TOL2 = max(abs(E2), abs(E1)) * _EPS
TOL3 = max(abs(E1), abs(E0)) * _EPS
converged = (ERR2 <= TOL2 and ERR3 <= TOL3)
if converged:
ABSERR = ERR2 + ERR3
RESULT = RES
break
if (I != 0):
E3 = EPSTAB[K1]
DELTA1 = E1 - E3
ERR1 = abs(DELTA1)
TOL1 = max(abs(E1), abs(E3)) * _EPS
converged = (ERR1 <= TOL1 or ERR2 <= TOL2 or
ERR3 <= TOL3)
if not converged:
SS = 1.0 / DELTA1 + 1.0 / DELTA2 - 1.0 / DELTA3
else:
converged = (ERR2 <= TOL2 or ERR3 <= TOL3)
if not converged:
SS = 1.0 / DELTA2 - 1.0 / DELTA3
EPSTAB[K1] = E1
if (converged or abs(SS * E1) <= 1e-04):
N = 2 * I
if (NRES == 0):
ABSERR = ERR2 + ERR3
RESULT = RES
else:
RESULT = RES3LA[min(NRES-1, 2)]
break
RES = E1 + 1.0 / SS
EPSTAB[K1] = RES
K1 = K1 - 2
if (NRES == 0):
ABSERR = ERR2 + abs(RES - E2) + ERR3
RESULT = RES
continue
ERROR = self._compute_error(RES3LA, NRES, RES)

if (ERROR > 10.0 * ABSERR):
continue
ABSERR = ERROR
RESULT = RES
else:
ERROR = self._compute_error(RES3LA, NRES, RES)

# 50
if (N == self.limexp - 1):
N = 2 * (self.limexp // 2) - 1
EPSTAB = self._shift_table(EPSTAB, N, NEWELM, NUM)
self._update_RES3LA(RES3LA, RESULT, NRES)

ABSERR = max(ABSERR, 10.0*_EPS * abs(RESULT))
NRES = NRES + 1

N += 1
self._n = N
self._nres = NRES
# EPSTAB[-3:] = RES3LA
self.ABSERR = ABSERR
return RESULT, ABSERR


def test_dea():
def linfun(i):
return np.linspace(0, np.pi/2., 2**i+1)
dea = Dea(limexp=11)
print('NO. PANELS TRAP. APPROX APPROX W/EA ABSERR')
for k in np.arange(10):
x = linfun(k)
val = np.trapz(np.sin(x), x)
vale, err = dea(val)
print('%5d %20.8f %20.8f %20.8f' % (len(x)-1, val, vale, err))


def test_epsal():
HUGE = 1.E+60
TINY = 1.E-60
ZERO = 0.E0
ONE = 1.E0
true_vals = [0.78539816, 0.94805945, 0.99945672]
E = []
for N, SOFN in enumerate([0.78539816, 0.94805945, 0.98711580]):
E.append(SOFN)
if N == 0:
ESTLIM = SOFN
else:
AUX2 = ZERO
for J in range(N, 0, -1):
AUX1 = AUX2
AUX2 = E[J-1]
DIFF = E[J] - AUX2
if (abs(DIFF) <= TINY):
E[J-1] = HUGE
else:
E[J-1] = AUX1 + ONE/DIFF

if (N % 2) == 0:
ESTLIM = E[0]
else:
ESTLIM = E[1]
print(ESTLIM, true_vals[N])


def dea3(v0, v1, v2, symmetric=False):
"""
Extrapolate a slowly convergent sequence
Parameters
----------
v0, v1, v2 : array-like
3 values of a convergent sequence to extrapolate
Returns
-------
result : array-like
extrapolated value
abserr : array-like
absolute error estimate
Description
-----------
DEA3 attempts to extrapolate nonlinearly to a better estimate
of the sequence's limiting value, thus improving the rate of
convergence. The routine is based on the epsilon algorithm of
P. Wynn, see [1]_.
Example
-------
# integrate sin(x) from 0 to pi/2
>>> import numpy as np
>>> import numdifftools as nd
>>> Ei= np.zeros(3)
>>> linfun = lambda i : np.linspace(0, np.pi/2., 2**(i+5)+1)
>>> 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])
>>> truErr = Ei-1.
>>> (truErr, err, En)
(array([ -2.00805680e-04, -5.01999079e-05, -1.25498825e-05]),
array([ 0.00020081]), array([ 1.]))
See also
--------
dea
Reference
---------
.. [1] C. Brezinski (1977)
"Acceleration de la convergence en analyse numerique",
"Lecture Notes in Math.", vol. 584,
Springer-Verlag, New York, 1977.
"""
E0, E1, E2 = np.atleast_1d(v0, v1, v2)
abs, max = np.abs, np.maximum # @ReservedAssignment
with warnings.catch_warnings():
warnings.simplefilter("ignore") # ignore division by zero and overflow
delta2, delta1 = E2 - E1, E1 - E0
err2, err1 = abs(delta2), abs(delta1)
tol2, tol1 = max(abs(E2), abs(E1)) * _EPS, max(abs(E1), abs(E0)) * _EPS
delta1[err1 < _TINY] = _TINY
delta2[err2 < _TINY] = _TINY # avoid division by zero and overflow
ss = 1.0 / delta2 - 1.0 / delta1 + _TINY
smalle2 = (abs(ss * E1) <= 1.0e-3)
converged = (err1 <= tol1) & (err2 <= tol2) | smalle2
result = np.where(converged, E2 * 1.0, E1 + 1.0 / ss)
abserr = err1 + err2 + np.where(converged, tol2 * 10, abs(result-E2))
if symmetric and len(result) > 1:
return result[:-1], abserr[1:]
return result, abserr


def fornberg_weights_all(x, x0, M=1):
'''
Return finite difference weights_and_points for derivatives
Expand Down Expand Up @@ -1635,7 +1872,8 @@ def main2():
num_extrap = 0
method = 'complex'
data = []
for name in ['exp', 'expm1', 'sin', 'cos', 'square']: # function_names[:-3]:
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):
Expand Down
Loading

0 comments on commit 2a88fb2

Please sign in to comment.