danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merge branch 'master' of https://github.com/wafo-project/pywafo.git

Browse Source
master
pbrod 9 years ago
parent c4e5bfd66a 52411de937
commit d04c28df0d
4 changed files with 2377 additions and 2362 deletions
Show Stats Download Patch File Download Diff File

11
wafo/definitions.py
Unescape Escape View File

@ -12,11 +12,12 @@ Examples
-------- --------
In order to view the documentation do the following in an ipython window: In order to view the documentation do the following in an ipython window:
>>> import wafo.definitions as wd import wafo.definitions as wd
>>> wd.crossings() wd.crossings()
or or
>>> wd.crossings?
wd.crossings?
""" """
@ -303,3 +304,7 @@ def waves():
findcross findcross
""" """
print(waves.__doc__) print(waves.__doc__)
if __name__ == '__main__':
import doctest
doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)

3
wafo/gaussian.py
Unescape Escape View File

@ -598,7 +598,8 @@ def prbnormndpc(rho, a, b, abserr=1e-4, relerr=1e-4, usesimpson=True,
''' '''
# Call fortran implementation # Call fortran implementation
val, err, ier = mvnprdmod.prbnormndpc(rho, a, b, abserr, relerr, usebreakpoints, usesimpson) # @UndefinedVariable @IgnorePep8 val, err, ier = mvnprdmod.prbnormndpc(rho, a, b, abserr, relerr,
usebreakpoints, usesimpson)
if ier > 0: if ier > 0:
warnings.warn('Abnormal termination ier = %d\n\n%s' % warnings.warn('Abnormal termination ier = %d\n\n%s' %

119
wafo/integrate.py
Unescape Escape View File

@ -5,8 +5,9 @@ from numpy import pi, sqrt, ones, zeros # @UnresolvedImport
from scipy import integrate as intg from scipy import integrate as intg
import scipy.special.orthogonal as ort import scipy.special.orthogonal as ort
from scipy import special as sp from scipy import special as sp
from .plotbackend import plotbackend as plt
from scipy.integrate import simps, trapz from scipy.integrate import simps, trapz
from .plotbackend import plotbackend as plt
from .demos import humps from .demos import humps
from .misc import dea3 from .misc import dea3
from .dctpack import dct from .dctpack import dct
@ -48,9 +49,9 @@ def clencurt(fun, a, b, n0=5, trace=False, args=()):
Example Example
------- -------
>>> import numpy as np >>> import numpy as np
>>> val,err = clencurt(np.exp,0,2) >>> val, err = clencurt(np.exp, 0, 2)
>>> abs(val-np.expm1(2))< err, err<1e-10 >>> np.allclose(val, np.expm1(2)), err[0] < 1e-10
(array([ True], dtype=bool), array([ True], dtype=bool)) (True, True)
See also See also
@ -103,24 +104,20 @@ def clencurt(fun, a, b, n0=5, trace=False, args=()):
f[0, :] = f[0, :] / 2 f[0, :] = f[0, :] / 2
f[n, :] = f[n, :] / 2 f[n, :] = f[n, :] / 2
# % x = cos(pi*0:n/n) # x = cos(pi*0:n/n)
# % f = f(x) # f = f(x)
# % #
# % N+1 # N+1
# % c(k) = (2/N) sum f''(n)*cos(pi*(2*k-2)*(n-1)/N), 1 <= k <= N/2+1. # c(k) = (2/N) sum f''(n)*cos(pi*(2*k-2)*(n-1)/N), 1 <= k <= N/2+1.
# % n=1 # n=1
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, :])
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
# 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)
# Q = (a-b).*sum( c(1:n/2+1,:)./repmat((s2-1).*(s2+1),1,Na))
abserr = (bf - af) * np.abs(c[n / 2, :]) abserr = (bf - af) * np.abs(c[n / 2, :])
@ -238,9 +235,9 @@ def h_roots(n, method='newton'):
Example Example
------- -------
>>> import numpy as np >>> import numpy as np
>>> [x,w] = h_roots(10) >>> x, w = h_roots(10)
>>> np.sum(x*w) >>> np.allclose(np.sum(x*w), -5.2516042729766621e-19)
-5.2516042729766621e-19 True
See also See also
-------- --------
@ -451,7 +448,7 @@ def la_roots(n, alpha=0, method='newton'):
>>> import numpy as np >>> import numpy as np
>>> [x,w] = h_roots(10) >>> [x,w] = h_roots(10)
>>> np.sum(x*w) >>> np.sum(x*w)
-5.2516042729766621e-19 1.3352627380516791e-17
See also See also
-------- --------
@ -555,9 +552,9 @@ def p_roots(n, method='newton', a=-1, b=1):
------- -------
Integral of exp(x) from a = 0 to b = 3 is: exp(3)-exp(0)= Integral of exp(x) from a = 0 to b = 3 is: exp(3)-exp(0)=
>>> import numpy as np >>> import numpy as np
>>> [x,w] = p_roots(11,a=0,b=3) >>> x, w = p_roots(11, a=0, b=3)
>>> np.sum(np.exp(x)*w) >>> np.allclose(np.sum(np.exp(x)*w), 19.085536923187668)
19.085536923187668 True
See also See also
-------- --------
@ -723,15 +720,22 @@ def qrule(n, wfun=1, alpha=0, beta=0):
Examples: Examples:
--------- ---------
>>> import numpy as np
# integral of x^2 from a = -1 to b = 1
>>> [bp,wf] = qrule(10) >>> [bp,wf] = qrule(10)
>>> sum(bp**2*wf) # integral of x^2 from a = -1 to b = 1 >>> np.allclose(sum(bp**2*wf), 0.66666666666666641)
0.66666666666666641 True
# integral of exp(-x.^2)*x.^2 from a = -inf to b = inf
>>> [bp,wf] = qrule(10,2) >>> [bp,wf] = qrule(10,2)
>>> sum(bp**2*wf) # integral of exp(-x.^2)*x.^2 from a = -inf to b = inf >>> np.allclose(sum(bp**2*wf), 0.88622692545275772)
0.88622692545275772 True
# integral of (x+1)*(1-x)^2 from a = -1 to b = 1
>>> [bp,wf] = qrule(10,4,1,2) >>> [bp,wf] = qrule(10,4,1,2)
>>> (bp*wf).sum() # integral of (x+1)*(1-x)^2 from a = -1 to b = 1 >>> np.allclose((bp*wf).sum(), 0.26666666666666755)
0.26666666666666755 True
See also See also
-------- --------
@ -841,23 +845,24 @@ class _Gaussq(object):
--------- ---------
integration of x**2 from 0 to 2 and from 1 to 4 integration of x**2 from 0 to 2 and from 1 to 4
>>> from scitools import numpyutils as npt >>> import numpy as np
>>> A = [0, 1]; B = [2,4] >>> A = [0, 1]
>>> fun = npt.wrap2callable('x**2') >>> B = [2, 4]
>>> [val1,err1] = gaussq(fun,A,B) >>> fun = lambda x: x**2
>>> val1 >>> val1, err1 = gaussq(fun,A,B)
array([ 2.6666667, 21. ]) >>> np.allclose(val1, [ 2.6666667, 21. ])
>>> err1 True
array([ 1.7763568e-15, 1.0658141e-14]) >>> np.allclose(err1, [ 1.7763568e-15, 1.0658141e-14])
True
Integration of x^2*exp(-x) from zero to infinity: Integration of x^2*exp(-x) from zero to infinity:
>>> fun2 = npt.wrap2callable('1') >>> fun2 = lambda x : np.ones(np.shape(x))
>>> val2, err2 = gaussq(fun2, 0, npt.inf, wfun=3, alpha=2) >>> val2, err2 = gaussq(fun2, 0, np.inf, wfun=3, alpha=2)
>>> val3, err3 = gaussq(lambda x: x**2,0, npt.inf, wfun=3, alpha=0) >>> val3, err3 = gaussq(lambda x: x**2,0, np.inf, wfun=3, alpha=0)
>>> val2, err2 >>> np.allclose(val2, 2), err2[0] < 1e-14
(array([ 2.]), array([ 6.6613381e-15])) (True, True)
>>> val3, err3 >>> np.allclose(val3, 2), err3[0] < 1e-14
(array([ 2.]), array([ 1.7763568e-15])) (True, True)
Integrate humps from 0 to 2 and from 1 to 4 Integrate humps from 0 to 2 and from 1 to 4
>>> val4, err4 = gaussq(humps,A,B) >>> val4, err4 = gaussq(humps,A,B)
@ -1024,23 +1029,29 @@ class _Quadgr(object):
-------- --------
>>> import numpy as np >>> import numpy as np
>>> Q, err = quadgr(np.log,0,1) >>> Q, err = quadgr(np.log,0,1)
>>> quadgr(np.exp,0,9999*1j*np.pi) >>> q, err = quadgr(np.exp,0,9999*1j*np.pi)
(-2.0000000000122662, 2.1933237448479304e-09) >>> np.allclose(q, -2.0000000000122662), err < 1.0e-08
(True, True)
>>> quadgr(lambda x: np.sqrt(4-x**2),0,2,1e-12) >>> q, err = quadgr(lambda x: np.sqrt(4-x**2), 0, 2, abseps=1e-12)
(3.1415926535897811, 1.5809575870662229e-13) >>> np.allclose(q, 3.1415926535897811), err < 1.0e-12
(True, True)
>>> quadgr(lambda x: x**-0.75,0,1) >>> q, err = quadgr(lambda x: x**-0.75, 0, 1)
(4.0000000000000266, 5.6843418860808015e-14) >>> np.allclose(q, 4), err < 1.e-13
(True, True)
>>> quadgr(lambda x: 1./np.sqrt(1-x**2),-1,1) >>> q, err = quadgr(lambda x: 1./np.sqrt(1-x**2), -1, 1)
(3.141596056985029, 6.2146261559092864e-06) >>> np.allclose(q, 3.141596056985029), err < 1.0e-05
(True, True)
>>> quadgr(lambda x: np.exp(-x**2),-np.inf,np.inf,1e-9) #% sqrt(pi) >>> q, err = quadgr(lambda x: np.exp(-x**2), -np.inf, np.inf, 1e-9)
(1.7724538509055152, 1.9722334876348668e-11) >>> np.allclose(q, np.sqrt(np.pi)), err < 1e-9
(True, True)
>>> quadgr(lambda x: np.cos(x)*np.exp(-x),0,np.inf,1e-9) >>> q, err = quadgr(lambda x: np.cos(x)*np.exp(-x), 0, np.inf, 1e-9)
(0.50000000000000044, 7.3296813063450372e-11) >>> np.allclose(q, 0.5), err < 1e-9
(True, True)
See also See also
-------- --------

