diff --git a/wafo/integrate.py b/wafo/integrate.py index 99ccc8d..02cca20 100644 --- a/wafo/integrate.py +++ b/wafo/integrate.py @@ -8,6 +8,8 @@ from scipy import special as sp from .plotbackend import plotbackend as plt from scipy.integrate import simps, trapz from .demos import humps +from .misc import dea3 +from .dctpack import dct # from pychebfun import Chebfun _EPS = np.finfo(float).eps @@ -18,79 +20,6 @@ __all__ = ['dea3', 'clencurt', 'romberg', '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=()): ''' 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 tmp = np.real(fft(f[:n, :], axis=0)) 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[n / 2, :] = c[n / 2, :] / 2 # % alternative call -# % c = dct(f) + # c2 = dct(f) c = c[0:n / 2 + 1, :] / ((s2 - 1) * (s2 + 1)) Q = (af - bf) * np.sum(c, axis=0) diff --git a/wafo/interpolate.py b/wafo/interpolate.py index 9f4eed2..a97292b 100644 --- a/wafo/interpolate.py +++ b/wafo/interpolate.py @@ -1102,10 +1102,10 @@ def test_smoothing_spline(): dy1 = pp1(x1) y01 = pp0(x1) # 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-') - plb.show('hold') + plt.plot(x, y, x1, y1, '.', x1, dy1, 'ro', x1, y01, 'r-') + plt.show('hold') pass # tck = interpolate.splrep(x, y, s=len(x)) @@ -1250,10 +1250,10 @@ def test_pp(): pp(1) pp(1.5) dpp = pp.derivative() - import pylab as plb - x = plb.linspace(-1, 3) - plb.plot(x, pp(x), x, dpp(x), '.') - plb.show() + import matplotlib.pyplot as plt + x = np.linspace(-1, 3) + plt.plot(x, pp(x), x, dpp(x), '.') + plt.show() def test_docstrings(): diff --git a/wafo/misc.py b/wafo/misc.py index c765aea..1d900f8 100644 --- a/wafo/misc.py +++ b/wafo/misc.py @@ -17,6 +17,7 @@ from scipy.special import gammaln from scipy.integrate import trapz, simps import warnings from time import strftime, gmtime +from numdifftools.extrapolation import dea3 from .plotbackend import plotbackend from collections import OrderedDict try: @@ -33,7 +34,7 @@ __all__ = ['now', 'spaceline', 'narg_smallest', 'args_flat', 'is_numlike', 'parse_kwargs', 'detrendma', 'ecross', 'findcross', 'findextrema', 'findpeaks', 'findrfc', 'rfcfilter', 'findtp', 'findtc', 'findoutliers', 'common_shape', 'argsreduce', 'stirlerr', - 'getshipchar', + 'getshipchar', 'dea3', 'betaloge', 'gravity', 'nextpow2', 'discretize', 'polar2cart', 'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere', @@ -2107,79 +2108,6 @@ def gravity(phi=45): 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): ''' Return next higher power of 2 diff --git a/wafo/polynomial.py b/wafo/polynomial.py index eb45614..dd8824a 100644 --- a/wafo/polynomial.py +++ b/wafo/polynomial.py @@ -724,7 +724,7 @@ def poly2str(p, variable='x'): N = len(coeffs) - 1 - for k in range(len(coeffs)): + for k in range(N+1): coefstr = '%.4g' % abs(coeffs[k]) if coefstr[-4:] == '0000': coefstr = coefstr[:-5] @@ -733,21 +733,18 @@ def poly2str(p, variable='x'): if coefstr != '0': newstr = '%s' % (coefstr,) else: - if k == 0: - newstr = '0' - else: - newstr = '' + newstr = '0' if k == 0 else '' elif power == 1: if coefstr == '0': newstr = '' - elif coefstr == 'b' or coefstr == '1': + elif coefstr in ['b', '1']: newstr = var else: newstr = '%s*%s' % (coefstr, var) else: if coefstr == '0': newstr = '' - elif coefstr == 'b' or coefstr == '1': + elif coefstr in ['b', '1']: newstr = '%s**%d' % (var, power,) else: newstr = '%s*%s**%d' % (coefstr, var, power) diff --git a/wafo/spectrum/core.py b/wafo/spectrum/core.py index 886f365..2f77fdc 100644 --- a/wafo/spectrum/core.py +++ b/wafo/spectrum/core.py @@ -1811,7 +1811,7 @@ class SpecData1D(PlotData): rind = Rind(**options) if (Nx > 1): # (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] else: # (M,m,TMm), (M,m,TMm)v (M,m,TMd)v or (M,M,Tdm)v