Removed duplicated dea3 from integrate and misc.py import from numdifftools.extrapolate.dea3 instead

pepified
master
Per A Brodtkorb 9 years ago
parent 990e2a5151
commit d0867b5c72

@ -8,6 +8,8 @@ from scipy import special as sp
from .plotbackend import plotbackend as plt from .plotbackend import plotbackend as plt
from scipy.integrate import simps, trapz from scipy.integrate import simps, trapz
from .demos import humps from .demos import humps
from .misc import dea3
from .dctpack import dct
# from pychebfun import Chebfun # from pychebfun import Chebfun
_EPS = np.finfo(float).eps _EPS = np.finfo(float).eps
@ -18,79 +20,6 @@ __all__ = ['dea3', 'clencurt', 'romberg',
'gaussq', 'richardson', 'quadgr', 'qdemo'] 'gaussq', 'richardson', 'quadgr', 'qdemo']
def dea3(v0, v1, v2):
'''
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
>>> Ei= np.zeros(3)
>>> linfun = lambda k : np.linspace(0, np.pi/2., 2.**(k+5)+1)
>>> for k in np.arange(3):
... x = linfun(k)
... Ei[k] = np.trapz(np.sin(x),x)
>>> En, err = dea3(Ei[0],Ei[1],Ei[2])
>>> En, err
(array([ 1.]), array([ 0.0002008]))
>>> TrueErr = Ei-1.
>>> TrueErr
array([ -2.0080568e-04, -5.0199908e-05, -1.2549882e-05])
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 = np.abs # @ReservedAssignment
max = np.maximum # @ReservedAssignment
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
with warnings.catch_warnings():
warnings.simplefilter("ignore") # ignore division by zero and overflow
ss = 1.0 / delta2 - 1.0 / delta1
smalle2 = (abs(ss * E1) <= 1.0e-3).ravel()
result = 1.0 * E2
abserr = err1 + err2 + abs(E2) * _EPS * 10.0
converged = (err1 <= tol1) & (err2 <= tol2).ravel() | smalle2
k4, = (1 - converged).nonzero()
if k4.size > 0:
result[k4] = E1[k4] + 1.0 / ss[k4]
abserr[k4] = err1[k4] + err2[k4] + abs(result[k4] - E2[k4])
return result, abserr
def clencurt(fun, a, b, n0=5, trace=False, args=()): def clencurt(fun, a, b, n0=5, trace=False, args=()):
''' '''
Numerical evaluation of an integral, Clenshaw-Curtis method. Numerical evaluation of an integral, Clenshaw-Curtis method.
@ -183,13 +112,11 @@ def clencurt(fun, a, b, n0=5, trace=False, args=()):
fft = np.fft.fft fft = np.fft.fft
tmp = np.real(fft(f[:n, :], axis=0)) tmp = np.real(fft(f[:n, :], axis=0))
c = 2 / n * (tmp[0:n / 2 + 1, :] + np.cos(np.pi * s2) * f[n, :]) c = 2 / n * (tmp[0:n / 2 + 1, :] + np.cos(np.pi * s2) * f[n, :])
# % old call
# % c = 2/n * cos(s2*s'*pi/n) * f
c[0, :] = c[0, :] / 2 c[0, :] = c[0, :] / 2
c[n / 2, :] = c[n / 2, :] / 2 c[n / 2, :] = c[n / 2, :] / 2
# % alternative call # % alternative call
# % c = dct(f) # c2 = dct(f)
c = c[0:n / 2 + 1, :] / ((s2 - 1) * (s2 + 1)) c = c[0:n / 2 + 1, :] / ((s2 - 1) * (s2 + 1))
Q = (af - bf) * np.sum(c, axis=0) Q = (af - bf) * np.sum(c, axis=0)

@ -1102,10 +1102,10 @@ def test_smoothing_spline():
dy1 = pp1(x1) dy1 = pp1(x1)
y01 = pp0(x1) y01 = pp0(x1)
# dy = y-y1 # dy = y-y1
import matplotlib.pyplot as plb import matplotlib.pyplot as plt
plb.plot(x, y, x1, y1, '.', x1, dy1, 'ro', x1, y01, 'r-') plt.plot(x, y, x1, y1, '.', x1, dy1, 'ro', x1, y01, 'r-')
plb.show('hold') plt.show('hold')
pass pass
# tck = interpolate.splrep(x, y, s=len(x)) # tck = interpolate.splrep(x, y, s=len(x))
@ -1250,10 +1250,10 @@ def test_pp():
pp(1) pp(1)
pp(1.5) pp(1.5)
dpp = pp.derivative() dpp = pp.derivative()
import pylab as plb import matplotlib.pyplot as plt
x = plb.linspace(-1, 3) x = np.linspace(-1, 3)
plb.plot(x, pp(x), x, dpp(x), '.') plt.plot(x, pp(x), x, dpp(x), '.')
plb.show() plt.show()
def test_docstrings(): def test_docstrings():

@ -17,6 +17,7 @@ from scipy.special import gammaln
from scipy.integrate import trapz, simps from scipy.integrate import trapz, simps
import warnings import warnings
from time import strftime, gmtime from time import strftime, gmtime
from numdifftools.extrapolation import dea3
from .plotbackend import plotbackend from .plotbackend import plotbackend
from collections import OrderedDict from collections import OrderedDict
try: try:
@ -33,7 +34,7 @@ __all__ = ['now', 'spaceline', 'narg_smallest', 'args_flat', 'is_numlike',
'parse_kwargs', 'detrendma', 'ecross', 'findcross', 'findextrema', 'parse_kwargs', 'detrendma', 'ecross', 'findcross', 'findextrema',
'findpeaks', 'findrfc', 'rfcfilter', 'findtp', 'findtc', 'findpeaks', 'findrfc', 'rfcfilter', 'findtp', 'findtc',
'findoutliers', 'common_shape', 'argsreduce', 'stirlerr', 'findoutliers', 'common_shape', 'argsreduce', 'stirlerr',
'getshipchar', 'getshipchar', 'dea3',
'betaloge', 'gravity', 'nextpow2', 'discretize', 'polar2cart', 'betaloge', 'gravity', 'nextpow2', 'discretize', 'polar2cart',
'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc',
'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere', 'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere',
@ -2107,79 +2108,6 @@ def gravity(phi=45):
0.0000059 * sin(2 * phir) ** 2.) 0.0000059 * sin(2 * phir) ** 2.)
def dea3(v0, v1, v2):
'''
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 k : np.linspace(0,np.pi/2.,2.**(k+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 = np.abs # @ReservedAssignment
max = np.maximum # @ReservedAssignment
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
with warnings.catch_warnings():
warnings.simplefilter("ignore") # ignore division by zero and overflow
ss = 1.0 / delta2 - 1.0 / delta1
smallE2 = (abs(ss * E1) <= 1.0e-3).ravel()
result = 1.0 * E2
abserr = err1 + err2 + abs(E2) * _EPS * 10.0
converged = (err1 <= tol1) & (err2 <= tol2).ravel() | smallE2
k4, = (1 - converged).nonzero()
if k4.size > 0:
result[k4] = E1[k4] + 1.0 / ss[k4]
abserr[k4] = err1[k4] + err2[k4] + abs(result[k4] - E2[k4])
return result, abserr
def nextpow2(x): def nextpow2(x):
''' '''
Return next higher power of 2 Return next higher power of 2

@ -724,7 +724,7 @@ def poly2str(p, variable='x'):
N = len(coeffs) - 1 N = len(coeffs) - 1
for k in range(len(coeffs)): for k in range(N+1):
coefstr = '%.4g' % abs(coeffs[k]) coefstr = '%.4g' % abs(coeffs[k])
if coefstr[-4:] == '0000': if coefstr[-4:] == '0000':
coefstr = coefstr[:-5] coefstr = coefstr[:-5]
@ -733,21 +733,18 @@ def poly2str(p, variable='x'):
if coefstr != '0': if coefstr != '0':
newstr = '%s' % (coefstr,) newstr = '%s' % (coefstr,)
else: else:
if k == 0: newstr = '0' if k == 0 else ''
newstr = '0'
else:
newstr = ''
elif power == 1: elif power == 1:
if coefstr == '0': if coefstr == '0':
newstr = '' newstr = ''
elif coefstr == 'b' or coefstr == '1': elif coefstr in ['b', '1']:
newstr = var newstr = var
else: else:
newstr = '%s*%s' % (coefstr, var) newstr = '%s*%s' % (coefstr, var)
else: else:
if coefstr == '0': if coefstr == '0':
newstr = '' newstr = ''
elif coefstr == 'b' or coefstr == '1': elif coefstr in ['b', '1']:
newstr = '%s**%d' % (var, power,) newstr = '%s**%d' % (var, power,)
else: else:
newstr = '%s*%s**%d' % (coefstr, var, power) newstr = '%s*%s**%d' % (coefstr, var, power)

@ -1811,7 +1811,7 @@ class SpecData1D(PlotData):
rind = Rind(**options) rind = Rind(**options)
if (Nx > 1): if (Nx > 1):
# (M,m) or (M,m)v distribution wanted # (M,m) or (M,m)v distribution wanted
if ((def_nr == 0 or def_nr == 2)): if def_nr in [0, 2]:
asize = [Nx1, Nx1] asize = [Nx1, Nx1]
else: else:
# (M,m,TMm), (M,m,TMm)v (M,m,TMd)v or (M,M,Tdm)v # (M,m,TMm), (M,m,TMm)v (M,m,TMd)v or (M,M,Tdm)v

Loading…
Cancel
Save