176
wafo/polynomial.py
Unescape Escape View File

@ -1,7 +1,7 @@
""" """
Extended functions to operate on polynomials Extended functions to operate on polynomials
""" """
# ------------------------------------------------------------------------- # ------------------------------------------------------------------------
# Name: polynomial # Name: polynomial
# Purpose: Functions to operate on polynomials. # Purpose: Functions to operate on polynomials.
# #
@ -15,21 +15,17 @@
# Created: 30.12.2008 # Created: 30.12.2008
# Copyright: (c) pab 2008 # Copyright: (c) pab 2008
# Licence: LGPL # Licence: LGPL
# ------------------------------------------------------------------------- # ------------------------------------------------------------------------
# !/usr/bin/env python # !/usr/bin/env python
from __future__ import absolute_import from __future__ import absolute_import
import warnings # @UnusedImport import warnings # @UnusedImport
from numpy.polynomial import polyutils as pu from numpy.polynomial import polyutils as pu
from .plotbackend import plotbackend as plt from .plotbackend import plotbackend as plt
import numpy as np import numpy as np
from numpy import (zeros, asarray, newaxis, arange, from numpy import (newaxis, arange, pi)
logical_or, any, pi, cos, round, diff, all, exp,
where, extract, linalg, sign, concatenate, floor,
linspace, sum, meshgrid)
from scipy.fftpack import dct, idct as _idct from scipy.fftpack import dct, idct as _idct
from numpy.lib.polynomial import * # @UnusedWildImport from numpy.lib.polynomial import * # @UnusedWildImport
from scipy.misc.common import pade # @UnresolvedImport from scipy.misc import pade # @UnresolvedImport
__all__ = np.lib.polynomial.__all__ __all__ = np.lib.polynomial.__all__
__all__ = __all__ + ['pade', 'padefit', 'polyreloc', 'polyrescl', 'polytrim', __all__ = __all__ + ['pade', 'padefit', 'polyreloc', 'polyrescl', 'polytrim',
'poly2hstr', 'poly2str', 'polyshift', 'polyishift', 'poly2hstr', 'poly2str', 'polyshift', 'polyishift',
@ -91,9 +87,8 @@ def polyint(p, m=1, k=None):
>>> np.polyder(P, 2)(0) >>> np.polyder(P, 2)(0)
0.0 0.0
>>> P = np.polyint(p, 3, k=[6, 5, 3]) >>> P = np.polyint(p, 3, k=[6, 5, 3])
>>> P.coefficients.tolist() >>> P
[0.016666666666666666, 0.041666666666666664, 0.16666666666666666, 3.0, poly1d([ 0.01666667, 0.04166667, 0.16666667, 3. , 5. , 3. ])
5.0, 3.0]
Note that 3 = 6 / 2!, and that the constants are given in the order of Note that 3 = 6 / 2!, and that the constants are given in the order of
integrations. Constant of the highest-order polynomial term comes first: integrations. Constant of the highest-order polynomial term comes first:
@ -110,7 +105,7 @@ def polyint(p, m=1, k=None):
if m < 0: if m < 0:
raise ValueError("Order of integral must be positive (see polyder)") raise ValueError("Order of integral must be positive (see polyder)")
if k is None: if k is None:
k = zeros(m, float) k = np.zeros(m, float)
k = np.atleast_1d(k) k = np.atleast_1d(k)
if len(k) == 1 and m > 1: if len(k) == 1 and m > 1:
k = k[0] * np.ones(m, float) k = k[0] * np.ones(m, float)
@ -118,7 +113,7 @@ def polyint(p, m=1, k=None):
raise ValueError( raise ValueError(
"k must be a scalar or a rank-1 array of length 1 or >m.") "k must be a scalar or a rank-1 array of length 1 or >m.")
truepoly = isinstance(p, poly1d) truepoly = isinstance(p, poly1d)
p = asarray(p) p = np.asarray(p)
if m == 0: if m == 0:
if truepoly: if truepoly:
return poly1d(p) return poly1d(p)
@ -195,7 +190,7 @@ def polyder(p, m=1):
if m < 0: if m < 0:
raise ValueError("Order of derivative must be positive (see polyint)") raise ValueError("Order of derivative must be positive (see polyint)")
truepoly = isinstance(p, poly1d) truepoly = isinstance(p, poly1d)
p = asarray(p) p = np.asarray(p)
if m == 0: if m == 0:
if truepoly: if truepoly:
return poly1d(p) return poly1d(p)
@ -271,12 +266,12 @@ def polydeg(x, y):
# required to take AIC noise into account and to ensure that this minimum # required to take AIC noise into account and to ensure that this minimum
# is a (likely) global minimum. # is a (likely) global minimum.
while nit < 3: while nit < 6:
p = orthofit(x, y, n) p = orthofit(x, y, n)
ys = orthoval(p, x) ys = orthoval(p, x)
# -- Akaike's Information Criterion # -- Akaike's Information Criterion
aic = (2 * (n + 1) * (1 + (n + 2) / (N - n - 2)) + aic = (2 * (n + 1) * (1 + (n + 2) / (N - n - 2)) +
N * (np.log(2 * pi * sum((ys - y) ** 2) / N) + 1)) N * (np.log(2 * pi * np.sum((ys - y) ** 2) / N) + 1))
if aic >= AIC: if aic >= AIC:
nit += 1 nit += 1
@ -444,10 +439,10 @@ def orthofit(x, y, n):
PL[1] = x - p[1, 1] PL[1] = x - p[1, 1]
for i in range(2, n + 1): for i in range(2, n + 1):
p[1, i] = np.dot(x, PL[i - 1] ** 2) / sum(PL[i - 1] ** 2) p[1, i] = np.dot(x, PL[i - 1] ** 2) / np.sum(PL[i - 1] ** 2)
p[2, i] = np.dot(x, PL[i - 2] * PL[i - 1]) / sum(PL[i - 2] ** 2) p[2, i] = np.dot(x, PL[i - 2] * PL[i - 1]) / np.sum(PL[i - 2] ** 2)
PL[i] = (x - p[1, i]) * PL[i - 1] - p[2, i] * PL[i - 2] PL[i] = (x - p[1, i]) * PL[i - 1] - p[2, i] * PL[i - 2]
p[0, :] = np.dot(PL, y) / sum(PL ** 2, axis=1) p[0, :] = np.dot(PL, y) / np.sum(PL ** 2, axis=1)
return p return p
# ys = np.dot(p[0, :], PL) # smoothed y # ys = np.dot(p[0, :], PL) # smoothed y
@ -594,11 +589,11 @@ def polytrim(p):
else: else:
r = np.atleast_1d(p).copy() r = np.atleast_1d(p).copy()
# Remove leading zeros # Remove leading zeros
is_not_lead_zeros = logical_or.accumulate(r != 0, axis=0) is_not_lead_zeros = np.logical_or.accumulate(r != 0, axis=0)
if r.ndim == 1: if r.ndim == 1:
r = r[is_not_lead_zeros] r = r[is_not_lead_zeros]
else: else:
is_not_lead_zeros = any(is_not_lead_zeros, axis=1) is_not_lead_zeros = np.any(is_not_lead_zeros, axis=1)
r = r[is_not_lead_zeros, :] r = r[is_not_lead_zeros, :]
return r return r
@ -952,8 +947,8 @@ def cheb2poly(ck, a=-1, b=1):
n = len(ck) n = len(ck)
b_Nmi = zeros(1) b_Nmi = np.zeros(1)
b_Nmip1 = zeros(1) b_Nmip1 = np.zeros(1)
y = np.r_[2 / (b - a), -(a + b) / (b - a)] y = np.r_[2 / (b - a), -(a + b) / (b - a)]
y2 = 2. * y y2 = 2. * y
@ -1003,7 +998,7 @@ def chebextr(n):
http://en.wikipedia.org/wiki/Chebyshev_nodes http://en.wikipedia.org/wiki/Chebyshev_nodes
http://en.wikipedia.org/wiki/Chebyshev_polynomials http://en.wikipedia.org/wiki/Chebyshev_polynomials
""" """
return - cos((pi * arange(n + 1)) / n) return - np.cos((pi * arange(n + 1)) / n)
def chebroot(n, kind=1): def chebroot(n, kind=1):
@ -1043,7 +1038,7 @@ def chebroot(n, kind=1):
""" """
if kind not in (1, 2): if kind not in (1, 2):
raise ValueError('kind must be 1 or 2') raise ValueError('kind must be 1 or 2')
return - cos(pi * (arange(n) + 0.5 * kind) / (n + kind - 1)) return - np.cos(pi * (arange(n) + 0.5 * kind) / (n + kind - 1))
def chebpoly(n, x=None, kind=1): def chebpoly(n, x=None, kind=1):
@ -1094,11 +1089,11 @@ def chebpoly(n, x=None, kind=1):
if n == 0: if n == 0:
p = np.ones(1) p = np.ones(1)
else: else:
p = round(pow(2, n - 2 + kind) * poly(chebroot(n, kind=kind))) p = np.round(pow(2, n - 2 + kind) * poly(chebroot(n, kind=kind)))
p[1::2] = 0 p[1::2] = 0
return p return p
else: # Evaluate polynomial in chebychev form else: # Evaluate polynomial in chebychev form
ck = zeros(n + 1) ck = np.zeros(n + 1)
ck[0] = 1. ck[0] = 1.
return _chebval(np.atleast_1d(x), ck, kind=kind) return _chebval(np.atleast_1d(x), ck, kind=kind)
@ -1137,12 +1132,12 @@ def chebfit(fun, n=10, a=-1, b=1, trace=False):
>>> a = 0; b = 2 >>> a = 0; b = 2
>>> ck = chebfit(np.exp,7,a,b); >>> ck = chebfit(np.exp,7,a,b);
>>> x = np.linspace(0,4); >>> x = np.linspace(0,4);
>>> h=plt.plot(x, np.exp(x), 'r', x, chebval(x,ck,a,b), 'g.')
>>> x1 = chebroot(9)*(b-a)/2+(b+a)/2 >>> x1 = chebroot(9)*(b-a)/2+(b+a)/2
>>> ck1 = chebfit(np.exp(x1)) >>> ck1 = chebfit(np.exp(x1))
>>> h = plt.plot(x,np.exp(x), 'r', x, chebval(x,ck1,a,b),'g.')
>>> plt.close() h=plt.plot(x, np.exp(x), 'r', x, chebval(x,ck,a,b), 'g.')
h = plt.plot(x,np.exp(x), 'r', x, chebval(x,ck1,a,b),'g.')
plt.close()
See also See also
-------- --------
@ -1283,7 +1278,7 @@ def idct(x, n=None):
Examples Examples
-------- --------
>>> import numpy as np >>> import numpy as np
>>> x = np.arange(5) >>> x = np.arange(5)*1.0
>>> np.abs(x-idct(dct(x)))<1e-14 >>> np.abs(x-idct(dct(x)))<1e-14
array([ True, True, True, True, True], dtype=bool) array([ True, True, True, True, True], dtype=bool)
>>> np.abs(x-dct(idct(x)))<1e-14 >>> np.abs(x-dct(idct(x)))<1e-14
@ -1320,7 +1315,7 @@ def _chebval(x, ck, kind=1):
http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html
""" """
n = len(ck) n = len(ck)
b_Nmi = zeros(x.shape) # b_(N-i) b_Nmi = np.zeros(x.shape) # b_(N-i)
b_Nmip1 = b_Nmi.copy() # b_(N-i+1) b_Nmip1 = b_Nmi.copy() # b_(N-i+1)
x2 = 2 * x x2 = 2 * x
# Clenshaw reccurence # Clenshaw reccurence
@ -1363,15 +1358,19 @@ def chebval(x, ck, a=-1, b=1, kind=1, fill=None):
>>> import matplotlib.pyplot as plt >>> import matplotlib.pyplot as plt
>>> x = np.linspace(-1,1) >>> x = np.linspace(-1,1)
>>> ck = np.zeros(5); ck[-1]=1 >>> ck = np.zeros(5); ck[-1]=1
>>> h = plt.plot(x,chebval(x,ck),x,chebpoly(4,x),'.') >>> y = chebval(x,ck)
>>> plt.close()
h = plt.plot(x, y, x, chebpoly(4,x),'.')
plt.close()
Fit exponential function: Fit exponential function:
>>> import matplotlib.pyplot as plt >>> import matplotlib.pyplot as plt
>>> ck = chebfit(np.exp,7,0,2) >>> ck = chebfit(np.exp,7,0,2)
>>> x = np.linspace(0,4); >>> x = np.linspace(0,4);
>>> h=plt.plot(x,chebval(x,ck,0,2),'g',x,np.exp(x)) >>> y2 = chebval(x,ck,0,2)
>>> plt.close()
h=plt.plot(x, y2, 'g', x, np.exp(x))
plt.close()
See also See also
-------- --------
@ -1383,14 +1382,14 @@ def chebval(x, ck, a=-1, b=1, kind=1, fill=None):
http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html
""" """
y = map_from_interval(atleast_1d(x), a, b) y = map_from_interval(np.atleast_1d(x), a, b)
if fill is None: if fill is None:
f = _chebval(y, ck, kind=kind) f = _chebval(y, ck, kind=kind)
else: else:
cond = (abs(y) <= 1) cond = (abs(y) <= 1)
f = where(cond, 0, fill) f = np.where(cond, 0, fill)
if any(cond): if np.any(cond):
yk = extract(cond, y) yk = np.extract(cond, y)
f[cond] = _chebval(yk, ck, kind=kind) f[cond] = _chebval(yk, ck, kind=kind)
return f return f
@ -1418,9 +1417,11 @@ def chebder(ck, a=-1, b=1):
>>> import matplotlib.pyplot as plt >>> import matplotlib.pyplot as plt
>>> ck = chebfit(np.exp,7,0,2) >>> ck = chebfit(np.exp,7,0,2)
>>> x = np.linspace(0,4) >>> x = np.linspace(0,4)
>>> ck2 = chebder(ck,0,2); >>> ck2 = chebder(ck,0,2)
>>> h = plt.plot(x,chebval(x,ck,0,2),'g',x,np.exp(x),'r') >>> y = chebval(x,ck2,0,2)
>>> plt.close()
h = plt.plot(x, y, 'g', x, np.exp(x), 'r')
plt.close()
See also See also
-------- --------
@ -1438,7 +1439,7 @@ def chebder(ck, a=-1, b=1):
""" """
n = len(ck) - 1 n = len(ck) - 1
cder = zeros(n, dtype=asarray(ck).dtype) cder = np.zeros(n, dtype=np.asarray(ck).dtype)
cder[0] = 2 * n * ck[0] cder[0] = 2 * n * ck[0]
cder[1] = 2 * (n - 1) * ck[1] cder[1] = 2 * (n - 1) * ck[1]
for j in range(2, n): for j in range(2, n):
@ -1470,8 +1471,10 @@ def chebint(ck, a=-1, b=1):
>>> ck = chebfit(np.exp,7,0,2) >>> ck = chebfit(np.exp,7,0,2)
>>> x = np.linspace(0,4) >>> x = np.linspace(0,4)
>>> ck2 = chebint(ck,0,2); >>> ck2 = chebint(ck,0,2);
>>> h=plt.plot(x,chebval(x,ck,0,2),'g',x,np.exp(x),'r.') >>> y =chebval(x,ck2,0,2)
>>> plt.close()
h=plt.plot(x,y,'g',x,np.exp(x),'r.')
plt.close()
See also See also
-------- --------
@ -1500,14 +1503,14 @@ def chebint(ck, a=-1, b=1):
n = len(ck) n = len(ck)
cint = zeros(n) cint = np.zeros(n)
con = 0.25 * (b - a) con = 0.25 * (b - a)
dif1 = diff(ck[-1::-2]) dif1 = np.diff(ck[-1::-2])
ix1 = np.r_[1:n - 1:2] ix1 = np.r_[1:n - 1:2]
cint[ix1] = -(con * dif1) / ix1 cint[ix1] = -(con * dif1) / ix1
if n > 3: if n > 3:
dif2 = diff(ck[-2::-2]) dif2 = np.diff(ck[-2::-2])
ix2 = np.r_[2:n - 1:2] ix2 = np.r_[2:n - 1:2]
cint[ix2] = -(con * dif2) / ix2 cint[ix2] = -(con * dif2) / ix2
cint = cint[::-1] cint = cint[::-1]
@ -1530,7 +1533,7 @@ class Cheb1d(object):
for key in ck.__dict__.keys(): for key in ck.__dict__.keys():
self.__dict__[key] = ck.__dict__[key] self.__dict__[key] = ck.__dict__[key]
return return
cki = trim_zeros(atleast_1d(ck), 'b') cki = trim_zeros(np.atleast_1d(ck), 'b')
if len(cki.shape) > 1: if len(cki.shape) > 1:
raise ValueError("Polynomial must be 1d only.") raise ValueError("Polynomial must be 1d only.")
self.__dict__['coeffs'] = cki self.__dict__['coeffs'] = cki
@ -1544,9 +1547,9 @@ class Cheb1d(object):
def __array__(self, t=None): def __array__(self, t=None):
if t: if t:
return asarray(self.coeffs, t) return np.asarray(self.coeffs, t)
else: else:
return asarray(self.coeffs) return np.asarray(self.coeffs)
def __repr__(self): def __repr__(self):
vals = repr(self.coeffs) vals = repr(self.coeffs)
@ -1590,11 +1593,11 @@ class Cheb1d(object):
def __eq__(self, other): def __eq__(self, other):
other = Cheb1d(other) other = Cheb1d(other)
return (all(self.coeffs == other.coeffs) and (self.a == other.a) and return (np.all(self.coeffs == other.coeffs) and (self.a == other.a) and
(self.b == other.b) and (self.kind == other.kind)) (self.b == other.b) and (self.kind == other.kind))
def __ne__(self, other): def __ne__(self, other):
return any(self.coeffs != other.coeffs) or (self.a != other.a) or ( return np.any(self.coeffs != other.coeffs) or (self.a != other.a) or (
self.b != other.b) or (self.kind != other.kind) self.b != other.b) or (self.kind != other.kind)
def __setattr__(self, key, val): def __setattr__(self, key, val):
@ -1631,8 +1634,8 @@ class Cheb1d(object):
if key < 0: if key < 0:
raise ValueError("Does not support negative powers.") raise ValueError("Does not support negative powers.")
if key > self.order: if key > self.order:
zr = zeros(key - self.order, self.coeffs.dtype) zr = np.zeros(key - self.order, self.coeffs.dtype)
self.__dict__['coeffs'] = concatenate((self.coeffs, zr)) self.__dict__['coeffs'] = np.concatenate((self.coeffs, zr))
self.__dict__['order'] = key self.__dict__['order'] = key
self.__dict__['coeffs'][key] = val self.__dict__['coeffs'][key] = val
return return
@ -1714,9 +1717,9 @@ def padefit(c, m=None):
poly1d([ 0.00277778, 0.03333333, 0.2 , 0.66666667, 1. ]) poly1d([ 0.00277778, 0.03333333, 0.2 , 0.66666667, 1. ])
poly1d([ 0.03333333, -0.33333333, 1. ]) poly1d([ 0.03333333, -0.33333333, 1. ])
>>> x = np.linspace(0,4); x = np.linspace(0,4)
>>> h = plt.plot(x,c(x),x,p(x)/q(x),'g-', x,np.exp(x),'r.') h = plt.plot(x,c(x),x,p(x)/q(x),'g-', x,np.exp(x),'r.')
>>> plt.close() plt.close()
See also See also
-------- --------
@ -1724,16 +1727,16 @@ def padefit(c, m=None):
""" """
if not m: if not m:
m = int(floor((len(c) - 1) * 0.5)) m = int(np.floor((len(c) - 1) * 0.5))
c = asarray(c) c = np.asarray(c)
return pade(c[::-1], m) return pade(c[::-1], m)
def test_pade(): def test_pade():
cof = array(([1.0, 1.0, 1.0 / 2, 1. / 6, 1. / 24])) cof = np.array(([1.0, 1.0, 1.0 / 2, 1. / 6, 1. / 24]))
p, q = pade(cof, 2) p, q = pade(cof, 2)
t = arange(0, 2, 0.1) t = arange(0, 2, 0.1)
assert(all(abs(p(t) / q(t) - exp(t)) < 0.3)) assert(np.all(abs(p(t) / q(t) - np.exp(t)) < 0.3))
def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True): def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True):
@ -1781,9 +1784,9 @@ def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True):
poly1d([ 0.01443847, 0.128842 , 0.55284547, 0.99999962]) poly1d([ 0.01443847, 0.128842 , 0.55284547, 0.99999962])
poly1d([-0.0049658 , 0.07610473, -0.44716929, 1. ]) poly1d([-0.0049658 , 0.07610473, -0.44716929, 1. ])
>>> x = np.linspace(0,4) x = np.linspace(0,4)
>>> h = plt.plot(x, polyval(c1,x)/polyval(c2,x),'g') h = plt.plot(x, polyval(c1,x)/polyval(c2,x),'g')
>>> h = plt.plot(x, np.exp(x), 'r') h = plt.plot(x, np.exp(x), 'r')
See also See also
-------- --------
@ -1822,9 +1825,8 @@ def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True):
fs = fun fs = fun
n = len(fs) n = len(fs)
if n < npt: if n < npt:
warnings.warn( warnings.warn('Check the result! Number of function values ' +
'Check the result! ' + 'should be at least: %d' % npt)
'Number of function values should be at least: %d' % npt)
if trace: if trace:
plt.plot(x, fs, '+') plt.plot(x, fs, '+')
@ -1833,11 +1835,11 @@ def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True):
ee = np.ones((npt)) ee = np.ones((npt))
mad = 0 mad = 0
u = zeros((npt, ncof)) u = np.zeros((npt, ncof))
for ix in range(MAXIT): for ix in range(MAXIT):
# Set up design matrix for least squares fit. # Set up design matrix for least squares fit.
pow1 = wt pow1 = wt
bb = pow1 * (fs + abs(mad) * sign(ee)) bb = pow1 * (fs + abs(mad) * np.sign(ee))
for jx in range(m + 1): for jx in range(m + 1):
u[:, jx] = pow1 u[:, jx] = pow1
@ -1848,8 +1850,8 @@ def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True):
pow1 = pow1 * x pow1 = pow1 * x
u[:, jx] = pow1 u[:, jx] = pow1
[u1, w, v] = linalg.svd(u, full_matrices=False) [u1, w, v] = np.linalg.svd(u, full_matrices=False)
cof = where(w == 0, 0.0, np.dot(bb, u1) / w) cof = np.where(w == 0, 0.0, np.dot(bb, u1) / w)
cof = np.dot(cof, v) cof = np.dot(cof, v)
# Tabulate the deviations and revise the weights # Tabulate the deviations and revise the weights
@ -1873,10 +1875,10 @@ def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True):
def main(): def main():
exp = np.exp
[c1, c2] = padefitlsq(exp, 3, 3, 0, 2) [c1, c2] = padefitlsq(exp, 3, 3, 0, 2)
x = linspace(0, 4) x = np.linspace(0, 4)
plt.plot(x, polyval(c1, x) / polyval(c2, x), 'g') plt.plot(x, polyval(c1, x) / polyval(c2, x), 'g')
plt.plot(x, exp(x), 'r') plt.plot(x, exp(x), 'r')
@ -1887,15 +1889,10 @@ def main():
_pr = polyreloc(p, 2) _pr = polyreloc(p, 2)
_pd = polyder(p) _pd = polyder(p)
_st = poly2str(p) _st = poly2str(p)
c = poly1d( # polynomial coeff exponential function:
1. / c = poly1d(1. / sp.gamma(np.r_[6 + 1:0:-1]))
sp.gamma(
np.r_[
6 +
1:0:-
1])) # polynomial coeff exponential function
[p, q] = padefit(c) [p, q] = padefit(c)
x = linspace(0, 4) x = np.linspace(0, 4)
plt.plot(x, c(x), x, p(x) / q(x), 'g-', x, exp(x), 'r.') plt.plot(x, c(x), x, p(x) / q(x), 'g-', x, exp(x), 'r.')
plt.close() plt.close()
x = arange(4) x = arange(4)
@ -1906,7 +1903,7 @@ def main():
b = 2 b = 2
ck = chebfit(exp, 6, a, b) ck = chebfit(exp, 6, a, b)
_t = chebval(0, ck, a, b) _t = chebval(0, ck, a, b)
x = linspace(0, 2, 6) x = np.linspace(0, 2, 6)
plt.plot(x, exp(x), 'r', x, chebval(x, ck, a, b), 'g.') plt.plot(x, exp(x), 'r', x, chebval(x, ck, a, b), 'g.')
# x1 = chebroot(9).'*(b-a)/2+(b+a)/2 ; # x1 = chebroot(9).'*(b-a)/2+(b+a)/2 ;
# ck1 =chebfit([x1 exp(x1)],9,a,b); # ck1 =chebfit([x1 exp(x1)],9,a,b);
@ -1932,7 +1929,7 @@ def test_polydeg():
n = polydeg(x, y) n = polydeg(x, y)
# n = 2 # n = 2
p = orthofit(x, y, n) p = orthofit(x, y, n)
xi = linspace(x.min(), x.max()) xi = np.linspace(x.min(), x.max())
ys0 = orthoval(p, x) ys0 = orthoval(p, x)
ys = orthoval(p, xi) ys = orthoval(p, xi)
@ -1947,7 +1944,8 @@ def test_polydeg():
def test_docstrings(): def test_docstrings():
import doctest import doctest
print('Testing docstrings in %s' % __file__) print('Testing docstrings in %s' % __file__)
doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE) options = doctest.NORMALIZE_WHITESPACE | doctest.ELLIPSIS
doctest.testmod(optionflags=options, verbose=False)
def chebvandernd(deg, *xi): def chebvandernd(deg, *xi):
@ -2280,7 +2278,7 @@ def test_chebfit2d():
n = 3 n = 3
xorder, yorder = n-1, n-1 xorder, yorder = n-1, n-1
x = chebroot(n=n, kind=1) x = chebroot(n=n, kind=1)
xgrid, ygrid = meshgrid(x, x) xgrid, ygrid = np.meshgrid(x, x)
def f(x, y): def f(x, y):
return np.exp(-x**2-6*y**2) return np.exp(-x**2-6*y**2)
@ -2314,6 +2312,6 @@ if __name__ == '__main__':
if False: # True: # if False: # True: #
main() main()
else: else:
test_chebfit2d() # test_chebfit2d()
# test_docstrings() test_docstrings()
# test_polydeg() # test_polydeg()

Write Preview
Loading…
Cancel
Save