diff --git a/.travis.yml b/.travis.yml index 274733d..51faa83 100644 --- a/.travis.yml +++ b/.travis.yml @@ -32,6 +32,7 @@ before_install: - conda config --add channels https://conda.anaconda.org/omnia - conda config --add channels https://conda.anaconda.org/pbrod - source activate condaenv + - sudo apt-get update - sudo apt-get install gfortran # Install packages install: diff --git a/wafo/dctpack.py b/wafo/dctpack.py index c7cb192..aa38bcc 100644 --- a/wafo/dctpack.py +++ b/wafo/dctpack.py @@ -432,7 +432,9 @@ def test_docstrings(): def test_commands(): import commands - commands.getstatusoutput('preprocess -DFORMAT=html -DDEVICE=screen tutorial.do.txt > tmp_preprocess__tutorial.do.txt') + commands.getstatusoutput('preprocess -DFORMAT=html -DDEVICE=screen ' + + 'tutorial.do.txt > ' + + 'tmp_preprocess__tutorial.do.txt') if __name__ == '__main__': diff --git a/wafo/kdetools.py b/wafo/kdetools.py index eb29cc9..ef4ae86 100644 --- a/wafo/kdetools.py +++ b/wafo/kdetools.py @@ -1841,8 +1841,8 @@ class Kernel(object): return self.hns(data) name = self.name[:4].lower() if name == 'epan': # Epanechnikov kernel - a = (8.0 * (d + 4.0) * (2 * sqrt(pi)) - ** d / sphere_volume(d)) ** (1. / (4.0 + d)) + a = (8.0 * (d + 4.0) * (2 * sqrt(pi)) ** d / + sphere_volume(d)) ** (1. / (4.0 + d)) elif name == 'biwe': # Bi-weight kernel a = 2.7779 if d > 2: @@ -3318,8 +3318,8 @@ def kde_demo4(N=50): f1.plot('b', label='hisj=%g' % kde1.hs) x = np.linspace(-4, 4) for loc in [-5, 5]: - plt.plot(x + loc, st.norm.pdf(x, 0, scale=1) - / 2, 'k:', label='True density') + plt.plot(x + loc, st.norm.pdf(x, 0, scale=1) / 2, 'k:', + label='True density') plt.legend() diff --git a/wafo/misc.py b/wafo/misc.py index b73e485..89b36ba 100644 --- a/wafo/misc.py +++ b/wafo/misc.py @@ -2,6 +2,7 @@ Misc ''' from __future__ import division +import collections import sys import fractions import numpy as np @@ -9,7 +10,7 @@ from numpy import ( meshgrid, abs, amax, any, logical_and, arange, linspace, atleast_1d, asarray, ceil, floor, frexp, hypot, - sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff, empty_like, + sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff, finfo, inf, pi, interp, isnan, isscalar, zeros, ones, linalg, r_, sign, unique, hstack, vstack, nonzero, where, extract) from scipy.special import gammaln, gamma, psi @@ -27,14 +28,187 @@ floatinfo = finfo(float) _TINY = np.finfo(float).tiny _EPS = np.finfo(float).eps -__all__ = [ - 'is_numlike', 'JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_select', - 'parse_kwargs', 'detrendma', 'ecross', 'findcross', - 'findextrema', 'findpeaks', 'findrfc', 'rfcfilter', 'findtp', 'findtc', - 'findoutliers', 'common_shape', 'argsreduce', - 'stirlerr', 'getshipchar', 'betaloge', 'gravity', 'nextpow2', - 'discretize', 'polar2cart', 'cart2polar', 'meshgrid', 'ndgrid', - 'trangood', 'tranproc', 'plot_histgrm', 'num2pistr', 'test_docstrings'] +__all__ = ['now', 'spaceline', 'narg_smallest', 'args_flat', 'is_numlike', + 'JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_select', + 'parse_kwargs', 'detrendma', 'ecross', 'findcross', 'findextrema', + 'findpeaks', 'findrfc', 'rfcfilter', 'findtp', 'findtc', + 'findoutliers', 'common_shape', 'argsreduce', 'stirlerr', + 'getshipchar', + 'betaloge', 'gravity', 'nextpow2', 'discretize', 'polar2cart', + 'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', + 'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere', + 'piecewise', + 'valarray'] + + +def valarray(shape, value=np.NaN, typecode=None): + """Return an array of all value. + """ + if typecode is None: + typecode = bool + out = ones(shape, dtype=typecode) * value + + if not isinstance(out, np.ndarray): + out = asarray(out) + return out + + +def piecewise(xi, condlist, funclist, fill_value=0.0, args=(), **kw): + """ + Evaluate a piecewise-defined function. + + Given a set of conditions and corresponding functions, evaluate each + function on the input data wherever its condition is true. + + Parameters + ---------- + xi : tuple + input arguments to the functions in funclist, i.e., (x0, x1,...., xn) + condlist : list of bool arrays + Each boolean array corresponds to a function in `funclist`. Wherever + `condlist[i]` is True, `funclist[i](x0,x1,...,xn)` is used as the + output value. Each boolean array in `condlist` selects a piece of `xi`, + and should therefore be of the same shape as `xi`. + + The length of `condlist` must correspond to that of `funclist`. + If one extra function is given, i.e. if + ``len(funclist) - len(condlist) == 1``, then that extra function + is the default value, used wherever all conditions are false. + funclist : list of callables, f(*(xi + args), **kw), or scalars + Each function is evaluated over `x` wherever its corresponding + condition is True. It should take an array as input and give an array + or a scalar value as output. If, instead of a callable, + a scalar is provided then a constant function (``lambda x: scalar``) is + assumed. + fill_value : scalar + fill value for out of range values. Default 0. + args : tuple, optional + Any further arguments given here passed to the functions + upon execution, i.e., if called ``piecewise(..., ..., args=(1, 'a'))``, + then each function is called as ``f(x0, x1,..., xn, 1, 'a')``. + kw : dict, optional + Keyword arguments used in calling `piecewise` are passed to the + functions upon execution, i.e., if called + ``piecewise(..., ..., lambda=1)``, then each function is called as + ``f(x0, x1,..., xn, lambda=1)``. + + Returns + ------- + out : ndarray + The output is the same shape and type as x and is found by + calling the functions in `funclist` on the appropriate portions of `x`, + as defined by the boolean arrays in `condlist`. Portions not covered + by any condition have undefined values. + + + See Also + -------- + choose, select, where + + Notes + ----- + This is similar to choose or select, except that functions are + evaluated on elements of `xi` that satisfy the corresponding condition from + `condlist`. + + The result is:: + + |-- + |funclist[0](x0[condlist[0]],x1[condlist[0]],...,xn[condlist[0]]) + out = |funclist[1](x0[condlist[1]],x1[condlist[1]],...,xn[condlist[1]]) + |... + |funclist[n2](x0[condlist[n2]],x1[condlist[n2]],...,xn[condlist[n2]]) + |-- + + Examples + -------- + Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``. + + >>> x = np.linspace(-2.5, 2.5, 6) + >>> piecewise(x, [x < 0, x >= 0], [-1, 1]) + array([-1., -1., -1., 1., 1., 1.]) + + Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for + ``x >= 0``. + + >>> piecewise((x,), [x < 0, x >= 0], [lambda x: -x, lambda x: x]) + array([ 2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) + + Define the absolute value, which is ``-x*y`` for ``x*y <0`` and ``x*y`` for + ``x*y >= 0`` + >>> X, Y = np.meshgrid(x, x) + >>> piecewise((X, Y), [X * Y < 0, ], + ... [lambda x, y: -x * y, lambda x, y: x * y]) + array([[ 6.25, 3.75, 1.25, 1.25, 3.75, 6.25], + [ 3.75, 2.25, 0.75, 0.75, 2.25, 3.75], + [ 1.25, 0.75, 0.25, 0.25, 0.75, 1.25], + [ 1.25, 0.75, 0.25, 0.25, 0.75, 1.25], + [ 3.75, 2.25, 0.75, 0.75, 2.25, 3.75], + [ 6.25, 3.75, 1.25, 1.25, 3.75, 6.25]]) + """ + def otherwise_condition(condlist): + return ~np.logical_or.reduce(condlist, axis=0) + + def check_shapes(condlist, funclist): + nc, nf = len(condlist), len(funclist) + if nc not in [nf-1, nf]: + raise ValueError("function list and condition list" + + " must be the same length") + + check_shapes(condlist, funclist) + if not isinstance(xi, tuple): + xi = (xi,) + + condlist = np.broadcast_arrays(*condlist) + if len(condlist) == len(funclist)-1: + condlist.append(otherwise_condition(condlist)) + + arrays = np.broadcast_arrays(*xi) + dtype = np.result_type(*arrays) + + out = valarray(arrays[0].shape, fill_value, dtype) + for cond, func in zip(condlist, funclist): + if isinstance(func, collections.Callable): + temp = tuple(np.extract(cond, arr) for arr in arrays) + args + np.place(out, cond, func(*temp, **kw)) + else: # func is a scalar value + np.place(out, cond, func) + return out + + +def lazywhere(cond, arrays, f, fillvalue=None, f2=None): + """ + np.where(cond, x, fillvalue) always evaluates x even where cond is False. + This one only evaluates f(arr1[cond], arr2[cond], ...). + For example, + >>> a, b = np.array([1, 2, 3, 4]), np.array([5, 6, 7, 8]) + >>> def f(a, b): + ... return a*b + >>> lazywhere(a > 2, (a, b), f, np.nan) + array([ nan, nan, 21., 32.]) + + Notice it assumes that all `arrays` are of the same shape, or can be + broadcasted together. + + """ + if fillvalue is None: + if f2 is None: + raise ValueError("One of (fillvalue, f2) must be given.") + else: + fillvalue = np.nan + else: + if f2 is not None: + raise ValueError("Only one of (fillvalue, f2) can be given.") + + arrays = np.broadcast_arrays(*arrays) + temp = tuple(np.extract(cond, arr) for arr in arrays) + out = valarray(np.shape(arrays[0]), value=fillvalue) + np.place(out, cond, f(*temp)) + if f2 is not None: + temp = tuple(np.extract(~cond, arr) for arr in arrays) + np.place(out, ~cond, f2(*temp)) + + return out def rotation_matrix(heading, pitch, roll): @@ -269,7 +443,7 @@ def index2sub(shape, index, order='C'): >>> a.ravel(order)[i] 23 >>> index2sub(shape, i, order=order) - (array([1]), array([2]), array([3])) + (1, 2, 3) See also -------- @@ -308,7 +482,7 @@ def sub2index(shape, *subscripts, **kwds): >>> a.ravel(order)[i] 23 >>> index2sub(shape, i, order=order) - (array([1]), array([2]), array([3])) + (1, 2, 3) See also -------- @@ -417,7 +591,8 @@ def sub_dict_select(somedict, somekeys): def parse_kwargs(options, **kwargs): - ''' Update options dict from keyword arguments if it exists in options + ''' + Update options dict from keyword arguments if the keyword exists in options Example >>> opt = dict(arg1=2, arg2=3) @@ -436,15 +611,6 @@ def parse_kwargs(options, **kwargs): return options -def testfun(*args, **kwargs): - opts = dict(opt1=1, opt2=2) - if (len(args) == 1 and len(kwargs) == 0 and type(args[0]) is str and - args[0].startswith('default')): - return opts - opts = parse_kwargs(opts, **kwargs) - return opts - - def detrendma(x, L): """ Removes a trend from data using a moving average @@ -484,7 +650,7 @@ def detrendma(x, L): if L != round(L): raise ValueError('L must be an integer') - x1 = atleast_1d(x) + x1 = np.atleast_1d(x) if x1.shape[0] == 1: x1 = x1.ravel() @@ -493,10 +659,10 @@ def detrendma(x, L): return x1 - x1.mean(axis=0) mn = x1[0:2 * L + 1].mean(axis=0) - y = empty_like(x1) + y = np.empty_like(x1) y[0:L] = x1[0:L] - mn - ix = r_[L:(n - L)] + ix = np.r_[L:(n - L)] trend = ((x1[ix + L] - x1[ix - L]) / (2 * L + 1)).cumsum(axis=0) + mn y[ix] = x1[ix] - trend y[n - L::] = x1[n - L::] - trend[-1] @@ -589,6 +755,13 @@ def _findcross(xn): return ind +def xor(a, b): + """ + Return True only when inputs differ. + """ + return a ^ b + + def findcross(x, v=0.0, kind=None): ''' Return indices to level v up and/or downcrossings of a vector @@ -627,10 +800,12 @@ def findcross(x, v=0.0, kind=None): >>> ind = wm.findcross(x,v) # all crossings >>> ind array([ 9, 25, 80, 97, 151, 168, 223, 239]) + >>> t0 = plt.plot(t,x,'.',t[ind],x[ind],'r.', t, ones(t.shape)*v) >>> ind2 = wm.findcross(x,v,'u') >>> ind2 array([ 9, 80, 151, 223]) + >>> t0 = plt.plot(t[ind2],x[ind2],'o') >>> plt.close('all') @@ -653,22 +828,19 @@ def findcross(x, v=0.0, kind=None): t_0 = int(xn[ind[0] + 1] < 0) ind = ind[t_0::2] elif kind in ('dw', 'uw', 'tw', 'cw'): - # make sure the first is a level v down-crossing if wdef=='dw' - # or make sure the first is a level v up-crossing if wdef=='uw' - # make sure the first is a level v down-crossing if wdef=='tw' - # or make sure the first is a level v up-crossing if - # wdef=='cw' - def xor(a, b): - return a ^ b + # make sure the first is a level v down-crossing + # if wdef=='dw' or wdef=='tw' + # or make sure the first is a level v up-crossing + # if wdef=='uw' or wdef=='cw' + first_is_down_crossing = int(xn[ind[0]] > xn[ind[0] + 1]) if xor(first_is_down_crossing, kind in ('dw', 'tw')): ind = ind[1::] - n_c = ind.size # number of level v crossings # make sure the number of troughs and crests are according to the # wavedef, i.e., make sure length(ind) is odd if dw or uw # and even if tw or cw - is_odd = mod(n_c, 2) + is_odd = mod(ind.size, 2) if xor(is_odd, kind in ('dw', 'uw')): ind = ind[:-1] else: @@ -696,6 +868,7 @@ def findextrema(x): >>> t = np.linspace(0,7*np.pi,250) >>> x = np.sin(t) >>> ind = wm.findextrema(x) + >>> a = plt.plot(t,x,'.',t[ind],x[ind],'r.') >>> plt.close('all') @@ -852,6 +1025,7 @@ def findrfc(tp, h=0.0, method='clib'): >>> x = np.sin(t)+0.1*np.sin(50*t) >>> ind = wm.findextrema(x) >>> ti, tp = t[ind], x[ind] + >>> a = plt.plot(t,x,'.',ti,tp,'r.') >>> ind1 = wm.findrfc(tp,0.3); ind1 array([ 0, 9, 32, 53, 74, 95, 116, 137]) @@ -957,15 +1131,14 @@ def findrfc(tp, h=0.0, method='clib'): def mctp2rfc(fmM, fMm=None): ''' - Return Rainflow matrix given a Markov matrix of a Markov chain - of turning points + Return Rainflow matrix given a Markov chain of turning points computes f_rfc = f_mM + F_mct(f_mM). Parameters ---------- fmM = the min2max Markov matrix, - fMm = the max2min Markov matrix, + fMm = the max2min Markov matrix, Returns ------- @@ -1020,8 +1193,7 @@ def mctp2rfc(fmM, fMm=None): SRA = RAA.sum() DRFC = SA - SRA - # ?? check - NT = min(mA[0] - sum(RAA[:, 0]), MA[0] - sum(RAA[0, :])) + NT = min(mA[0] - sum(RAA[:, 0]), MA[0] - sum(RAA[0, :])) # check! NT = max(NT, 0) # ??check if NT > 1e-6 * max(MA[0], mA[0]): @@ -1118,11 +1290,11 @@ def rfcfilter(x, h, method=0): Examples: --------- # 1. Filtered signal y is the turning points of x. - >>> import wafo.data + >>> import wafo.data as data >>> import wafo.misc as wm - >>> x = wafo.data.sea() + >>> x = data.sea() >>> y = wm.rfcfilter(x[:,1], h=0, method=1) - >>> np.all(np.abs(y[0:5]-np.array([-1.2004945 , 0.83950546, -0.09049454, + >>> np.all(np.abs(y[0:5]-np.array([-1.2004945 , 0.83950546, -0.09049454, ... -0.02049454, -0.09049454]))<1e-7) True >>> y.shape @@ -1132,7 +1304,7 @@ def rfcfilter(x, h, method=0): >>> y1 = wm.rfcfilter(x[:,1], h=0.5) >>> y1.shape (863,) - >>> np.all(np.abs(y1[0:5]-np.array([-1.2004945 , 0.83950546, -0.43049454, + >>> np.all(np.abs(y1[0:5]-np.array([-1.2004945 , 0.83950546, -0.43049454, ... 0.34950546, -0.51049454]))<1e-7) True >>> ind = wm.findtp(x[:,1], h=0.5) @@ -1306,8 +1478,6 @@ def findtp(x, h=0.0, kind=None): ind = ind[ind1] if kind in ('mw', 'Mw'): - def xor(a, b): - return a ^ b # make sure that the first is a Max if wdef == 'Mw' # or make sure that the first is a min if wdef == 'mw' first_is_max = (x[ind[0]] > x[ind[1]]) @@ -1351,7 +1521,6 @@ def findtc(x_in, v=None, kind=None): Example: -------- - >>> import wafo.data >>> import pylab as plt >>> import wafo.misc as wm >>> t = np.linspace(0,30,500).reshape((-1,1)) @@ -1454,7 +1623,6 @@ def findoutliers(x, zcrit=0.0, dcrit=None, ddcrit=None, verbose=False): Examples -------- >>> import numpy as np - >>> import wafo >>> import wafo.misc as wm >>> t = np.linspace(0,30,500).reshape((-1,1)) >>> xx = np.hstack((t, np.cos(t))) @@ -1958,8 +2126,9 @@ def betaloge(z, w): ''' # y = gammaln(z)+gammaln(w)-gammaln(z+w) zpw = z + w - return (stirlerr(z) + stirlerr(w) + 0.5 * log(2 * pi) + (w - 0.5) * log(w) - + (z - 0.5) * log(z) - stirlerr(zpw) - (zpw - 0.5) * log(zpw)) + return (stirlerr(z) + stirlerr(w) + 0.5 * log(2 * pi) + + (w - 0.5) * log(w) + (z - 0.5) * log(z) - stirlerr(zpw) - + (zpw - 0.5) * log(zpw)) # stirlings approximation: # (-(zpw-0.5).*log(zpw) +(w-0.5).*log(w)+(z-0.5).*log(z) +0.5*log(2*pi)) @@ -1995,7 +2164,6 @@ def gravity(phi=45): ... 9.79640552, 9.80629387]))<1.e-7 array([ True, True, True, True, True], dtype=bool) - See also -------- wdensity @@ -2088,623 +2256,99 @@ def dea3(v0, v1, v2): return result, abserr -def hyp2f1_taylor(a, b, c, z, tol=1e-13, itermax=500): - a, b, c, z = np.broadcast_arrays(*np.atleast_1d(a, b, c, z)) - shape = a.shape - ak, bk, ck, zk = [d.ravel() for d in (a, b, c, z)] - ajm1 = np.ones(ak.shape) - bjm2 = 0.5 * np.ones(ak.shape) - bjm1 = np.ones(ak.shape) - hout = np.zeros(ak.shape) - k0 = np.arange(len(ak)) - for j in range(0, itermax): - aj = ajm1 * (ak + j) * (bk + j) / (ck + j) * zk / (j + 1) - bj = bjm1 + aj - h, err = dea3(bjm2, bjm1, bj) - k = np.flatnonzero(err > tol * np.abs(h)) - hout[k0] = h - if len(k) == 0: - break - k0 = k0[k] - ak, bk, ck, zk = ak[k], bk[k], ck[k], zk[k] - ajm1 = aj[k] - bjm2 = bjm1[k] - bjm1 = bj[k] +def nextpow2(x): + ''' + Return next higher power of 2 + + Example + ------- + >>> import wafo.misc as wm + >>> wm.nextpow2(10) + 4 + >>> wm.nextpow2(np.arange(5)) + 3 + ''' + t = isscalar(x) or len(x) + if (t > 1): + f, n = frexp(t) else: - warnings.warn(('Reached %d limit! \n' + - '#%d values did not converge! Max error=%g') % - (j, len(k), np.max(err))) - return hout.reshape(shape) + f, n = frexp(abs(x)) + if (f == 0.5): + n = n - 1 + return n -def hyp2f1(a, b, c, z, rho=0.5): - e1 = gammaln(a) - e2 = gammaln(b) - e3 = gammaln(c) - e4 = gammaln(b - a) - e5 = gammaln(a - b) - e6 = gammaln(c - a) - e7 = gammaln(c - b) - e8 = gammaln(c - a - b) - e9 = gammaln(a + b - c) - _cmab = c - a - b - # ~(np.round(cmab) == cmab & cmab <= 0) - if abs(z) <= rho: - h = hyp2f1_taylor(a, b, c, z, 1e-15) - elif abs(1 - z) <= rho: # % Require that |arg(1-z)| 10: - break - xjm2 = xjm1 - xjm1 = xj + Example + ------- + >>> import wafo.misc as wm + >>> import numpy as np + >>> import pylab as plt + >>> x,y = wm.discretize(np.cos, 0, np.pi) + >>> xa,ya = wm.discretize(np.cos, 0, np.pi, method='adaptive') + >>> t = plt.plot(x, y, xa, ya, 'r.') + + plt.show() + >>> plt.close('all') + + ''' + if method.startswith('a'): + return _discretize_adaptive(fun, a, b, tol, n) else: - warnings.warn('Reached %d limit' % j) - return h + return _discretize_linear(fun, a, b, tol, n) -def hygfz(A, B, C, Z): - ''' Return hypergeometric function for a complex argument, F(a,b,c,z) +def _discretize_linear(fun, a, b, tol=0.005, n=5): + ''' + Automatic discretization of function, linear gridding + ''' + x = linspace(a, b, n) + y = fun(x) - Parameters - ---------- - a, b, c: - parameters where c <> 0,-1,-2,... - z :--- Complex argument + err0 = inf + err = 10000 + nmax = 2 ** 20 + while (err != err0 and err > tol and n < nmax): + err0 = err + x0 = x + y0 = y + n = 2 * (n - 1) + 1 + x = linspace(a, b, n) + y = fun(x) + y00 = interp(x, x0, y0) + err = 0.5 * amax(abs((y00 - y) / (abs(y00 + y) + _TINY))) + return x, y + + +def _discretize_adaptive(fun, a, b, tol=0.005, n=5): ''' - X = np.real(Z) - Y = np.imag(Z) - EPS = 1.0e-15 - L0 = C == np.round(C) and C < 0.0e0 - L1 = abs(1.0 - X) < EPS and Y == 0.0 and C - A - B <= 0.0 - L2 = abs(Z + 1.0) < EPS and abs(C - A + B - 1.0) < EPS - L3 = A == np.round(A) and A < 0.0 - L4 = B == np.round(B) and B < 0.0 - L5 = C - A == np.round(C - A) and C - A <= 0.0 - L6 = C - B == np.round(C - B) and C - B <= 0.0 - AA = A - BB = B - A0 = abs(Z) - if (A0 > 0.95): - EPS = 1.0e-8 - PI = 3.141592653589793 - EL = .5772156649015329 - if (L0 or L1): - # 'The hypergeometric series is divergent' - return np.inf - - NM = 0 - if (A0 == 0.0 or A == 0.0 or B == 0.0): - ZHF = 1.0 - elif (Z == 1.0 and C - A - B > 0.0): - GC = gamma(C) - GCAB = gamma(C - A - B) - GCA = gamma(C - A) - GCB = gamma(C - B) - ZHF = GC * GCAB / (GCA * GCB) - elif L2: - G0 = sqrt(PI) * 2.0 ** (-A) - G1 = gamma(C) - G2 = gamma(1.0 + A / 2.0 - B) - G3 = gamma(0.5 + 0.5 * A) - ZHF = G0 * G1 / (G2 * G3) - elif L3 or L4: - if (L3): - NM = int(np.round(abs(A))) - if (L4): - NM = int(np.round(abs(B))) - ZHF = 1.0 - ZR = 1.0 - for K in range(NM): - ZR = ZR * (A + K) * (B + K) / ((K + 1.) * (C + K)) * Z - ZHF = ZHF + ZR - elif L5 or L6: - if (L5): - NM = np.round(abs(C - A)) - if (L6): - NM = np.round(abs(C - B)) - ZHF = 1.0 + 0j - ZR = 1.0 + 0j - for K in range(NM): - ZR *= (C - A + K) * (C - B + K) / ((K + 1.) * (C + K)) * Z - ZHF = ZHF + ZR - ZHF = (1.0 - Z) ** (C - A - B) * ZHF - elif (A0 <= 1.0): - if (X < 0.0): - Z1 = Z / (Z - 1.0) - if (C > A and B < A and B > 0.0): - A = BB - B = AA - - ZC0 = 1.0 / ((1.0 - Z) ** A) - ZHF = 1.0 + 0j - ZR0 = 1.0 + 0j - ZW = 0 - for K in range(500): - ZR0 *= (A + K) * (C - B + K) / ((K + 1.0) * (C + K)) * Z1 - ZHF += ZR0 - if (abs(ZHF - ZW) < abs(ZHF) * EPS): - break - ZW = ZHF - ZHF = ZC0 * ZHF - elif (A0 >= 0.90): - ZW = 0.0 - GM = 0.0 - MCAB = np.round(C - A - B) - if (abs(C - A - B - MCAB) < EPS): - M = int(np.round(C - A - B)) - GA = gamma(A) - GB = gamma(B) - GC = gamma(C) - GAM = gamma(A + M) - GBM = gamma(B + M) - PA = psi(A) - PB = psi(B) - if (M != 0): - GM = 1.0 - for j in range(1, abs(M)): - GM *= j - RM = 1.0 - for j in range(1, abs(M) + 1): # DO 35 J=1,abs(M) - RM *= j - ZF0 = 1.0 - ZR0 = 1.0 - ZR1 = 1.0 - SP0 = 0.0 - SP = 0.0 - if (M >= 0): - ZC0 = GM * GC / (GAM * GBM) - ZC1 = -GC * (Z - 1.0) ** M / (GA * GB * RM) - for K in range(1, M): - ZR0 = ZR0 * \ - (A + K - 1.) * (B + K - 1.) / \ - (K * (K - M)) * (1. - Z) - ZF0 = ZF0 + ZR0 - for K in range(M): - SP0 = SP0 + 1.0 / \ - (A + K) + 1.0 / (B + K) - 1. / (K + 1.) - ZF1 = PA + PB + SP0 + 2.0 * EL + np.log(1.0 - Z) - for K in range(1, 501): - SP = SP + \ - (1.0 - A) / (K * (A + K - 1.0)) + ( - 1.0 - B) / (K * (B + K - 1.0)) - SM = 0.0 - for J in range(1, M): - SM += (1.0 - A) / ( - (J + K) * (A + J + K - 1.0)) + \ - 1.0 / (B + J + K - 1.0) - - ZP = PA + PB + 2.0 * EL + SP + SM + np.log(1.0 - Z) - ZR1 = ZR1 * \ - (A + M + K - 1.0) * (B + M + K - 1.0) / ( - K * (M + K)) * (1.0 - Z) - ZF1 = ZF1 + ZR1 * ZP - if (abs(ZF1 - ZW) < abs(ZF1) * EPS): - break - ZW = ZF1 - ZHF = ZF0 * ZC0 + ZF1 * ZC1 - elif (M < 0): - M = -M - ZC0 = GM * GC / (GA * GB * (1.0 - Z) ** M) - ZC1 = -(-1) ** M * GC / (GAM * GBM * RM) - for K in range(1, M): - ZR0 = ZR0 * \ - (A - M + K - 1.0) * (B - M + K - 1.0) / ( - K * (K - M)) * (1.0 - Z) - ZF0 = ZF0 + ZR0 - for K in range(1, M + 1): - SP0 = SP0 + 1.0 / K - ZF1 = PA + PB - SP0 + 2.0 * EL + np.log(1.0 - Z) - for K in range(1, 501): - SP = SP + \ - (1.0 - A) / (K * (A + K - 1.0)) + ( - 1.0 - B) / (K * (B + K - 1.0)) - SM = 0.0 - for J in range(1, M + 1): - SM = SM + 1.0 / (J + K) - ZP = PA + PB + 2.0 * EL + SP - SM + np.log(1.0 - Z) - ZR1 = ZR1 * \ - (A + K - 1.) * (B + K - 1.) / \ - (K * (M + K)) * (1. - Z) - ZF1 = ZF1 + ZR1 * ZP - if (abs(ZF1 - ZW) < abs(ZF1) * EPS): - break - ZW = ZF1 - ZHF = ZF0 * ZC0 + ZF1 * ZC1 - else: - GA = gamma(A) - GB = gamma(B) - GC = gamma(C) - GCA = gamma(C - A) - GCB = gamma(C - B) - GCAB = gamma(C - A - B) - GABC = gamma(A + B - C) - ZC0 = GC * GCAB / (GCA * GCB) - ZC1 = GC * GABC / (GA * GB) * (1.0 - Z) ** (C - A - B) - ZHF = 0 + 0j - ZR0 = ZC0 - ZR1 = ZC1 - for K in range(1, 501): - ZR0 = ZR0 * \ - (A + K - 1.) * (B + K - 1.) / \ - (K * (A + B - C + K)) * (1. - Z) - ZR1 = ZR1 * \ - (C - A + K - 1.0) * (C - B + K - 1.0) / ( - K * (C - A - B + K)) * (1.0 - Z) - ZHF = ZHF + ZR0 + ZR1 - if (abs(ZHF - ZW) < abs(ZHF) * EPS): - break - ZW = ZHF - ZHF = ZHF + ZC0 + ZC1 - else: - ZW = 0.0 - Z00 = 1.0 # + 0j - if (C - A < A and C - B < B): - Z00 = (1.0 - Z) ** (C - A - B) - A = C - A - B = C - B - ZHF = 1.0 - ZR = 1.0 - for K in range(1, 501): - ZR = ZR * \ - (A + K - 1.0) * (B + K - 1.0) / (K * (C + K - 1.0)) * Z - ZHF = ZHF + ZR - if (abs(ZHF - ZW) <= abs(ZHF) * EPS): - break - ZW = ZHF - ZHF = Z00 * ZHF - elif (A0 > 1.0): - MAB = np.round(A - B) - if (abs(A - B - MAB) < EPS and A0 <= 1.1): - B = B + EPS - if (abs(A - B - MAB) > EPS): - GA = gamma(A) - GB = gamma(B) - GC = gamma(C) - GAB = gamma(A - B) - GBA = gamma(B - A) - GCA = gamma(C - A) - GCB = gamma(C - B) - ZC0 = GC * GBA / (GCA * GB * (-Z) ** A) - ZC1 = GC * GAB / (GCB * GA * (-Z) ** B) - ZR0 = ZC0 - ZR1 = ZC1 - ZHF = 0.0 + 0j - for K in range(1, 501): - ZR0 = ZR0 * (A + K - 1.0) * (A - C + K) / ((A - B + K) * K * Z) - ZR1 = ZR1 * (B + K - 1.0) * (B - C + K) / ((B - A + K) * K * Z) - ZHF = ZHF + ZR0 + ZR1 - if (abs((ZHF - ZW) / ZHF) <= EPS): - break - ZW = ZHF - ZHF = ZHF + ZC0 + ZC1 - else: - if (A - B < 0.0): - A = BB - B = AA - CA = C - A - CB = C - B - NCA = np.round(CA) - NCB = np.round(CB) - if (abs(CA - NCA) < EPS or abs(CB - NCB) < EPS): - C = C + EPS - GA = gamma(A) - GC = gamma(C) - GCB = gamma(C - B) - PA = psi(A) - PCA = psi(C - A) - PAC = psi(A - C) - MAB = np.round(A - B + EPS) - ZC0 = GC / (GA * (-Z) ** B) - GM = gamma(A - B) - ZF0 = GM / GCB * ZC0 - ZR = ZC0 - for K in range(1, MAB): - ZR = ZR * (B + K - 1.0) / (K * Z) - T0 = A - B - K - G0 = gamma(T0) - GCBK = gamma(C - B - K) - ZF0 = ZF0 + ZR * G0 / GCBK - if (MAB == 0): - ZF0 = 0.0 + 0j - ZC1 = GC / (GA * GCB * (-Z) ** A) - SP = -2.0 * EL - PA - PCA - for J in range(1, MAB + 1): - SP = SP + 1.0 / J - ZP0 = SP + np.log(-Z) - SQ = 1.0 - for J in range(1, MAB + 1): - SQ = SQ * (B + J - 1.0) * (B - C + J) / J - ZF1 = (SQ * ZP0) * ZC1 - ZR = ZC1 - RK1 = 1.0 - SJ1 = 0.0 - W0 = 0.0 - for K in range(1, 10001): - ZR = ZR / Z - RK1 = RK1 * (B + K - 1.0) * (B - C + K) / (K * K) - RK2 = RK1 - for J in range(K + 1, K + MAB + 1): - RK2 = RK2 * (B + J - 1.0) * (B - C + J) / J - SJ1 = SJ1 + \ - (A - 1.0) / (K * (A + K - 1.0)) + \ - (A - C - 1.0) / (K * (A - C + K - 1.0)) - SJ2 = SJ1 - for J in range(K + 1, K + MAB + 1): - SJ2 = SJ2 + 1.0 / J - ZP = -2.0 * EL - PA - PAC + SJ2 - 1.0 / \ - (K + A - C) - PI / np.tan(PI * (K + A - C)) + np.log(-Z) - ZF1 = ZF1 + RK2 * ZR * ZP - WS = abs(ZF1) - if (abs((WS - W0) / WS) < EPS): - break - W0 = WS - ZHF = ZF0 + ZF1 - A = AA - B = BB - if (K > 150): - warnings.warn('Warning! You should check the accuracy') - return ZHF - -# def hypgf(a, b, c, x, abseps=0, releps=1e-13, kmax=10000): -# '''HYPGF Hypergeometric function F(a,b,c,x) -# -# CALL: [y ,abserr] = hypgf(a,b,c,x,abseps,releps) -# -# y = F(a,b,c,x) -# abserr = absolute error estimate -# a,b,c,x = input parameters -# abseps = requested absolute error -# releps = requested relative error -# -# HYPGF calculates one solution to Gauss's hypergeometric differential -# equation: -# -# x*(1-x)Y''(x)+[c-(a+b+1)*x]*Y'(x)-a*b*Y(x) = 0 -# where -# F(a,b,c,x) = Y1(x) = 1 + a*b*x/c + a*(a+1)*b*(b+1)*x^2/(c*(c+1))+.... -# -# -# Many elementary functions are special cases of F(a,b,c,x): -# 1/(1-x) = F(1,1,1,x) = F(1,b,b,x) = F(a,1,a,x) -# (1+x)^n = F(-n,b,b,-x) -# atan(x) = x*F(.5,1,1.5,-x^2) -# asin(x) = x*F(.5,.5,1.5,x^2) -# log(x) = x*F(1,1,2,-x) -# log(1+x)-log(1-x) = 2*x*F(.5,1,1.5,x^2) -# -# NOTE: only real x, abs(x) < 1 and c~=0,-1,-2,... are allowed. -# -# Examples: -# x = linspace(-.99,.99)'; -# [Sn1,err1] = hypgf(1,1,1,x) -# plot(x,abs(Sn1-1./(1-x)),'b',x,err1,'r'),set(gca,'yscale','log') -# [Sn2,err2] = hypgf(.5,.5,1.5,x.^2); -# plot(x,abs(x.*Sn2-asin(x)),'b',x,abs(x.*err2),'r') -# set(gca,'yscale','log') -# -# -# Reference: -# --------- -# Kreyszig, Erwin (1988) -# Advanced engineering mathematics -# John Wiley & Sons, sixth edition, pp 204. -# ''' -# csize = common_shape(x, a, b, c) -# kmin = 2 -# fsum = np.zeros(csize) -# delta = np.zeros(csize) -# err = np.zeros(csize) -# -# ok = ~((np.round(c) == c & c <= 0) | np.abs(x) > 1) -# if np.any(~ok): -# warnings.warn('HYPGF', 'Illegal input: c = 0,-1,-2,... or abs(x)>1') -# fsum[~ok] = np.NaN -# err[~ok] = np.NaN -# -# k0=find(ok & abs(x)==1); -# if any(k0) -# cmab = c(k0)-a(k0)-b(k0); -# fsum(k0) = exp(gammaln(c(k0))+gammaln(cmab)-... -# gammaln(c(k0)-a(k0))-gammaln(c(k0)-b(k0))); -# err(k0) = eps; -# k00 = find(real(cmab)<=0); -# if any(k00) -# err(k0(k00)) = nan; -# fsum(k0(k00)) = nan; -# end -# end -# k=find(ok & abs(x)<1); -# if any(k), -# delta(k) = ones(size(k)); -# fsum(k) = delta(k); -# -# k1 = k; -# E = cell(1,3); -# E{3} = fsum(k); -# converge = 'n'; -# for ix=0:Kmax-1, -# delta(k1) = delta(k1).*((a(k1)+ix)./(ix+1)).*((b(k1)+ix)./(c(k1)+ ix)).*x(k1); @IgnorePep8 -# fsum(k1) = fsum(k1)+delta(k1); -# -# E(1:2) = E(2:3); -# E{3} = fsum(k1); -# -# if ix>Kmin -# if useDEA, -# [Sn, err(k1)] = dea3(E{:}); -# k00 = find((abs(err(k1))) <= max(absEps,abs(relEps.*fsum(k1)))); -# if any(k00) -# fsum(k1(k00)) = Sn(k00); -# end -# if (ix==Kmax-1) -# fsum(k1) = Sn; -# end -# k0 = (find((abs(err(k1))) > max(absEps,abs(relEps.*fsum(k1))))); -# if any(k0),% compute more terms -# %nk=length(k0);%# of values we have to compute again -# E{2} = E{2}(k0); -# E{3} = E{3}(k0); -# else -# converge='y'; -# break; -# end -# else -# err(k1) = 10*abs(delta(k1)); -# k0 = (find((abs(err(k1))) > max(absEps,abs(relEps.* ... -# fsum(k1))))); -# if any(k0),% compute more terms -# %nk=length(k0);%# of values we have to compute again -# else -# converge='y'; -# break; -# end -# end -# k1 = k1(k0); -# end -# end -# if ~strncmpi(converge,'y',1) -# disp(sprintf('#%d values did not converge',length(k1))) -# end -# end -# %ix -# return - - -def nextpow2(x): - ''' - Return next higher power of 2 - - Example - ------- - >>> import wafo.misc as wm - >>> wm.nextpow2(10) - 4 - >>> wm.nextpow2(np.arange(5)) - 3 - ''' - t = isscalar(x) or len(x) - if (t > 1): - f, n = frexp(t) - else: - f, n = frexp(abs(x)) - - if (f == 0.5): - n = n - 1 - return n - - -def discretize(fun, a, b, tol=0.005, n=5, method='linear'): - ''' - Automatic discretization of function - - Parameters - ---------- - fun : callable - function to discretize - a,b : real scalars - evaluation limits - tol : real, scalar - absoute error tolerance - n : scalar integer - number of values - method : string - defining method of gridding, options are 'linear' and 'adaptive' - - Returns - ------- - x : discretized values - y : fun(x) - - Example - ------- - >>> import wafo.misc as wm - >>> import numpy as np - >>> import pylab as plt - >>> x,y = wm.discretize(np.cos, 0, np.pi) - >>> xa,ya = wm.discretize(np.cos, 0, np.pi, method='adaptive') - >>> t = plt.plot(x, y, xa, ya, 'r.') - - plt.show() - >>> plt.close('all') - - ''' - if method.startswith('a'): - return _discretize_adaptive(fun, a, b, tol, n) - else: - return _discretize_linear(fun, a, b, tol, n) - - -def _discretize_linear(fun, a, b, tol=0.005, n=5): - ''' - Automatic discretization of function, linear gridding - ''' - x = linspace(a, b, n) - y = fun(x) - - err0 = inf - err = 10000 - nmax = 2 ** 20 - while (err != err0 and err > tol and n < nmax): - err0 = err - x0 = x - y0 = y - n = 2 * (n - 1) + 1 - x = linspace(a, b, n) - y = fun(x) - y00 = interp(x, x0, y0) - err = 0.5 * amax(abs((y00 - y) / (abs(y00 + y) + _TINY))) - return x, y - - -def _discretize_adaptive(fun, a, b, tol=0.005, n=5): - ''' - Automatic discretization of function, adaptive gridding. - ''' - n += (mod(n, 2) == 0) # make sure n is odd - x = linspace(a, b, n) - fx = fun(x) + Automatic discretization of function, adaptive gridding. + ''' + n += (mod(n, 2) == 0) # make sure n is odd + x = linspace(a, b, n) + fx = fun(x) n2 = (n - 1) / 2 erri = hstack((zeros((n2, 1)), ones((n2, 1)))).ravel() @@ -2860,9 +2504,9 @@ def trangood(x, f, min_n=None, min_x=None, max_x=None, max_n=inf): L = float(xn - x0) if ((nf < min_n) or (max_n < nf) or any(abs(ddx) > 10 * _EPS * (L))): # pab 07.01.2001: Always choose the stepsize df so that - # it is an exactly representable number. - # This is important when calculating numerical derivatives and is - # accomplished by the following. + # it is an exactly representable number. + # This is important when calculating numerical derivatives and is + # accomplished by the following. dx = L / (min(min_n, max_n) - 1) dx = (dx + 2.) - 2. xi = arange(x0, xn + dx / 2., dx) @@ -3103,12 +2747,12 @@ def plot_histgrm(data, bins=None, range=None, # @ReservedAssignment >>> import wafo.misc as wm >>> import wafo.stats as ws >>> R = ws.weibull_min.rvs(2,loc=0,scale=2, size=100) + >>> h0 = wm.plot_histgrm(R, 20, normed=True) >>> x = np.linspace(-3,16,200) >>> h1 = plt.plot(x,ws.weibull_min.pdf(x,2,0,2),'r') >>> plt.close('all') - See also -------- wafo.misc.good_bins @@ -3150,8 +2794,7 @@ def num2pistr(x, n=3): Example >>> import wafo.misc as wm - >>> t = wm.num2pistr(np.pi*3/4) - >>> t=='3\\pi/4' + >>> wm.num2pistr(np.pi*3/4)=='3\\pi/4' True ''' @@ -3181,9 +2824,8 @@ def fourier(data, t=None, T=None, m=None, n=None, method='trapz'): vector or matrix of row vectors with data points shape p x n. t : array-like vector with n values indexed from 1 to N. - T : real scalar + T : real scalar, (default T = t[-1]-t[0]) primitive period of signal, i.e., smallest period. - (default T = t[-1]-t[0] m : scalar integer defines no of harmonics desired (default M = N) n : scalar integer @@ -3191,93 +2833,621 @@ def fourier(data, t=None, T=None, m=None, n=None, method='trapz'): method : string integration method used - Returns - ------- - a,b = Fourier coefficients size m x p + Returns + ------- + a,b = Fourier coefficients size m x p + + FOURIER finds the coefficients for a Fourier series representation + of the signal x(t) (given in digital form). It is assumed the signal + is periodic over T. N is the number of data points, and M-1 is the + number of coefficients. + + The signal can be estimated by using M-1 harmonics by: + M-1 + x[i] = 0.5*a[0] + sum (a[n]*c[n,i] + b[n]*s[n,i]) + n=1 + where + c[n,i] = cos(2*pi*(n-1)*t[i]/T) + s[n,i] = sin(2*pi*(n-1)*t[i]/T) + + Note that a[0] is the "dc value". + Remaining values are a[1], a[2], ... , a[M-1]. + + Example + ------- + >>> import wafo.misc as wm + >>> import numpy as np + >>> T = 2*np.pi + >>> t = np.linspace(0,4*T) + >>> x = np.sin(t) + >>> a, b = wm.fourier(x, t, T=T, m=5) + >>> np.abs(a.ravel())<1e-12 + array([ True, True, True, True, True], dtype=bool) + >>> np.abs(b.ravel()-np.array([ 0., 4., 0., 0., 0.]))<1e-12 + array([ True, True, True, True, True], dtype=bool) + + See also + -------- + fft + ''' + x = np.atleast_2d(data) + p, n = x.shape + if t is None: + t = np.arange(n) + else: + t = np.atleast_1d(t) + + n = len(t) if n is None else n + m = n if n is None else m + T = t[-1] - t[0] if T is None else T + + if method.startswith('trapz'): + intfun = trapz + elif method.startswith('simp'): + intfun = simps + + # Define the vectors for computing the Fourier coefficients + t.shape = (1, -1) + a = zeros((m, p)) + b = zeros((m, p)) + a[0] = intfun(x, t, axis=-1) + + # Compute M-1 more coefficients + tmp = 2 * pi * t / T + # tmp = 2*pi*(0:N-1).'/(N-1); + for i in range(1, m): + a[i] = intfun(x * cos(i * tmp), t, axis=-1) + b[i] = intfun(x * sin(i * tmp), t, axis=-1) + + a = a / pi + b = b / pi + + # Alternative: faster for large M, but gives different results than above. +# nper = diff(t([1 end]))/T; %No of periods given +# if nper == round(nper): +# N1 = n/nper +# else: +# N1 = n +# +# +# +# Fourier coefficients by fft +# Fcof1 = 2*ifft(x(1:N1,:),[],1); +# Pcor = [1; exp(sqrt(-1)*(1:M-1).'*t(1))] # correction term to get +# # the correct integration limits +# Fcof = Fcof1(1:M,:).*Pcor(:,ones(1,P)); +# a = real(Fcof(1:M,:)); +# b = imag(Fcof(1:M,:)); + + return a, b + + +def hyp2f1_taylor(a, b, c, z, tol=1e-13, itermax=500): + a, b, c, z = np.broadcast_arrays(*np.atleast_1d(a, b, c, z)) + shape = a.shape + ak, bk, ck, zk = [d.ravel() for d in (a, b, c, z)] + ajm1 = np.ones(ak.shape) + bjm2 = 0.5 * np.ones(ak.shape) + bjm1 = np.ones(ak.shape) + hout = np.zeros(ak.shape) + k0 = np.arange(len(ak)) + for j in range(0, itermax): + aj = ajm1 * (ak + j) * (bk + j) / (ck + j) * zk / (j + 1) + bj = bjm1 + aj + h, err = dea3(bjm2, bjm1, bj) + k = np.flatnonzero(err > tol * np.abs(h)) + hout[k0] = h + if len(k) == 0: + break + k0 = k0[k] + ak, bk, ck, zk = ak[k], bk[k], ck[k], zk[k] + ajm1 = aj[k] + bjm2 = bjm1[k] + bjm1 = bj[k] + else: + warnings.warn(('Reached %d limit! \n' + + '#%d values did not converge! Max error=%g') % + (j, len(k), np.max(err))) + return hout.reshape(shape) + + +def hyp2f1(a, b, c, z, rho=0.5): + e1 = gammaln(a) + e2 = gammaln(b) + e3 = gammaln(c) + e4 = gammaln(b - a) + e5 = gammaln(a - b) + + e6 = gammaln(c - a) + e7 = gammaln(c - b) + e8 = gammaln(c - a - b) + e9 = gammaln(a + b - c) + _cmab = c - a - b + # ~(np.round(cmab) == cmab & cmab <= 0) + if abs(z) <= rho: + h = hyp2f1_taylor(a, b, c, z, 1e-15) + elif abs(1 - z) <= rho: # % Require that |arg(1-z)| 10: + break + xjm2 = xjm1 + xjm1 = xj + else: + warnings.warn('Reached %d limit' % j) + return h - Note that a[0] is the "dc value". - Remaining values are a[1], a[2], ... , a[M-1]. - Example - ------- - >>> import wafo.misc as wm - >>> import numpy as np - >>> T = 2*np.pi - >>> t = np.linspace(0,4*T) - >>> x = np.sin(t) - >>> a, b = wm.fourier(x, t, T=T, m=5) - >>> np.abs(a.ravel())<1e-12 - array([ True, True, True, True, True], dtype=bool) - >>> np.abs(b.ravel()-np.array([ 0., 4., 0., 0., 0.]))<1e-12 - array([ True, True, True, True, True], dtype=bool) +def hygfz(A, B, C, Z): + ''' Return hypergeometric function for a complex argument, F(a,b,c,z) - See also - -------- - fft + Parameters + ---------- + a, b, c: + parameters where c <> 0,-1,-2,... + z :--- Complex argument ''' - x = np.atleast_2d(data) - p, n = x.shape - if t is None: - t = np.arange(n) - else: - t = np.atleast_1d(t) - - n = len(t) if n is None else n - m = n if n is None else m - T = t[-1] - t[0] if T is None else T - - if method.startswith('trapz'): - intfun = trapz - elif method.startswith('simp'): - intfun = simps + X = np.real(Z) + Y = np.imag(Z) + EPS = 1.0e-15 + L0 = C == np.round(C) and C < 0.0e0 + L1 = abs(1.0 - X) < EPS and Y == 0.0 and C - A - B <= 0.0 + L2 = abs(Z + 1.0) < EPS and abs(C - A + B - 1.0) < EPS + L3 = A == np.round(A) and A < 0.0 + L4 = B == np.round(B) and B < 0.0 + L5 = C - A == np.round(C - A) and C - A <= 0.0 + L6 = C - B == np.round(C - B) and C - B <= 0.0 + AA = A + BB = B + A0 = abs(Z) + if (A0 > 0.95): + EPS = 1.0e-8 + PI = 3.141592653589793 + EL = .5772156649015329 + if (L0 or L1): + # 'The hypergeometric series is divergent' + return np.inf - # Define the vectors for computing the Fourier coefficients - t.shape = (1, -1) - a = zeros((m, p)) - b = zeros((m, p)) - a[0] = intfun(x, t, axis=-1) + NM = 0 + if (A0 == 0.0 or A == 0.0 or B == 0.0): + ZHF = 1.0 + elif (Z == 1.0 and C - A - B > 0.0): + GC = gamma(C) + GCAB = gamma(C - A - B) + GCA = gamma(C - A) + GCB = gamma(C - B) + ZHF = GC * GCAB / (GCA * GCB) + elif L2: + G0 = sqrt(PI) * 2.0 ** (-A) + G1 = gamma(C) + G2 = gamma(1.0 + A / 2.0 - B) + G3 = gamma(0.5 + 0.5 * A) + ZHF = G0 * G1 / (G2 * G3) + elif L3 or L4: + if (L3): + NM = int(np.round(abs(A))) + if (L4): + NM = int(np.round(abs(B))) + ZHF = 1.0 + ZR = 1.0 + for K in range(NM): + ZR = ZR * (A + K) * (B + K) / ((K + 1.) * (C + K)) * Z + ZHF = ZHF + ZR + elif L5 or L6: + if (L5): + NM = np.round(abs(C - A)) + if (L6): + NM = np.round(abs(C - B)) + ZHF = 1.0 + 0j + ZR = 1.0 + 0j + for K in range(NM): + ZR *= (C - A + K) * (C - B + K) / ((K + 1.) * (C + K)) * Z + ZHF = ZHF + ZR + ZHF = (1.0 - Z) ** (C - A - B) * ZHF + elif (A0 <= 1.0): + if (X < 0.0): + Z1 = Z / (Z - 1.0) + if (C > A and B < A and B > 0.0): + A = BB + B = AA - # Compute M-1 more coefficients - tmp = 2 * pi * t / T - # tmp = 2*pi*(0:N-1).'/(N-1); - for i in range(1, m): - a[i] = intfun(x * cos(i * tmp), t, axis=-1) - b[i] = intfun(x * sin(i * tmp), t, axis=-1) + ZC0 = 1.0 / ((1.0 - Z) ** A) + ZHF = 1.0 + 0j + ZR0 = 1.0 + 0j + ZW = 0 + for K in range(500): + ZR0 *= (A + K) * (C - B + K) / ((K + 1.0) * (C + K)) * Z1 + ZHF += ZR0 + if (abs(ZHF - ZW) < abs(ZHF) * EPS): + break + ZW = ZHF + ZHF = ZC0 * ZHF + elif (A0 >= 0.90): + ZW = 0.0 + GM = 0.0 + MCAB = np.round(C - A - B) + if (abs(C - A - B - MCAB) < EPS): + M = int(np.round(C - A - B)) + GA = gamma(A) + GB = gamma(B) + GC = gamma(C) + GAM = gamma(A + M) + GBM = gamma(B + M) + PA = psi(A) + PB = psi(B) + if (M != 0): + GM = 1.0 + for j in range(1, abs(M)): + GM *= j + RM = 1.0 + for j in range(1, abs(M) + 1): # DO 35 J=1,abs(M) + RM *= j + ZF0 = 1.0 + ZR0 = 1.0 + ZR1 = 1.0 + SP0 = 0.0 + SP = 0.0 + if (M >= 0): + ZC0 = GM * GC / (GAM * GBM) + ZC1 = -GC * (Z - 1.0) ** M / (GA * GB * RM) + for K in range(1, M): + ZR0 = ZR0 * \ + (A + K - 1.) * (B + K - 1.) / \ + (K * (K - M)) * (1. - Z) + ZF0 = ZF0 + ZR0 + for K in range(M): + SP0 = SP0 + 1.0 / \ + (A + K) + 1.0 / (B + K) - 1. / (K + 1.) + ZF1 = PA + PB + SP0 + 2.0 * EL + np.log(1.0 - Z) + for K in range(1, 501): + SP = SP + \ + (1.0 - A) / (K * (A + K - 1.0)) + ( + 1.0 - B) / (K * (B + K - 1.0)) + SM = 0.0 + for J in range(1, M): + SM += (1.0 - A) / ( + (J + K) * (A + J + K - 1.0)) + \ + 1.0 / (B + J + K - 1.0) - a = a / pi - b = b / pi + ZP = PA + PB + 2.0 * EL + SP + SM + np.log(1.0 - Z) + ZR1 = ZR1 * \ + (A + M + K - 1.0) * (B + M + K - 1.0) / ( + K * (M + K)) * (1.0 - Z) + ZF1 = ZF1 + ZR1 * ZP + if (abs(ZF1 - ZW) < abs(ZF1) * EPS): + break + ZW = ZF1 + ZHF = ZF0 * ZC0 + ZF1 * ZC1 + elif (M < 0): + M = -M + ZC0 = GM * GC / (GA * GB * (1.0 - Z) ** M) + ZC1 = -(-1) ** M * GC / (GAM * GBM * RM) + for K in range(1, M): + ZR0 = ZR0 * \ + (A - M + K - 1.0) * (B - M + K - 1.0) / ( + K * (K - M)) * (1.0 - Z) + ZF0 = ZF0 + ZR0 + for K in range(1, M + 1): + SP0 = SP0 + 1.0 / K + ZF1 = PA + PB - SP0 + 2.0 * EL + np.log(1.0 - Z) + for K in range(1, 501): + SP = SP + \ + (1.0 - A) / (K * (A + K - 1.0)) + ( + 1.0 - B) / (K * (B + K - 1.0)) + SM = 0.0 + for J in range(1, M + 1): + SM = SM + 1.0 / (J + K) + ZP = PA + PB + 2.0 * EL + SP - SM + np.log(1.0 - Z) + ZR1 = ZR1 * \ + (A + K - 1.) * (B + K - 1.) / \ + (K * (M + K)) * (1. - Z) + ZF1 = ZF1 + ZR1 * ZP + if (abs(ZF1 - ZW) < abs(ZF1) * EPS): + break + ZW = ZF1 + ZHF = ZF0 * ZC0 + ZF1 * ZC1 + else: + GA = gamma(A) + GB = gamma(B) + GC = gamma(C) + GCA = gamma(C - A) + GCB = gamma(C - B) + GCAB = gamma(C - A - B) + GABC = gamma(A + B - C) + ZC0 = GC * GCAB / (GCA * GCB) + ZC1 = GC * GABC / (GA * GB) * (1.0 - Z) ** (C - A - B) + ZHF = 0 + 0j + ZR0 = ZC0 + ZR1 = ZC1 + for K in range(1, 501): + ZR0 = ZR0 * \ + (A + K - 1.) * (B + K - 1.) / \ + (K * (A + B - C + K)) * (1. - Z) + ZR1 = ZR1 * \ + (C - A + K - 1.0) * (C - B + K - 1.0) / ( + K * (C - A - B + K)) * (1.0 - Z) + ZHF = ZHF + ZR0 + ZR1 + if (abs(ZHF - ZW) < abs(ZHF) * EPS): + break + ZW = ZHF + ZHF = ZHF + ZC0 + ZC1 + else: + ZW = 0.0 + Z00 = 1.0 # + 0j + if (C - A < A and C - B < B): + Z00 = (1.0 - Z) ** (C - A - B) + A = C - A + B = C - B + ZHF = 1.0 + ZR = 1.0 + for K in range(1, 501): + ZR = ZR * \ + (A + K - 1.0) * (B + K - 1.0) / (K * (C + K - 1.0)) * Z + ZHF = ZHF + ZR + if (abs(ZHF - ZW) <= abs(ZHF) * EPS): + break + ZW = ZHF + ZHF = Z00 * ZHF + elif (A0 > 1.0): + MAB = np.round(A - B) + if (abs(A - B - MAB) < EPS and A0 <= 1.1): + B = B + EPS + if (abs(A - B - MAB) > EPS): + GA = gamma(A) + GB = gamma(B) + GC = gamma(C) + GAB = gamma(A - B) + GBA = gamma(B - A) + GCA = gamma(C - A) + GCB = gamma(C - B) + ZC0 = GC * GBA / (GCA * GB * (-Z) ** A) + ZC1 = GC * GAB / (GCB * GA * (-Z) ** B) + ZR0 = ZC0 + ZR1 = ZC1 + ZHF = 0.0 + 0j + for K in range(1, 501): + ZR0 = ZR0 * (A + K - 1.0) * (A - C + K) / ((A - B + K) * K * Z) + ZR1 = ZR1 * (B + K - 1.0) * (B - C + K) / ((B - A + K) * K * Z) + ZHF = ZHF + ZR0 + ZR1 + if (abs((ZHF - ZW) / ZHF) <= EPS): + break + ZW = ZHF + ZHF = ZHF + ZC0 + ZC1 + else: + if (A - B < 0.0): + A = BB + B = AA + CA = C - A + CB = C - B + NCA = np.round(CA) + NCB = np.round(CB) + if (abs(CA - NCA) < EPS or abs(CB - NCB) < EPS): + C = C + EPS + GA = gamma(A) + GC = gamma(C) + GCB = gamma(C - B) + PA = psi(A) + PCA = psi(C - A) + PAC = psi(A - C) + MAB = np.round(A - B + EPS) + ZC0 = GC / (GA * (-Z) ** B) + GM = gamma(A - B) + ZF0 = GM / GCB * ZC0 + ZR = ZC0 + for K in range(1, MAB): + ZR = ZR * (B + K - 1.0) / (K * Z) + T0 = A - B - K + G0 = gamma(T0) + GCBK = gamma(C - B - K) + ZF0 = ZF0 + ZR * G0 / GCBK + if (MAB == 0): + ZF0 = 0.0 + 0j + ZC1 = GC / (GA * GCB * (-Z) ** A) + SP = -2.0 * EL - PA - PCA + for J in range(1, MAB + 1): + SP = SP + 1.0 / J + ZP0 = SP + np.log(-Z) + SQ = 1.0 + for J in range(1, MAB + 1): + SQ = SQ * (B + J - 1.0) * (B - C + J) / J + ZF1 = (SQ * ZP0) * ZC1 + ZR = ZC1 + RK1 = 1.0 + SJ1 = 0.0 + W0 = 0.0 + for K in range(1, 10001): + ZR = ZR / Z + RK1 = RK1 * (B + K - 1.0) * (B - C + K) / (K * K) + RK2 = RK1 + for J in range(K + 1, K + MAB + 1): + RK2 = RK2 * (B + J - 1.0) * (B - C + J) / J + SJ1 = SJ1 + \ + (A - 1.0) / (K * (A + K - 1.0)) + \ + (A - C - 1.0) / (K * (A - C + K - 1.0)) + SJ2 = SJ1 + for J in range(K + 1, K + MAB + 1): + SJ2 = SJ2 + 1.0 / J + ZP = -2.0 * EL - PA - PAC + SJ2 - 1.0 / \ + (K + A - C) - PI / np.tan(PI * (K + A - C)) + np.log(-Z) + ZF1 = ZF1 + RK2 * ZR * ZP + WS = abs(ZF1) + if (abs((WS - W0) / WS) < EPS): + break + W0 = WS + ZHF = ZF0 + ZF1 + A = AA + B = BB + if (K > 150): + warnings.warn('Warning! You should check the accuracy') + return ZHF - # Alternative: faster for large M, but gives different results than above. -# nper = diff(t([1 end]))/T; %No of periods given -# if nper == round(nper): -# N1 = n/nper -# else: -# N1 = n +# def hypgf(a, b, c, x, abseps=0, releps=1e-13, kmax=10000): +# '''HYPGF Hypergeometric function F(a,b,c,x) # +# CALL: [y ,abserr] = hypgf(a,b,c,x,abseps,releps) # +# y = F(a,b,c,x) +# abserr = absolute error estimate +# a,b,c,x = input parameters +# abseps = requested absolute error +# releps = requested relative error # -# Fourier coefficients by fft -# Fcof1 = 2*ifft(x(1:N1,:),[],1); -# Pcor = [1; exp(sqrt(-1)*(1:M-1).'*t(1))]; % correction term to get -# % the correct integration limits -# Fcof = Fcof1(1:M,:).*Pcor(:,ones(1,P)); -# a = real(Fcof(1:M,:)); -# b = imag(Fcof(1:M,:)); +# HYPGF calculates one solution to Gauss's hypergeometric differential +# equation: +# +# x*(1-x)Y''(x)+[c-(a+b+1)*x]*Y'(x)-a*b*Y(x) = 0 +# where +# F(a,b,c,x) = Y1(x) = 1 + a*b*x/c + a*(a+1)*b*(b+1)*x^2/(c*(c+1))+.... +# +# +# Many elementary functions are special cases of F(a,b,c,x): +# 1/(1-x) = F(1,1,1,x) = F(1,b,b,x) = F(a,1,a,x) +# (1+x)^n = F(-n,b,b,-x) +# atan(x) = x*F(.5,1,1.5,-x^2) +# asin(x) = x*F(.5,.5,1.5,x^2) +# log(x) = x*F(1,1,2,-x) +# log(1+x)-log(1-x) = 2*x*F(.5,1,1.5,x^2) +# +# NOTE: only real x, abs(x) < 1 and c~=0,-1,-2,... are allowed. +# +# Examples: +# x = linspace(-.99,.99)'; +# [Sn1,err1] = hypgf(1,1,1,x) +# plot(x,abs(Sn1-1./(1-x)),'b',x,err1,'r'),set(gca,'yscale','log') +# [Sn2,err2] = hypgf(.5,.5,1.5,x.^2); +# plot(x,abs(x.*Sn2-asin(x)),'b',x,abs(x.*err2),'r') +# set(gca,'yscale','log') +# +# +# Reference: +# --------- +# Kreyszig, Erwin (1988) +# Advanced engineering mathematics +# John Wiley & Sons, sixth edition, pp 204. +# ''' +# csize = common_shape(x, a, b, c) +# kmin = 2 +# fsum = np.zeros(csize) +# delta = np.zeros(csize) +# err = np.zeros(csize) +# +# ok = ~((np.round(c) == c & c <= 0) | np.abs(x) > 1) +# if np.any(~ok): +# warnings.warn('HYPGF', 'Illegal input: c = 0,-1,-2,... or abs(x)>1') +# fsum[~ok] = np.NaN +# err[~ok] = np.NaN +# +# k0=find(ok & abs(x)==1); +# if any(k0) +# cmab = c(k0)-a(k0)-b(k0); +# fsum(k0) = exp(gammaln(c(k0))+gammaln(cmab)-... +# gammaln(c(k0)-a(k0))-gammaln(c(k0)-b(k0))); +# err(k0) = eps; +# k00 = find(real(cmab)<=0); +# if any(k00) +# err(k0(k00)) = nan; +# fsum(k0(k00)) = nan; +# end +# end +# k=find(ok & abs(x)<1); +# if any(k), +# delta(k) = ones(size(k)); +# fsum(k) = delta(k); +# +# k1 = k; +# E = cell(1,3); +# E{3} = fsum(k); +# converge = 'n'; +# for ix=0:Kmax-1, +# delta(k1) = delta(k1).*((a(k1)+ix)./(ix+1)).*((b(k1)+ix)./(c(k1)+ ix)).*x(k1); @IgnorePep8 +# fsum(k1) = fsum(k1)+delta(k1); +# +# E(1:2) = E(2:3); +# E{3} = fsum(k1); +# +# if ix>Kmin +# if useDEA, +# [Sn, err(k1)] = dea3(E{:}); +# k00 = find((abs(err(k1))) <= max(absEps,abs(relEps.*fsum(k1)))); +# if any(k00) +# fsum(k1(k00)) = Sn(k00); +# end +# if (ix==Kmax-1) +# fsum(k1) = Sn; +# end +# k0 = (find((abs(err(k1))) > max(absEps,abs(relEps.*fsum(k1))))); +# if any(k0),% compute more terms +# %nk=length(k0);%# of values we have to compute again +# E{2} = E{2}(k0); +# E{3} = E{3}(k0); +# else +# converge='y'; +# break; +# end +# else +# err(k1) = 10*abs(delta(k1)); +# k0 = (find((abs(err(k1))) > max(absEps,abs(relEps.* ... +# fsum(k1))))); +# if any(k0),% compute more terms +# %nk=length(k0);%# of values we have to compute again +# else +# converge='y'; +# break; +# end +# end +# k1 = k1(k0); +# end +# end +# if ~strncmpi(converge,'y',1) +# disp(sprintf('#%d values did not converge',length(k1))) +# end +# end +# %ix +# return - return a, b + +def _test_find_cross(): + t = findcross([0, 0, 1, -1, 1], 0) # @UnusedVariable def real_main0(): diff --git a/wafo/numpy_utils.py b/wafo/numpy_utils.py index 2bcd2c4..206fcee 100644 --- a/wafo/numpy_utils.py +++ b/wafo/numpy_utils.py @@ -6,12 +6,13 @@ import collections import sys import fractions import numpy as np -from numpy import (abs, amax, any, logical_and, arange, linspace, atleast_1d, - array, asarray, broadcast_arrays, ceil, floor, frexp, hypot, - sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff, finfo, - inf, pi, interp, isnan, isscalar, zeros, ones, linalg, - r_, sign, unique, hstack, vstack, nonzero, where, extract, - meshgrid) +from numpy import ( + meshgrid, + abs, amax, any, logical_and, arange, linspace, atleast_1d, + asarray, ceil, floor, frexp, hypot, + sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff, + finfo, inf, pi, interp, isnan, isscalar, zeros, ones, linalg, + r_, sign, unique, hstack, vstack, nonzero, where, extract) from scipy.special import gammaln from scipy.integrate import trapz, simps import warnings @@ -22,6 +23,8 @@ try: except: clib = None floatinfo = finfo(float) +_TINY = np.finfo(float).tiny +_EPS = np.finfo(float).eps __all__ = ['now', 'spaceline', 'narg_smallest', 'args_flat', 'is_numlike', @@ -29,9 +32,10 @@ __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', - 'betaloge', 'gravity', 'nextpow2', 'discretize', 'pol2cart', - 'cart2pol', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', + 'betaloge', 'gravity', 'nextpow2', 'discretize', 'polar2cart', + 'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere', + 'piecewise', 'valarray'] @@ -209,23 +213,24 @@ def rotation_matrix(heading, pitch, roll): ''' Examples + >>> import numpy as np >>> rotation_matrix(heading=0, pitch=0, roll=0) array([[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 1.]]) >>> np.all(np.abs(rotation_matrix(heading=180, pitch=0, roll=0)- - ... array([[ -1.000000e+00, -1.224647e-16, 0.000000e+00], + ... np.array([[ -1.000000e+00, -1.224647e-16, 0.000000e+00], ... [ 1.224647e-16, -1.000000e+00, 0.000000e+00], ... [ -0.000000e+00, 0.000000e+00, 1.000000e+00]]))<1e-7) True >>> np.all(np.abs(rotation_matrix(heading=0, pitch=180, roll=0)- - ... array([[ -1.000000e+00, 0.000000e+00, 1.224647e-16], + ... np.array([[ -1.000000e+00, 0.000000e+00, 1.224647e-16], ... [ -0.000000e+00, 1.000000e+00, 0.000000e+00], ... [ -1.224647e-16, -0.000000e+00, -1.000000e+00]]))<1e-7) True >>> np.all(np.abs(rotation_matrix(heading=0, pitch=0, roll=180)- - ... array([[ 1.000000e+00, 0.000000e+00, 0.000000e+00], + ... np.array([[ 1.000000e+00, 0.000000e+00, 0.000000e+00], ... [ 0.000000e+00, -1.000000e+00, -1.224647e-16], ... [ -0.000000e+00, 1.224647e-16, -1.000000e+00]]))<1e-7) True @@ -235,7 +240,7 @@ def rotation_matrix(heading, pitch, roll): deg2rad = np.pi / 180 H = heading * deg2rad P = pitch * deg2rad - R = roll * deg2rad # Convert to radians + R = roll * deg2rad # Convert to radians data.put(0, cos(H) * cos(P)) data.put(1, cos(H) * sin(P) * sin(R) - sin(H) * cos(R)) @@ -327,13 +332,12 @@ def spaceline(start_point, stop_point, num=10): [ 2.5 , 0. , 0. ], [ 2.75, 0. , 0. ], [ 3. , 0. , 0. ]]) - ''' num = int(num) e1, e2 = np.atleast_1d(start_point, stop_point) e2m1 = e2 - e1 length = np.sqrt((e2m1 ** 2).sum()) - # length = sqrt((E2[0]-E1(1))^2 + (E2(2)-E1(2))^2 + (E2(3)-E1(3))^2); + # length = sqrt((E2[0]-E1(1))^2 + (E2(2)-E1(2))^2 + (E2(3)-E1(3))^2) C = e2m1 / length delta = length / float(num - 1) return np.array([e1 + n * delta * C for n in range(num)]) @@ -400,37 +404,14 @@ def args_flat(*args): 'POS array must be of shape N x 3!') return pos, None elif nargin == 3: - x, y, z = broadcast_arrays(*args[:3]) + x, y, z = np.broadcast_arrays(*args[:3]) c_shape = x.shape return np.vstack((x.ravel(), y.ravel(), z.ravel())).T, c_shape else: raise ValueError('Number of arguments must be 1 or 3!') -def _check_and_adjust_shape(shape, nsub=None): - s = np.atleast_1d(shape) - ndim = len(s) - if ndim < 1: - raise ValueError('Shape vector must have at least 1 element.') - ndim = len(s) - if nsub is None: - nsub = ndim - if ndim <= nsub: # add trailing singleton dimensions - s = np.hstack([s, np.ones(nsub - ndim, dtype=int)]) - else: # Adjust for linear indexing on last element - s = np.hstack([s[:nsub - 1], np.prod(s[nsub - 1:])]) - return s - - -def _sub2index_factor(shape, order='C'): - ''' Return multiplier needed for calculating linear index - from multiple subscripts. - ''' - step = 1 if order == 'F' else -1 # C order - return np.hstack([1, np.cumprod(shape[::step][:-1])])[::step] - - -def index2sub(shape, index, nsub=None, order='C'): +def index2sub(shape, index, order='C'): ''' Returns Multiple subscripts from linear index. @@ -440,8 +421,6 @@ def index2sub(shape, index, nsub=None, order='C'): shape of array index : linear index into array - nsub : int optional - Number of subscripts returned. default nsub=len(shape) order : {'C','F'}, optional The order of the linear index. 'C' means C (row-major) order. @@ -468,18 +447,7 @@ def index2sub(shape, index, nsub=None, order='C'): -------- sub2index ''' - ndx = np.atleast_1d(index) - s = _check_and_adjust_shape(shape, nsub) - k = _sub2index_factor(s, order) - n = len(s) - step = -1 if order == 'F' else 1 # C order - subscripts = [0, ] * n - for i in range(n)[::step]: - vi = np.remainder(ndx, k[i]) - subscript = np.array((ndx - vi) / k[i], dtype=int) - subscripts[i] = subscript - ndx = vi - return tuple(subscripts) + return np.unravel_index(index, shape, order=order) def sub2index(shape, *subscripts, **kwds): @@ -518,21 +486,7 @@ def sub2index(shape, *subscripts, **kwds): -------- index2sub ''' - nsub = len(subscripts) - s = _check_and_adjust_shape(shape, nsub) - k = _sub2index_factor(s, **kwds) - - ndx = 0 - s0 = np.shape(subscripts[0]) - for i, subscript in enumerate(subscripts): - np.testing.assert_equal(s0, np.shape(subscript), - 'The subscripts vectors must all ' + - 'be of the same shape.') - if (np.any(subscript < 0)) or (np.any(s[i] <= subscript)): - raise IndexError('Out of range subscript.') - - ndx = ndx + k[i] * subscript - return ndx + return np.ravel_multi_index(subscripts, shape, **kwds) def is_numlike(obj): @@ -683,12 +637,14 @@ def detrendma(x, L): Examples -------- + >>> import utilities.numpy_utils as wm >>> import pylab as plt >>> exp = plt.exp; cos = plt.cos; randn = plt.randn >>> x = plt.linspace(0,1,200) >>> y = exp(x)+cos(5*2*pi*x)+1e-1*randn(x.size) - >>> y0 = detrendma(y,20); tr = y-y0 + >>> y0 = wm.detrendma(y,20); tr = y-y0 >>> h = plt.plot(x, y, x, y0, 'r', x, exp(x), 'k', x, tr, 'm') + >>> plt.close('all') See also @@ -744,7 +700,7 @@ def ecross(t, f, ind, v=0): Example ------- >>> from matplotlib import pylab as plt - >>> import utilities.numpy_utils as wm + >>> import wafo.misc as wm >>> ones = np.ones >>> t = np.linspace(0,7*np.pi,250) >>> x = np.sin(t) @@ -752,8 +708,8 @@ def ecross(t, f, ind, v=0): >>> ind array([ 9, 25, 80, 97, 151, 168, 223, 239]) >>> t0 = wm.ecross(t,x,ind,0.75) - >>> np.abs(t0 - np.array([ 0.84910514, 2.2933879, 7.13205663, 8.57630119, - ... 13.41484739, 14.85909194, 19.69776067, 21.14204343]))<1e-7 + >>> np.abs(t0 - np.array([0.84910514, 2.2933879 , 7.13205663, 8.57630119, + ... 13.41484739, 14.85909194, 19.69776067, 21.14204343]))<1e-7 array([ True, True, True, True, True, True, True, True], dtype=bool) >>> a = plt.plot(t, x, '.', t[ind], x[ind], 'r.', t, ones(t.shape)*0.75, @@ -767,8 +723,8 @@ def ecross(t, f, ind, v=0): # Tested on: Python 2.5 # revised pab Feb2004 # By pab 18.06.2001 - return t[ind] + (v - f[ind]) * (t[ind + 1] - t[ind]) / \ - (f[ind + 1] - f[ind]) + return (t[ind] + (v - f[ind]) * (t[ind + 1] - t[ind]) / + (f[ind + 1] - f[ind])) def _findcross(xn): @@ -806,6 +762,13 @@ def _findcross(xn): return ind +def xor(a, b): + """ + Return True only when inputs differ. + """ + return a ^ b + + def findcross(x, v=0.0, kind=None): ''' Return indices to level v up and/or downcrossings of a vector @@ -872,20 +835,19 @@ def findcross(x, v=0.0, kind=None): t_0 = int(xn[ind[0] + 1] < 0) ind = ind[t_0::2] elif kind in ('dw', 'uw', 'tw', 'cw'): - # make sure that the first is a level v down-crossing - # if kind=='dw' or kind=='tw' - # or that the first is a level v up-crossing - # if kind=='uw' or kind=='cw' - xor = lambda a, b: a ^ b + # make sure the first is a level v down-crossing + # if wdef=='dw' or wdef=='tw' + # or make sure the first is a level v up-crossing + # if wdef=='uw' or wdef=='cw' + first_is_down_crossing = int(xn[ind[0]] > xn[ind[0] + 1]) if xor(first_is_down_crossing, kind in ('dw', 'tw')): ind = ind[1::] - n_c = ind.size # number of level v crossings # make sure the number of troughs and crests are according to the # wavedef, i.e., make sure length(ind) is odd if dw or uw # and even if tw or cw - is_odd = mod(n_c, 2) + is_odd = mod(ind.size, 2) if xor(is_odd, kind in ('dw', 'uw')): ind = ind[:-1] else: @@ -1041,7 +1003,7 @@ def findrfc_astm(tp): return sig_rfc -def findrfc(tp, hmin=0.0, method='clib'): +def findrfc(tp, h=0.0, method='clib'): ''' Return indices to rainflow cycles of a sequence of TP. @@ -1064,10 +1026,10 @@ def findrfc(tp, hmin=0.0, method='clib'): Example: -------- - >>> import pylab as plt + >>> import matplotlib.pyplot as plt >>> import utilities.numpy_utils as wm - >>> t = plt.linspace(0,7*np.pi,250) - >>> x = plt.sin(t)+0.1*np.sin(50*t) + >>> t = np.linspace(0,7*np.pi,250) + >>> x = np.sin(t)+0.1*np.sin(50*t) >>> ind = wm.findextrema(x) >>> ti, tp = t[ind], x[ind] @@ -1127,7 +1089,7 @@ def findrfc(tp, hmin=0.0, method='clib'): Tmi = Tstart + 2 * j j -= 1 if (xminus >= xplus): - if (y[2 * i + 1] - xminus >= hmin): + if (y[2 * i + 1] - xminus >= h): ind[ix] = Tmi ix += 1 ind[ix] = (Tstart + 2 * i + 1) @@ -1143,7 +1105,7 @@ def findrfc(tp, hmin=0.0, method='clib'): Tpl = (Tstart + 2 * j + 2) j += 1 else: - if ((y[2 * i + 1] - xminus) >= hmin): + if ((y[2 * i + 1] - xminus) >= h): ind[ix] = Tmi ix += 1 ind[ix] = (Tstart + 2 * i + 1) @@ -1154,12 +1116,12 @@ def findrfc(tp, hmin=0.0, method='clib'): # goto L180 # L170: if (xplus <= xminus): - if ((y[2 * i + 1] - xminus) >= hmin): + if ((y[2 * i + 1] - xminus) >= h): ind[ix] = Tmi ix += 1 ind[ix] = (Tstart + 2 * i + 1) ix += 1 - elif ((y[2 * i + 1] - xplus) >= hmin): + elif ((y[2 * i + 1] - xplus) >= h): ind[ix] = (Tstart + 2 * i + 1) ix += 1 ind[ix] = Tpl @@ -1169,7 +1131,7 @@ def findrfc(tp, hmin=0.0, method='clib'): # iy=i # /* for i */ else: - ind, ix = clib.findrfc(y, hmin) + ind, ix = clib.findrfc(y, h) return np.sort(ind[:ix]) @@ -1271,13 +1233,8 @@ def mctp2rfc(fmM, fMm=None): fx = NN * (A / (1 - B * A) * e) else: rh = np.eye(A.shape[0]) - np.dot(B, A) - fx = np.dot( - NN, - np.dot( - A, - linalg.solve( - rh, - e))) # least squares + # least squares + fx = np.dot(NN, np.dot(A, linalg.solve(rh, e))) # end # end f_rfc[N - 1 - k, k - i] = fx + DRFC @@ -1339,19 +1296,29 @@ def rfcfilter(x, h, method=0): Examples: --------- # 1. Filtered signal y is the turning points of x. - >>> import utilities.data + >>> import utilities.data as data >>> import utilities.numpy_utils as wm - >>> x = utilities.data.sea() + >>> x = data.sea() >>> y = wm.rfcfilter(x[:,1], h=0, method=1) >>> np.all(np.abs(y[0:5]-np.array([-1.2004945 , 0.83950546, -0.09049454, ... -0.02049454, -0.09049454]))<1e-7) True + >>> y.shape + (2172,) # 2. This removes all rainflow cycles with range less than 0.5. >>> y1 = wm.rfcfilter(x[:,1], h=0.5) + >>> y1.shape + (863,) >>> np.all(np.abs(y1[0:5]-np.array([-1.2004945 , 0.83950546, -0.43049454, - ... 0.34950546, -0.51049454]))<1e-7) + ... 0.34950546, -0.51049454]))<1e-7) True + >>> ind = wm.findtp(x[:,1], h=0.5) + >>> y2 = x[ind,1] + >>> y2[0:5] + array([-1.2004945 , 0.83950546, -0.43049454, 0.34950546, -0.51049454]) + >>> y2[-5::] + array([ 0.83950546, -0.64049454, 0.65950546, -1.0004945 , 0.91950546]) See also -------- @@ -1364,21 +1331,27 @@ def rfcfilter(x, h, method=0): j = 0 t0 = 0 y0 = y[t0] - z0 = 0 + + def aleb(a, b): + return a <= b + + def altb(a, b): + return a < b + if method == 0: - cmpfun1 = lambda a, b: a <= b - cmpfun2 = lambda a, b: a < b + cmpfun1 = aleb + cmpfun2 = altb else: - cmpfun1 = lambda a, b: a < b - cmpfun2 = lambda a, b: a <= b + cmpfun1 = altb + cmpfun2 = aleb # The rainflow filter for tim1, yi in enumerate(y[1::]): fpi = y0 + h fmi = y0 - h ti = tim1 + 1 - #yi = y[ti] + # yi = y[ti] if z0 == 0: if cmpfun1(yi, fmi): @@ -1390,7 +1363,7 @@ def rfcfilter(x, h, method=0): t1, y1 = (t0, y0) if z1 == 0 else (ti, yi) else: if (((z0 == +1) & cmpfun1(yi, fmi)) | - ((z0 == -1) & cmpfun2(yi, fpi))): + ((z0 == -1) & cmpfun2(yi, fpi))): z1 = -1 elif (((z0 == +1) & cmpfun2(fmi, yi)) | ((z0 == -1) & cmpfun1(fpi, yi))): @@ -1417,7 +1390,7 @@ def rfcfilter(x, h, method=0): t0, y0, z0 = t1, y1, z1 # end - #% Update y if last y0 is greater than (or equal) threshold + # Update y if last y0 is greater than (or equal) threshold if cmpfun1(h, abs(y0 - y[t[j]])): j += 1 t[j] = t0 @@ -1456,19 +1429,20 @@ def findtp(x, h=0.0, kind=None): >>> import pylab as plt >>> import utilities.numpy_utils as wm >>> t = np.linspace(0,30,500).reshape((-1,1)) - >>> x = np.hstack((t, np.cos(t))) - >>> x1 = x[0:200,:] + >>> x = np.hstack((t, np.cos(t) + 0.3 * np.sin(5*t))) + >>> x1 = x[0:100,:] >>> itp = wm.findtp(x1[:,1],0,'Mw') >>> itph = wm.findtp(x1[:,1],0.3,'Mw') >>> tp = x1[itp,:] >>> tph = x1[itph,:] - >>> a = plt.plot(x1[:,0],x1[:,1],tp[:,0],tp[:,1],'ro', - ... tph[:,1],tph[:,1],'k.') + >>> a = plt.plot(x1[:,0],x1[:,1], + ... tp[:,0],tp[:,1],'ro', + ... tph[:,0],tph[:,1],'k.') >>> plt.close('all') >>> itp - array([105, 157, 199]) + array([ 5, 18, 24, 38, 46, 57, 70, 76, 91, 98, 99]) >>> itph - array([105]) + array([91]) See also --------- @@ -1487,9 +1461,9 @@ def findtp(x, h=0.0, kind=None): return None # In order to get the exact up-crossing intensity from rfc by - # mm2lc(tp2mm(rfc)) we have to add the indices - # to the last value (and also the first if the - # sequence of turning points does not start with a minimum). + # mm2lc(tp2mm(rfc)) we have to add the indices to the last value + # (and also the first if the sequence of turning points does not start + # with a minimum). if kind == 'astm': # the Nieslony approach always put the first loading point as the first @@ -1510,7 +1484,6 @@ def findtp(x, h=0.0, kind=None): ind = ind[ind1] if kind in ('mw', 'Mw'): - xor = lambda a, b: a ^ b # make sure that the first is a Max if wdef == 'Mw' # or make sure that the first is a min if wdef == 'mw' first_is_max = (x[ind[0]] > x[ind[1]]) @@ -1582,11 +1555,8 @@ def findtc(x_in, v=None, kind=None): return zeros(0, dtype=np.int), zeros(0, dtype=np.int) # determine the number of trough2crest (or crest2trough) cycles - isodd = mod(n_c, 2) - if isodd: - n_tc = int((n_c - 1) / 2) - else: - n_tc = int((n_c - 2) / 2) + is_even = mod(n_c + 1, 2) + n_tc = int((n_c - 1 - is_even) / 2) # allocate variables before the loop increases the speed ind = zeros(n_c - 1, dtype=np.int) @@ -1604,16 +1574,16 @@ def findtc(x_in, v=None, kind=None): # trough ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]].argmin() - else: # %%%% the first is a up-crossing + else: # the first is a up-crossing for i in xrange(n_tc): - # trough + # crest j = 2 * i ind[j] = x[v_ind[j] + 1:v_ind[j + 1] + 1].argmax() - # crest + # trough ind[j + 1] = x[v_ind[j + 1] + 1:v_ind[j + 2] + 1].argmin() if (2 * n_tc + 1 < n_c) and (kind in (None, 'cw')): - # trough + # crest ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]].argmax() return v_ind[:n_c - 1] + ind + 1, v_ind @@ -1786,9 +1756,8 @@ def findoutliers(x, zcrit=0.0, dcrit=None, ddcrit=None, verbose=False): if zcrit == 0.: print('Found %d consecutive equal values' % indz.size) else: - print( - 'Found %d consecutive values less than %g apart.' % - (indz.size, zcrit)) + print('Found %d consecutive values less than %g apart.' % + (indz.size, zcrit)) indg = ones(xn.size, dtype=bool) if ind.size > 1: @@ -1937,7 +1906,7 @@ def stirlerr(n): Example ------- >>> import utilities.numpy_utils as wm - >>> np.abs(wm.stirlerr(2)-array([ 0.0413407]))<1e-7 + >>> np.abs(wm.stirlerr(2)- 0.0413407)<1e-7 array([ True], dtype=bool) See also @@ -1949,11 +1918,10 @@ def stirlerr(n): ----------- Catherine Loader (2000). Fast and Accurate Computation of Binomial Probabilities - - + ''' - S0 = 0.083333333333333333333 # /* 1/12 */ + S0 = 0.083333333333333333333 # /* 1/12 */ S1 = 0.00277777777777777777778 # /* 1/360 */ S2 = 0.00079365079365079365079365 # /* 1/1260 */ S3 = 0.000595238095238095238095238 # /* 1/1680 */ @@ -2016,9 +1984,9 @@ def binomln(z, w): # log(n!) = stirlerr(n) + log( sqrt(2*pi*n)*(n/exp(1))**n ) # y = gammaln(z+1)-gammaln(w+1)-gammaln(z-w+1) zpw = z - w - return (stirlerr(z + 1) - stirlerr(w + 1) - 0.5 * log(2 * pi) - - (w + 0.5) * log1p(w) + (z + 0.5) * log1p(z) - stirlerr(zpw + 1) - - (zpw + 0.5) * log1p(zpw) + 1) + return (stirlerr(z + 1) - stirlerr(w + 1) - 0.5 * log(2 * pi) - + (w + 0.5) * log1p(w) + (z + 0.5) * log1p(z) - stirlerr(zpw + 1) - + (zpw + 0.5) * log1p(zpw) + 1) def betaloge(z, w): @@ -2052,8 +2020,9 @@ def betaloge(z, w): ''' # y = gammaln(z)+gammaln(w)-gammaln(z+w) zpw = z + w - return (stirlerr(z) + stirlerr(w) + 0.5 * log(2 * pi) + (w - 0.5) * log(w) - + (z - 0.5) * log(z) - stirlerr(zpw) - (zpw - 0.5) * log(zpw)) + return (stirlerr(z) + stirlerr(w) + 0.5 * log(2 * pi) + + (w - 0.5) * log(w) + (z - 0.5) * log(z) - stirlerr(zpw) - + (zpw - 0.5) * log(zpw)) # stirlings approximation: # (-(zpw-0.5).*log(zpw) +(w-0.5).*log(w)+(z-0.5).*log(z) +0.5*log(2*pi)) @@ -2085,7 +2054,7 @@ def gravity(phi=45): >>> import utilities.numpy_utils as wm >>> import numpy as np >>> phi = np.linspace(0,45,5) - >>> np.abs(wm.gravity(phi)-array([ 9.78049 , 9.78245014, 9.78803583, + >>> np.abs(wm.gravity(phi)-np.array([ 9.78049 , 9.78245014, 9.78803583, ... 9.79640552, 9.80629387]))<1.e-7 array([ True, True, True, True, True], dtype=bool) @@ -2104,8 +2073,81 @@ def gravity(phi=45): ''' phir = phi * pi / 180. # change from degrees to radians - return 9.78049 * \ - (1. + 0.0052884 * sin(phir) ** 2. - 0.0000059 * sin(2 * phir) ** 2.) + return 9.78049 * (1. + 0.0052884 * sin(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): @@ -2176,8 +2218,6 @@ def _discretize_linear(fun, a, b, tol=0.005, n=5): ''' Automatic discretization of function, linear gridding ''' - tiny = floatinfo.tiny - x = linspace(a, b, n) y = fun(x) @@ -2192,7 +2232,7 @@ def _discretize_linear(fun, a, b, tol=0.005, n=5): x = linspace(a, b, n) y = fun(x) y00 = interp(x, x0, y0) - err = 0.5 * amax(abs((y00 - y) / (abs(y00 + y) + tiny))) + err = 0.5 * amax(abs((y00 - y) / (abs(y00 + y) + _TINY))) return x, y @@ -2200,7 +2240,6 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5): ''' Automatic discretization of function, adaptive gridding. ''' - tiny = floatinfo.tiny n += (mod(n, 2) == 0) # make sure n is odd x = linspace(a, b, n) fx = fun(x) @@ -2222,7 +2261,7 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5): fy = fun(y) fy0 = interp(y, x, fx) - erri = 0.5 * (abs((fy0 - fy) / (abs(fy0 + fy) + tiny))) + erri = 0.5 * (abs((fy0 - fy) / (abs(fy0 + fy) + _TINY))) err = erri.max() @@ -2281,7 +2320,6 @@ def cart2polar(x, y, z=None): return t, r else: return t, r, z - cart2pol = cart2polar @@ -2358,8 +2396,7 @@ def trangood(x, f, min_n=None, min_x=None, max_x=None, max_n=inf): xn = xo[-1] x0 = xo[0] L = float(xn - x0) - eps = floatinfo.eps - if ((nf < min_n) or (max_n < nf) or any(abs(ddx) > 10 * eps * (L))): + if ((nf < min_n) or (max_n < nf) or any(abs(ddx) > 10 * _EPS * (L))): # pab 07.01.2001: Always choose the stepsize df so that # it is an exactly representable number. # This is important when calculating numerical derivatives and is @@ -2439,8 +2476,6 @@ def tranproc(x, f, x0, *xi): -------- trangood. """ - - eps = floatinfo.eps xo, fo, x0 = atleast_1d(x, f, x0) xi = atleast_1d(*xi) if not isinstance(xi, list): @@ -2464,10 +2499,11 @@ def tranproc(x, f, x0, *xi): if N > 0: y = [y0] hn = xo[1] - xo[0] - if hn ** N < sqrt(eps): - warnings.warn('Numerical problems may occur for the derivatives' + - ' in tranproc. The sampling of the transformation' + - ' may be too small.') + if hn ** N < sqrt(_EPS): + msg = ('Numerical problems may occur for the derivatives in ' + + 'tranproc.\n' + + 'The sampling of the transformation may be too small.') + warnings.warn(msg) # Transform X with the derivatives of f. fxder = zeros((N, x0.size)) @@ -2504,8 +2540,8 @@ def tranproc(x, f, x0, *xi): fxder[1] * (3. * xi[1] ** 2. + 4. * xi[0] * xi[1])) y.append(y4) if N > 4: - warnings.warn('Transformation of derivatives of' + - ' order>4 not supported.') + warnings.warn('Transformation of derivatives of ' + + 'order>4 not supported.') return y # y0,y1,y2,y3,y4 @@ -2517,10 +2553,11 @@ def good_bins(data=None, range=None, num_bins=None, # @ReservedAssignment ---------- data : array-like the data - range : (float, float), (default data.min(), data.max()) - minimum and maximum range of bins - num_bins : scalar integer, (default depending on num_data=len(data)) + range : (float, float) + minimum and maximum range of bins (default data.min(), data.max()) + num_bins : scalar integer approximate number of bins wanted + (default depending on num_data=len(data)) odd : bool placement of bins (0 or 1) (default 0) loose : bool @@ -2621,14 +2658,14 @@ def plot_histgrm(data, bins=None, range=None, # @ReservedAssignment if bins is None: bins = np.ceil(4 * np.sqrt(np.sqrt(len(x)))) - # , new=True) @ReservedAssignment - y, limits = np.histogram(data, bins=bins, normed=normed, weights=weights) + bin_, limits = np.histogram( + data, bins=bins, normed=normed, weights=weights) limits.shape = (-1, 1) xx = limits.repeat(3, axis=1) xx.shape = (-1,) xx = xx[1:-1] - y.shape = (-1, 1) - yy = y.repeat(3, axis=1) + bin_.shape = (-1, 1) + yy = bin_.repeat(3, axis=1) # yy[0,0] = 0.0 # pdf yy[:, 0] = 0.0 # histogram yy.shape = (-1,) @@ -2719,8 +2756,10 @@ def fourier(data, t=None, T=None, m=None, n=None, method='trapz'): >>> t = np.linspace(0,4*T) >>> x = np.sin(t) >>> a, b = wm.fourier(x, t, T=T, m=5) - >>> (np.round(a.ravel()), np.round(b.ravel())) - (array([ 0., 0., 0., 0., 0.]), array([ 0., 4., 0., 0., 0.])) + >>> np.abs(a.ravel())<1e-12 + array([ True, True, True, True, True], dtype=bool) + >>> np.abs(b.ravel()-np.array([ 0., 4., 0., 0., 0.]))<1e-12 + array([ True, True, True, True, True], dtype=bool) See also -------- @@ -2778,93 +2817,9 @@ def fourier(data, t=None, T=None, m=None, n=None, method='trapz'): return a, b -def _test_find_cross(): - t = findcross([0, 0, 1, -1, 1], 0) # @UnusedVariable - - -def _test_common_shape(): - - A = ones((4, 1)) - B = 2 - C = ones((1, 5)) * 5 - common_shape(A, B, C) - - common_shape(A, B, C, shape=(3, 4, 1)) - - A = ones((4, 1)) - B = 2 - C = ones((1, 5)) * 5 - common_shape(A, B, C, shape=(4, 5)) - - -def _test_meshgrid(): - x = array([-1, -0.5, 1, 4, 5], float) - y = array([0, -2, -5], float) - xv, yv = meshgrid(x, y, sparse=False) - print(xv) - print(yv) - xv, yv = meshgrid(x, y, sparse=True) # make sparse output arrays - print(xv) - print(yv) - print(meshgrid(0, 1, 5, sparse=True)) # just a 3D point - print(meshgrid([0, 1, 5], sparse=True)) # just a 3D point - xv, yv = meshgrid(y, y) - yv[0, 0] = 10 - print(xv) - print(yv) -# >>> xv -# array([[ 0. , 0.5, 1. ]]) -# >>> yv -# array([[ 0.], -# [ 1.]]) -# array([[-1. , -0.5, 1. , 4. , 5. ], -# [-1. , -0.5, 1. , 4. , 5. ], -# [-1. , -0.5, 1. , 4. , 5. ]]) -## -# array([[ 0., 0., 0., 0., 0.], -# [-2., -2., -2., -2., -2.], -# [-5., -5., -5., -5., -5.]]) - - -def _test_tranproc(): - # import utilitiestransform.models as wtm - tr = lambda x: x # wtm.TrHermite() - x = linspace(-5, 5, 501) - g = tr(x) - _gder = tranproc(x, g, x, ones(g.size)) - pass - # >>> gder(:,1) = g(:,1) - # >>> plot(g(:,1),[g(:,2),gder(:,2)]) - # >>> plot(g(:,1),pdfnorm(g(:,2)).*gder(:,2),g(:,1),pdfnorm(g(:,1))) - # >>> legend('Transformed model','Gaussian model') - -def _test_detrend(): - import pylab as plt - cos = np.cos - randn = np.random.randn - x = linspace(0, 1, 200) - y = exp(x) + cos(5 * 2 * pi * x) + 1e-1 * randn(x.size) - y0 = detrendma(y, 20) - tr = y - y0 - plt.plot(x, y, x, y0, 'r', x, exp(x), 'k', x, tr, 'm') -def _test_extrema(): - import pylab as plt - from pylab import plot - t = plt.linspace(0, 7 * pi, 250) - x = plt.sin(t) + 0.1 * sin(50 * t) - ind = findextrema(x) - ti, tp = t[ind], x[ind] - plot(t, x, '.', ti, tp, 'r.') - _ind1 = findrfc(tp, 0.3) - - -def _test_discretize(): - import pylab as plt - x, y = discretize(cos, 0, pi) - plt.plot(x, y) plt.show() plt.close('all') @@ -2915,14 +2870,14 @@ def profile_main1(): import pstats prof = cProfile.Profile() prof = prof.runctx("real_main()", globals(), locals()) - print "
"
+    print("
")
     stats = pstats.Stats(prof)
     stats.sort_stats("time")  # Or cumulative
     stats.print_stats(80)  # 80 = how many to print
     # The rest is optional.
     # stats.print_callees()
     # stats.print_callers()
-    print "
"
+    print("
") main = profile_main1 diff --git a/wafo/padua.py b/wafo/padua.py index d43f7a4..afa69d9 100644 --- a/wafo/padua.py +++ b/wafo/padua.py @@ -274,7 +274,7 @@ class _ExampleFunctions(object): arg_z = 1. / arg_z return np.exp(-arg_z) - def __call__(self, x, y, id=0): + def __call__(self, x, y, id=0): # @ReservedAssignment s = self test_function = [s.franke, s.half_sphere, s.poly_degree20, s.exp_fun1, s.exp_fun100, s.cos30, s.constant, s.exp_xy, s.runge, @@ -410,7 +410,7 @@ def paduavals2coefs(f): else: # dct = @(c) chebtech2.coeffs2vals(c); - C = np.rot90(dct(dct(G.T).T)) #, axis=1) + C = np.rot90(dct(dct(G.T).T)) # , axis=1) C[0] = .5 * C[0] C[:, 1] = .5 * C[:, 1] diff --git a/wafo/polynomial.py b/wafo/polynomial.py index 1a1efa3..bb35f82 100644 --- a/wafo/polynomial.py +++ b/wafo/polynomial.py @@ -17,14 +17,14 @@ # Licence: LGPL # ------------------------------------------------------------------------- # !/usr/bin/env python -import warnings +import warnings # @UnusedImport from numpy.polynomial import polyutils as pu from plotbackend import plotbackend as plt import numpy as np -from numpy import (zeros, ones, zeros_like, array, asarray, newaxis, arange, - logical_or, any, pi, cos, round, diff, all, exp, atleast_1d, - where, extract, linalg, sign, concatenate, floor, isreal, - conj, remainder, linspace, sum, meshgrid, hstack) +from numpy import (zeros, asarray, newaxis, arange, + 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 numpy.lib.polynomial import * # @UnusedWildImport @@ -110,9 +110,9 @@ def polyint(p, m=1, k=None): raise ValueError("Order of integral must be positive (see polyder)") if k is None: k = zeros(m, float) - k = atleast_1d(k) + k = np.atleast_1d(k) if len(k) == 1 and m > 1: - k = k[0] * ones(m, float) + k = k[0] * np.ones(m, float) if len(k) < m: raise ValueError( "k must be a scalar or a rank-1 array of length 1 or >m.") @@ -127,7 +127,7 @@ def polyint(p, m=1, k=None): if p.ndim > 1: ix = ix[..., newaxis] pieces = p.shape[-1] - k0 = k[0] * ones((1, pieces), dtype=int) + k0 = k[0] * np.ones((1, pieces), dtype=int) else: k0 = [k[0]] y = np.concatenate((p.__truediv__(ix), k0), axis=0) @@ -260,8 +260,8 @@ def polydeg(x, y): # developed in a series of orthogonal polynomials. ys = np.ones((N,)) * y.mean() # correction for small sample sizes - AIC = 2 + N * \ - (np.log(2 * pi * ((ys - y) ** 2).sum() / N) + 1) + 4 / (N - 2) + logsum2 = (np.log(2 * pi * ((ys - y) ** 2).sum() / N) + 1) + AIC = 2 + N * logsum2 + 4 / (N - 2) n = 1 nit = 0 @@ -495,7 +495,7 @@ def polyreloc(p, x, y=0.0): """ truepoly = isinstance(p, poly1d) - r = atleast_1d(p).copy() + r = np.atleast_1d(p).copy() n = r.shape[0] # Relocate polynomial using Horner's algorithm @@ -548,7 +548,7 @@ def polyrescl(p, x, y=1.0): """ truepoly = isinstance(p, poly1d) - r = atleast_1d(p) + r = np.atleast_1d(p) n = r.shape[0] xscale = (float(x) ** arange(1 - n, 1)) @@ -591,7 +591,7 @@ def polytrim(p): if truepoly: return p else: - r = atleast_1d(p).copy() + r = np.atleast_1d(p).copy() # Remove leading zeros is_not_lead_zeros = logical_or.accumulate(r != 0, axis=0) if r.ndim == 1: @@ -630,7 +630,7 @@ def poly2hstr(p, variable='x'): """ var = variable - coefs = polytrim(atleast_1d(p)) + coefs = polytrim(np.atleast_1d(p)) order = len(coefs) - 1 # Order of polynomial. s = '' # Initialize output string. ix = 1 @@ -719,7 +719,7 @@ def poly2str(p, variable='x'): var = variable # Remove leading zeros - coeffs = polytrim(atleast_1d(p)) + coeffs = polytrim(np.atleast_1d(p)) N = len(coeffs) - 1 @@ -1094,7 +1094,7 @@ def chebpoly(n, x=None, kind=1): """ if x is None: # Calculate coefficients. if n == 0: - p = ones(1) + p = np.ones(1) else: p = round(pow(2, n - 2 + kind) * poly(chebroot(n, kind=kind))) p[1::2] = 0 @@ -1102,7 +1102,7 @@ def chebpoly(n, x=None, kind=1): else: # Evaluate polynomial in chebychev form ck = zeros(n + 1) ck[0] = 1. - return _chebval(atleast_1d(x), ck, kind=kind) + return _chebval(np.atleast_1d(x), ck, kind=kind) def chebfit(fun, n=10, a=-1, b=1, trace=False): @@ -1188,7 +1188,8 @@ def chebfit_dct(f, n=(10, ), domain=None): Fit Chebyshev series to N-dimensional function so that f(x1, x2,..., xn) can be approximated by: - .. math:: f(x_1, x_2,...,x_n) = \\sum_{i,j,...k} c_i T_i(x_1)*...*c_k T_k(x_n) , + .. math:: f(x_1, x_2,...,x_n) = + \\sum_{i,j,...k} c_i T_i(x_1)*...*c_k T_k(x_n) , where Tk is the k'th Chebyshev polynomial of the first kind. @@ -1251,10 +1252,10 @@ def chebfit_dct(f, n=(10, ), domain=None): if hasattr(f, '__call__'): if domain is None: - domain = (-1,1) * len(n) + domain = (-1, 1) * len(n) domain = np.atleast_2d(domain).reshape((-1, 2)) xi = [map_to_interval(chebroot(ni), d[0], d[1]) - for ni, d in zip(n, domain)] + for ni, d in zip(n, domain)] Xi = np.meshgrid(*xi) ck = f(*Xi) else: @@ -1263,7 +1264,7 @@ def chebfit_dct(f, n=(10, ), domain=None): ndim = len(n) for i in range(ndim): - ck = dct(ck[...,::-1]) + ck = dct(ck[..., ::-1]) ck[..., 0] = ck[..., 0] / 2. if i < ndim-1: ck = np.rollaxis(ck, axis=-1) @@ -1591,8 +1592,8 @@ class Cheb1d(object): def __eq__(self, other): other = Cheb1d(other) - return (all(self.coeffs == other.coeffs) and (self.a == other.a) - and (self.b == other.b) and (self.kind == other.kind)) + return (all(self.coeffs == other.coeffs) and (self.a == other.a) and + (self.b == other.b) and (self.kind == other.kind)) def __ne__(self, other): return any(self.coeffs != other.coeffs) or (self.a != other.a) or ( @@ -1830,8 +1831,8 @@ def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True): if trace: plt.plot(x, fs, '+') - wt = ones((npt)) - ee = ones((npt)) + wt = np.ones((npt)) + ee = np.ones((npt)) mad = 0 u = zeros((npt, ncof)) @@ -1989,37 +1990,40 @@ def chebvandernd(deg, *xi): ------- vander : ndarray The shape of the returned matrix is ``x1.shape + (order,)``, where - :math:`order = (deg[0]+1)*(deg([1]+1)*...*(deg[n-1]+1)`. The dtype will - be the same as the converted `x1`, `x2`, ... `xn`. + :math:`order = (deg[0]+1)*(deg([1]+1)*...*(deg[n-1]+1)`. The dtype + will be the same as the converted `x1`, `x2`, ... `xn`. See Also -------- chebvander, chebvalnd, chebfitnd """ ideg = [int(d) for d in deg] - is_valid = np.array([id == d and id >= 0 for id, d in zip(ideg, deg)]) + is_valid = np.array([di == d and di >= 0 for di, d in zip(ideg, deg)]) if np.any(is_valid != 1): raise ValueError("degrees must be non-negative integers") ndim = len(xi) - if len(ideg)!=ndim: - raise ValueError('length of deg must be the same as number of dimensions') + if len(ideg) != ndim: + msg = 'length of deg must be the same as number of dimensions' + raise ValueError(msg) xi = np.array(xi, copy=0) + 0.0 chebvander = np.polynomial.chebyshev.chebvander shape0 = xi[0].shape s0 = (1,) * ndim vxi = [chebvander(x, d).reshape(shape0 + s0[:i] + (-1,) + s0[i + 1::]) - for i, (d, x) in enumerate(zip(ideg, xi))] + for i, (d, x) in enumerate(zip(ideg, xi))] v = reduce(np.multiply, vxi) return v.reshape(v.shape[:-ndim] + (-1,)) + def chebfitnd(xi, f, deg, rcond=None, full=False, w=None): """ Least squares fit of Chebyshev series to N-dimensional data. Return the coefficients of a Chebyshev series of degree `deg` that is the - least squares fit to the data values `f` given at points `x1`, `x2`,..., `xn` + least squares fit to the data values `f` given at points + `x1`, `x2`,..., `xn` The fitted polynomial(s) are in the form .. math:: p(x,y) = c_00 + c_11 * T_1(x)*T_1(y) + ..c_ij * T_i(x)*T_j(y). @@ -2033,7 +2037,8 @@ def chebfitnd(xi, f, deg, rcond=None, full=False, w=None): f : array_like function values at the sample points ``(x1[i], x2[i], ..., xn[i])``. deg : list - Degrees of the fitting series in the x1, x2, ..., xn directions, respectively. + Degrees of the fitting series in the x1, x2, ..., xn directions, + respectively. rcond : float, optional Relative condition number of the fit. Singular values smaller than this relative to the largest singular value will be ignored. The @@ -2103,7 +2108,7 @@ def chebfitnd(xi, f, deg, rcond=None, full=False, w=None): Examples -------- """ - xi_ = np.array(xi, copy=0) + 0.0 + # xi = np.array(xi, copy=0) + 0.0 z = np.array(f) degrees = np.asarray(deg, dtype=int) orders = degrees + 1 @@ -2112,7 +2117,7 @@ def chebfitnd(xi, f, deg, rcond=None, full=False, w=None): ndims = np.array([x.ndim for x in xi]) ndim = len(ndims) sizes = np.array([x.size for x in xi]) - if np.any(ndims!=ndim) or z.ndim!=ndim: + if np.any(ndims != ndim) or z.ndim != ndim: raise TypeError("expected %dD array for x1, x2,...,xn and f" % ndim) if np.any(sizes == 0): raise TypeError("expected non-empty vector for xi") @@ -2148,13 +2153,15 @@ def chebfitnd(xi, f, deg, rcond=None, full=False, w=None): else: return c + def chebvalnd(c, *xi): """ Evaluate a N-D Chebyshev series at points (x1, x2, ..., xn). This function returns the values: - .. math:: p(x1,x2,...,xn) = \\sum_{i,j,...,k} c_{i,j,...,k} * T_i(x1) * T_j(x2)*...* T_k(xn) + .. math:: p(x1,x2,...,xn) = + \\sum_{i,j,...,k} c_{i,j,...,k} * T_i(x1) * T_j(x2)*...* T_k(xn) The parameters `x1`, `x2`, ...., `xn` are converted to arrays only if they are tuples or a lists, otherwise they are treated as a scalars and @@ -2171,14 +2178,14 @@ def chebvalnd(c, *xi): c : array_like Array of coefficients ordered so that the coefficient of the term of multi-degree i,j,...,k is contained in ``c[i,j,...,k]``. If `c` has - dimension greater than N the remaining indices enumerate multiple sets of - coefficients. + dimension greater than N the remaining indices enumerate multiple sets + of coefficients. x1, x2,..., xn : array_like, compatible object The N dimensional series is evaluated at the points `(x1, x2,...,xn)`, where `x1`, `x2`,..., `xn` must have the same shape. - If any of `x1`, `x2`, ..., `xn` is a list or tuple, it is first converted - to an ndarray, otherwise it is left unchanged and if it isn't an - ndarray it is treated as a scalar. + If any of `x1`, `x2`, ..., `xn` is a list or tuple, it is first + converted to an ndarray, otherwise it is left unchanged and if it isn't + an ndarray it is treated as a scalar. Returns ------- @@ -2194,12 +2201,13 @@ def chebvalnd(c, *xi): xi = np.array(xi, copy=0) except: raise ValueError('x, y, z are incompatible') - chebval = np.polynomial.chebyshev.chebval + chebval = np.polynomial.chebyshev.chebval c = chebval(xi[0], c) for x in xi[1:]: c = chebval(x, c, tensor=False) return c + def chebgridnd(c, *xi): """ Evaluate a N-D Chebyshev series on the Cartesian product of x1, x2,..., xn. @@ -2212,8 +2220,8 @@ def chebgridnd(c, *xi): `a` from `x1`, `b` from `x2`, and so on. The resulting points form a grid with `x1` in the first dimension, `x2` in the second, and so on. - The parameters `x1`, `x2`, ... and `xn` are converted to arrays only if they - are tuples or a lists, otherwise they are treated as a scalars. In + The parameters `x1`, `x2`, ... and `xn` are converted to arrays only if + they are tuples or a lists, otherwise they are treated as a scalars. In either case, either `x1`, `x2`,... and `xn` or their elements must support multiplication and addition both with themselves and with the elements of `c`. @@ -2247,26 +2255,25 @@ def chebgridnd(c, *xi): -------- chebval, chebvalnd, chebfitnd """ - chebval = np.polynomial.chebyshev.chebval + chebval = np.polynomial.chebyshev.chebval for x in xi: c = chebval(x, c) return c -def test_chebfit1d(): - n = 63 - x = chebroot(n=64, kind=1) +def test_chebfit1d(): def f(x): return np.exp(-x**2) - z = f(x) + # x = chebroot(n=64, kind=1) + # z = f(x) c = chebfit(f, n=64)[::-1] xi = np.linspace(-1, 1, 151) zi = np.polynomial.chebyshev.chebval(xi, c) - #plt.plot(xi, zi,'.', xi, f(xi)) + # plt.plot(xi, zi,'.', xi, f(xi)) plt.semilogy(xi, np.abs(zi-f(xi))) plt.show('hold') @@ -2276,28 +2283,30 @@ def test_chebfit2d(): xorder, yorder = n-1, n-1 x = chebroot(n=n, kind=1) xgrid, ygrid = meshgrid(x, x) + def f(x, y): return np.exp(-x**2-6*y**2) zgrid = f(xgrid, ygrid) - #v2d = np.polynomial.chebyshev.chebvander2d(xgrid, ygrid, [xorder,yorder]).reshape((-1, (xorder+1)*(yorder+1))) - #coeff, residuals, rank, s = np.linalg.lstsq(v2d, zgrid.ravel()) - #doeff = coeff.reshape(xorder+1,yorder+1) - dcoeff2 = chebfitnd((xgrid, ygrid), zgrid, [xorder,yorder]) - dcoeff = chebfit_dct(f, n=(xorder+1,yorder+1)) + # v2d = np.polynomial.chebyshev.chebvander2d(xgrid, ygrid, + # [xorder,yorder]).reshape((-1, (xorder+1)*(yorder+1))) + # coeff, residuals, rank, s = np.linalg.lstsq(v2d, zgrid.ravel()) + # doeff = coeff.reshape(xorder+1,yorder+1) + _dcoeff2 = chebfitnd((xgrid, ygrid), zgrid, [xorder, yorder]) + dcoeff = chebfit_dct(f, n=(xorder+1, yorder+1)) xi = np.linspace(-1, 1, 151) - Xi,Yi = np.meshgrid(xi, xi) + Xi, Yi = np.meshgrid(xi, xi) Zi = f(Xi, Yi) zzi = chebvalnd(dcoeff, Xi, Yi) - devi = Zi - zzi + _devi = Zi - zzi # plot residuals - #zz = np.polynomial.chebyshev.chebval2d(xgrid, ygrid, dcoeff) + # zz = np.polynomial.chebyshev.chebval2d(xgrid, ygrid, dcoeff) zz = chebvalnd(dcoeff, xgrid, ygrid) dev = zgrid - zz - #plt.spy(np.abs(dcoeff)>1e-13) + # plt.spy(np.abs(dcoeff)>1e-13) plt.contourf(xgrid, ygrid, np.abs(dev)) - #plt.contourf(Xi, Yi, np.abs(devi)) + # plt.contourf(Xi, Yi, np.abs(devi)) plt.colorbar() # plt.semilogy(np.abs(devi.ravel())) plt.show('hold') diff --git a/wafo/polynomial_old.py b/wafo/polynomial_old.py deleted file mode 100644 index 4f78206..0000000 --- a/wafo/polynomial_old.py +++ /dev/null @@ -1,1207 +0,0 @@ -#------------------------------------------------------------------------------- -# Name: polynomial -# Purpose: Functions to operate on polynomials. -# -# Author: pab -# polyXXX functions are based on functions found in the matlab toolbox polyutil written by -# Author: Peter J. Acklam -# E-mail: pjacklam@online.no -# WWW URL: http://home.online.no/~pjacklam -# -# Created: 30.12.2008 -# Copyright: (c) pab 2008 -# Licence: LGPL -#------------------------------------------------------------------------------- -#!/usr/bin/env python - -""" - Extended functions to operate on polynomials -""" -import warnings -import numpy as np -from numpy.lib.polynomial import * -__all__ = np.lib.polynomial.__all__ -__all__ = __all__ + ['polyreloc', 'polyrescl', 'polytrim', 'poly2hstr', 'poly2str', - 'polyshift', 'polyishift', 'map_from_intervall', 'map_to_intervall', - 'cheb2poly', 'chebextr', 'chebroot', 'chebpoly', 'chebfit', 'chebval', - 'chebder', 'chebint', 'Cheb1d', 'dct', 'idct'] - -def polyreloc(p,x,y=0.0): - """ - Relocate polynomial - - The polynomial `p` is relocated by "moving" it `x` - units along the x-axis and `y` units along the y-axis. - So the polynomial `r` is relative to the point (x,y) as - the polynomial `p` is relative to the point (0,0). - - Parameters - ---------- - p : array-like, poly1d - vector or matrix of column vectors of polynomial coefficients to relocate. - (Polynomial coefficients are in decreasing order.) - x : scalar - distance to relocate P along x-axis - y : scalar - distance to relocate P along y-axis (default 0) - - Returns - ------- - r : ndarray, poly1d - vector/matrix/poly1d of relocated polynomial coefficients. - - See also - -------- - polyrescl - - Example - ------- - >>> import numpy as np - >>> p = np.arange(6); p.shape = (2,-1) - >>> np.polyval(p,0) - array([3, 4, 5]) - >>> np.polyval(p,1) - array([3, 5, 7]) - >>> r = polyreloc(p,-1) # move to the left along x-axis - >>> np.polyval(r,-1) # = polyval(p,0) - array([3, 4, 5]) - >>> np.polyval(r,0) # = polyval(p,1) - array([3, 5, 7]) - """ - - truepoly = isinstance(p, poly1d) - r = np.atleast_1d(p).copy() - n = r.shape[0] - - # Relocate polynomial using Horner's algorithm - for ii in range(n,1,-1): - for i in range(1,ii): - r[i] = r[i] - x*r[i-1] - r[-1] = r[-1] + y - if r.ndim>1 and r.shape[-1]==1: - r.shape = (r.size,) - if truepoly: - r = poly1d(r) - return r - -def polyrescl(p, x, y=1.0): - """ - Rescale polynomial. - - Parameters - ---------- - p : array-like, poly1d - vector or matrix of column vectors of polynomial coefficients to rescale. - (Polynomial coefficients are in decreasing order.) - x,y : scalars - defining the factors to rescale the polynomial `p` in - x-direction and y-direction, respectively. - - Returns - ------- - r : ndarray, poly1d - vector/matrix/poly1d of rescaled polynomial coefficients. - - See also - -------- - polyreloc - - Example - ------- - >>> import numpy as np - >>> p = np.arange(6); p.shape = (2,-1) - >>> np.polyval(p,0) - array([3, 4, 5]) - >>> np.polyval(p,1) - array([3, 5, 7]) - >>> r = polyrescl(p,2) # scale by 2 along x-axis - >>> np.polyval(r,0) # = polyval(p,0) - array([ 3., 4., 5.]) - >>> np.polyval(r,2) # = polyval(p,1) - array([ 3., 5., 7.]) - """ - - truepoly = isinstance(p, poly1d) - r = np.atleast_1d(p).copy() - n = r.shape[0] - - xscale =(float(x)**np.arange(1-n , 1)) - if r.ndim==1: - q = y*r*xscale - else: - q = y*r*xscale[:,np.newaxis] - if truepoly: - q = poly1d(q) - return q - -def polytrim(p): - """ - Trim polynomial by stripping off leading zeros. - - Parameters - ---------- - p : array-like, poly1d - vector of polynomial coefficients in decreasing order. - - Returns - ------- - r : ndarray, poly1d - vector/matrix/poly1d of trimmed polynomial coefficients. - - Example - ------- - >>> p = [0,1,2] - >>> polytrim(p) - array([1, 2]) - """ - - truepoly = isinstance(p, poly1d) - if truepoly: - return p - else: - r = np.atleast_1d(p).copy() - # Remove leading zeros - is_not_lead_zeros =np.logical_or.accumulate(r != 0,axis=0) - if r.ndim==1: - r = r[is_not_lead_zeros] - else: - is_not_lead_zeros = np.any(is_not_lead_zeros,axis=1) - r = r[is_not_lead_zeros,:] - return r - -def poly2hstr(p, variable='x' ): - """ - Return polynomial as a Horner represented string. - - Parameters - ---------- - p : array-like poly1d - vector of polynomial coefficients in decreasing order. - variable : string - display character for variable - - Returns - ------- - p_str : string - consisting of the polynomial coefficients in the vector P multiplied - by powers of the given `variable`. - - Examples - -------- - >>> poly2hstr([1, 1, 2], 's' ) - '(s + 1)*s + 2' - - See also - -------- - poly2str - """ - var = variable - - coefs = polytrim(np.atleast_1d(p)) - order = len(coefs)-1 # Order of polynomial. - s = '' # Initialize output string. - ix = 1; - for expon in range(order,-1,-1): - coef = coefs[order-expon] - #% There is no point in adding a zero term (except if it's the only - #% term, but we'll take care of that later). - if coef == 0: - ix = ix+1 - else: - #% Append exponent if necessary. - if ix>1: - exponstr = '%.0f' % ix - s = '%s**%s' % (s,exponstr); - ix = 1 - #% Is it the first term? - isfirst = s == '' - - # We need the coefficient only if it is different from 1 or -1 or - # when it is the constant term. - needcoef = (( abs(coef) != 1 ) | ( expon == 0 ) & isfirst) | 1-isfirst - - # We need the variable except in the constant term. - needvar = ( expon != 0 ) - - #% Add sign, but we don't need a leading plus-sign. - if isfirst: - if coef < 0: - s = '-' # % Unary minus. - else: - if coef < 0: - s = '%s - ' % s # % Binary minus (subtraction). - else: - s = '%s + ' % s # % Binary plus (addition). - - - #% Append the coefficient if it is different from one or when it is - #% the constant term. - if needcoef: - coefstr = '%.20g' % abs(coef) - s = '%s%s' % (s,coefstr) - - #% Append variable if necessary. - if needvar: - #% Append a multiplication sign if necessary. - if needcoef: - if 1-isfirst: - s = '(%s)' % s - s = '%s*' % s - s = '%s%s' % ( s, var ) - - #% Now treat the special case where the polynomial is zero. - if s=='': - s = '0' - return s - -def poly2str(p,variable='x'): - """ - Return polynomial as a string. - - Parameters - ---------- - p : array-like poly1d - vector of polynomial coefficients in decreasing order. - variable : string - display character for variable - - Returns - ------- - p_str : string - consisting of the polynomial coefficients in the vector P multiplied - by powers of the given `variable`. - - See also - -------- - poly2hstr - - Examples - -------- - >>> poly2str([1, 1, 2], 's' ) - 's**2 + s + 2' - """ - thestr = "0" - var = variable - - # Remove leading zeros - coeffs = polytrim(np.atleast_1d(p)) - - N = len(coeffs)-1 - - for k in range(len(coeffs)): - coefstr = '%.4g' % abs(coeffs[k]) - if coefstr[-4:] == '0000': - coefstr = coefstr[:-5] - power = (N-k) - if power == 0: - if coefstr != '0': - newstr = '%s' % (coefstr,) - else: - if k == 0: - newstr = '0' - else: - newstr = '' - elif power == 1: - if coefstr == '0': - newstr = '' - elif coefstr == 'b' or coefstr == '1': - newstr = var - else: - newstr = '%s*%s' % (coefstr, var) - else: - if coefstr == '0': - newstr = '' - elif coefstr == 'b' or coefstr == '1': - newstr = '%s**%d' % (var, power,) - else: - newstr = '%s*%s**%d' % (coefstr, var, power) - - if k > 0: - if newstr != '': - if coeffs[k] < 0: - thestr = "%s - %s" % (thestr, newstr) - else: - thestr = "%s + %s" % (thestr, newstr) - elif (k == 0) and (newstr != '') and (coeffs[k] < 0): - thestr = "-%s" % (newstr,) - else: - thestr = newstr - return thestr - -def polyshift(py,a=-1,b=1): - """ - Polynomial coefficient shift - - Polyshift shift the polynomial coefficients by a variable shift: - - Y = 2*(X-.5*(b+a))/(b-a) - - i.e., the interval -1 <= Y <= 1 is mapped to the interval a <= X <= b - - Parameters - ---------- - py : array-like - polynomial coefficients for the variable y. - a,b : scalars - lower and upper limit. - - Returns - ------- - px : ndarray - polynomial coefficients for the variable x. - - See also - -------- - polyishift - - Example - ------- - >>> py = [1, 0] - >>> px = polyshift(py,0,5) - >>> polyval(px,[0, 2.5, 5]) #% This is the same as the line below - array([-1., 0., 1.]) - >>> polyval(py,[-1, 0, 1 ]) - array([-1, 0, 1]) - """ - - if (a==-1) & (b ==1): - return py - L = b-a - return polyishift(py,-(2.+b+a)/L,(2.-b-a)/L) - -def polyishift(px,a=-1,b=1): - """ - Inverse polynomial coefficient shift - - Polyishift does the inverse of Polyshift, - shift the polynomial coefficients by a variable shift: - - Y = 2*(X-.5*(b+a)/(b-a) - - i.e., the interval a <= X <= b is mapped to the interval -1 <= Y <= 1 - - Parameters - ---------- - px : array-like - polynomial coefficients for the variable x. - a,b : scalars - lower and upper limit. - - Returns - ------- - py : ndarray - polynomial coefficients for the variable y. - - See also - -------- - polyishift - - Example - ------- - >>> px = [1, 0] - >>> py = polyishift(px,0,5); - >>> polyval(px,[0, 2.5, 5]) #% This is the same as the line below - array([ 0. , 2.5, 5. ]) - >>> polyval(py,[-1, 0, 1]) - array([ 0. , 2.5, 5. ]) - """ - if (a==-1) & (b ==1): - return px - L = b-a - xscale = 2./L - xloc = -float(a+b)/L - return polyreloc( polyrescl(px,xscale),xloc) - -def map_from_interval(x,a,b) : - """F(x), where F: [a,b] -> [-1,1].""" - return (x - (b + a)/2.0)*(2.0/(b - a)) - -def map_to_interval(x,a,b) : - """F(x), where F: [-1,1] -> [a,b].""" - return (x*(b - a) + (b + a))/2.0 - -def poly2cheb(p,a=-1,b=1): - """ - Convert polynomial coefficients into Chebyshev coefficients - - Parameters - ---------- - p : array-like - polynomial coefficients - a,b : real scalars - lower and upper limits (Default -1,1) - - Returns - ------- - ck : ndarray - Chebychef coefficients - - POLY2CHEB do the inverse of CHEB2POLY: given a vector of polynomial - coefficients AK, returns an equivalent vector of Chebyshev - coefficients CK. - - This is useful for economization of power series. - The steps for doing so: - 1. Convert polynomial coefficients to Chebychev coefficients, CK. - 2. Truncate the CK series to a smaller number of terms, using the - coefficient of the first neglected Chebychev polynomial as an error - estimate. - 3 Convert back to a polynomial by CHEB2POLY - - See also - -------- - cheb2poly - chebval - chebfit - - Examples - -------- - >>> import numpy as np - >>> p = np.arange(5) - >>> ck = poly2cheb(p) - >>> cheb2poly(ck) - array([ 1., 2., 3., 4.]) - - Reference - --------- - William H. Press, Saul Teukolsky, - William T. Wetterling and Brian P. Flannery (1997) - "Numerical recipes in Fortran 77", Vol. 1, pp 184-194 - """ - f = poly1d(p) - n = len(f.coeffs) - return chebfit(f,n,a,b) - -def cheb2poly(ck,a=-1,b=1): - """ - Converts Chebyshev coefficients to polynomial coefficients - - Parameters - ---------- - ck : array-like - Chebychef coefficients - a,b : real, scalars - lower and upper limits (Default -1,1) - - Returns - ------- - p : ndarray - polynomial coefficients - - It is not advised to do this for len(ck)>10 due to numerical cancellations. - - See also - -------- - chebval - chebfit - - Examples - -------- - >>> import numpy as np - >>> p = np.arange(5) - >>> ck = poly2cheb(p) - >>> cheb2poly(ck) - array([ 1., 2., 3., 4.]) - - - References - ---------- - http://en.wikipedia.org/wiki/Chebyshev_polynomials - http://en.wikipedia.org/wiki/Chebyshev_form - http://en.wikipedia.org/wiki/Clenshaw_algorithm - """ - - n = len(ck) - - b_Nmi = np.zeros(1) - b_Nmip1 = np.zeros(1) - y = np.r_[2/(b-a), -(a+b)/(b-a)] - y2 = 2.*y - - # Clenshaw recurence - for ix in range(n-1,0,-1): - tmp = b_Nmi - b_Nmi = polymul(y2,b_Nmi) # polynomial multiplication - nb = len(b_Nmip1) - b_Nmip1[-1] = b_Nmip1[-1]-ck[ix] - b_Nmi[-nb::] = b_Nmi[-nb::]-b_Nmip1 - b_Nmip1 = tmp - - p = polymul(y,b_Nmi) # polynomial multiplication - nb = len(b_Nmip1) - b_Nmip1[-1] = b_Nmip1[-1]-ck[0] - p[-nb::] = p[-nb::]-b_Nmip1 - return polytrim(p) - -def chebextr(n): - """ - Return roots of derivative of Chebychev polynomial of the first kind. - - All local extreme values of the polynomial are either -1 or 1. So, - CHEBPOLY( N, CHEBEXTR(N) ) ) return the same as (-1).^(N:-1:0) - except for the numerical noise in the former. - - Because the extreme values of Chebychev polynomials of the first - kind are either -1 or 1, their roots are often used as starting - values for the nodes in minimax approximations. - - - Parameters - ---------- - n : scalar, integer - degree of Chebychev polynomial. - - Examples - -------- - >>> x = chebextr(4) - >>> chebpoly(4,x) - array([ 1., -1., 1., -1., 1.]) - - - Reference - --------- - http://en.wikipedia.org/wiki/Chebyshev_nodes - http://en.wikipedia.org/wiki/Chebyshev_polynomials - """ - return -np.cos((np.pi*np.arange(n+1))/n); - -def chebroot(n,kind=1): - """ - Return roots of Chebychev polynomial of the first or second kind. - - The roots of the Chebychev polynomial of the first kind form a particularly - good set of nodes for polynomial interpolation because the resulting - interpolation polynomial minimizes the problem of Runge's phenomenon. - - Parameters - ---------- - n : scalar, integer - degree of Chebychev polynomial. - kind: 1 or 2, optional - kind of Chebychev polynomial (default 1) - - Examples - -------- - >>> import numpy as np - >>> x = chebroot(3) - >>> np.abs(chebpoly(3,x))<1e-15 - array([ True, True, True], dtype=bool) - >>> chebpoly(3) - array([ 4., 0., -3., 0.]) - >>> x2 = chebroot(4,kind=2) - >>> np.abs(chebpoly(4,x2,kind=2))<1e-15 - array([ True, True, True, True], dtype=bool) - >>> chebpoly(4,kind=2) - array([ 16., 0., -12., 0., 1.]) - - - Reference - --------- - http://en.wikipedia.org/wiki/Chebyshev_nodes - http://en.wikipedia.org/wiki/Chebyshev_polynomials - """ - if kind not in (1,2): - raise ValueError('kind must be 1 or 2') - return -np.cos(np.pi*(np.arange(n)+0.5*kind)/(n+kind-1)); - - -def chebpoly(n, x=None, kind=1): - """ - Return Chebyshev polynomial of the first or second kind. - - These polynomials are orthogonal on the interval [-1,1], with - respect to the weight function w(x) = (1-x^2)^(-1/2+kind-1). - - chebpoly(n) returns the coefficients of the Chebychev polynomial of degree N. - chebpoly(n,x) returns the Chebychev polynomial of degree N evaluated in X. - - Parameters - ---------- - n : integer, scalar - degree of Chebychev polynomial. - x : array-like, optional - evaluation points - kind: 1 or 2, optional - kind of Chebychev polynomial (default 1) - - Returns - ------- - p : ndarray - polynomial coefficients if x is None. - Chebyshev polynomial evaluated at x otherwise - - Examples - -------- - >>> import numpy as np - >>> x = chebroot(3) - >>> np.abs(chebpoly(3,x))<1e-15 - array([ True, True, True], dtype=bool) - >>> chebpoly(3) - array([ 4., 0., -3., 0.]) - >>> x2 = chebroot(4,kind=2) - >>> np.abs(chebpoly(4,x2,kind=2))<1e-15 - array([ True, True, True, True], dtype=bool) - >>> chebpoly(4,kind=2) - array([ 16., 0., -12., 0., 1.]) - - - Reference - --------- - http://en.wikipedia.org/wiki/Chebyshev_polynomials - """ - if x is None: # Calculate coefficients. - if n == 0: - p = np.ones(1) - else: - p = np.round( pow(2,n-2+kind) * np.poly( chebroot(n,kind=kind) ) ) - p[1::2] = 0; - return p - else: # Evaluate polynomial in chebychev form - ck = np.zeros(n+1) - ck[n] = 1. - return _chebval(np.atleast_1d(x),ck,kind=kind) - -def chebfit(fun,n=10,a=-1,b=1,trace=False): - """ - Computes the Chebyshevs coefficients - - so that f(x) can be approximated by: - - n-1 - f(x) = sum ck*Tk(x) - k=0 - - where Tk is the k'th Chebyshev polynomial of the first kind. - - Parameters - ---------- - fun : callable - function to approximate - n : integer, scalar, optional - number of base points (abscissas). Default n=10 (maximum 50) - a,b : real, scalars, optional - integration limits - - Returns - ------- - ck : ndarray - polynomial coefficients in Chebychev form. - - Examples - -------- - Fit exp(x) - - >>> import pylab as pb - >>> a = 0; b = 2 - >>> ck = chebfit(pb.exp,7,a,b); - >>> x = pb.linspace(0,4); - >>> h=pb.plot(x,pb.exp(x),'r',x,chebval(x,ck,a,b),'g.') - >>> x1 = chebroot(9)*(b-a)/2+(b+a)/2 - >>> ck1 = chebfit(pb.exp(x1)) - >>> h=pb.plot(x,pb.exp(x),'r',x,chebval(x,ck1,a,b),'g.') - - >>> pb.close() - - See also - -------- - chebval - - Reference - --------- - http://en.wikipedia.org/wiki/Chebyshev_nodes - http://mathworld.wolfram.com/ChebyshevApproximationFormula.html - - W. Fraser (1965) - "A Survey of Methods of Computing Minimax and Near-Minimax Polynomial - Approximations for Functions of a Single Independent Variable" - Journal of the ACM (JACM), Vol. 12 , Issue 3, pp 295 - 314 - - William H. Press, Saul Teukolsky, - William T. Wetterling and Brian P. Flannery (1997) - "Numerical recipes in Fortran 77", Vol. 1, pp 184-194 - """ - - if (n>50): - warnings.warn('CHEBFIT should only be used for n<50') - - if hasattr(fun,'__call__'): - x = map_to_interval(chebroot(n),a,b) - f = fun(x); - if trace: - import pylab as plb - plb.plot(x,f,'+') - else: - f = fun - n = len(f) - #raise ValueError('Function must be callable!') - # N-1 - # c(k) = (2/N) sum w(n) f(n)*cos(pi*k*(2n+1)/(2N)), 0 <= k < N. - # n=0 - # - # w(0) = 0.5, w(n)=1 for n>0 - ck = dct(f[::-1])/n - ck[0] = ck[0]/2. - return ck - -def dct(x,n=None): - """ - Discrete Cosine Transform - - N-1 - y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. - n=0 - - Examples - -------- - >>> import numpy as np - >>> x = np.arange(5) - >>> np.abs(x-idct(dct(x)))<1e-14 - array([ True, True, True, True, True], dtype=bool) - >>> np.abs(x-dct(idct(x)))<1e-14 - array([ True, True, True, True, True], dtype=bool) - - Reference - --------- - http://en.wikipedia.org/wiki/Discrete_cosine_transform - http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ - """ - fft = np.fft.fft - x = np.atleast_1d(x) - - if n is None: - n = x.shape[-1] - - if x.shape[-1]0 - - Examples - -------- - >>> import numpy as np - >>> x = np.arange(5) - >>> np.abs(x-idct(dct(x)))<1e-14 - array([ True, True, True, True, True], dtype=bool) - >>> np.abs(x-dct(idct(x)))<1e-14 - array([ True, True, True, True, True], dtype=bool) - - Reference - --------- - http://en.wikipedia.org/wiki/Discrete_cosine_transform - http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ - """ - - ifft = np.fft.ifft - x = np.atleast_1d(x) - - if n is None: - n = x.shape[-1] - - w = np.exp(1j * np.arange(n) * np.pi/(2*n)) - - if x.shape[-1]>> import pylab as pb - >>> x = pb.linspace(-1,1) - >>> ck = pb.zeros(5); ck[-1]=1 - >>> h = pb.plot(x,chebval(x,ck),x,chebpoly(4,x),'.') - >>> pb.close() - - Fit exponential function: - >>> import pylab as pb - >>> ck = chebfit(pb.exp,7,0,2) - >>> x = pb.linspace(0,4); - >>> h=pb.plot(x,chebval(x,ck,0,2),'g',x,pb.exp(x)) - >>> pb.close() - - See also - -------- - chebfit - - References - ---------- - http://en.wikipedia.org/wiki/Clenshaw_algorithm - http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html - """ - - y = map_from_interval(np.atleast_1d(x),a,b) - if fill is None: - f = _chebval(y,ck,kind=kind) - else: - cond = (np.abs(y)<=1) - f = np.where(cond,0,fill) - if np.any(cond): - yk = np.extract(cond,y) - f[cond] = _chebval(yk,ck,kind=kind) - return f - - -def chebder(ck,a=-1,b=1): - """ - Differentiate Chebyshev polynomial - - Parameters - ---------- - ck : array-like - polynomial coefficients in Chebyshev form of function to differentiate - a,b : real, scalars - limits for polynomial(Default -1,1) - - Return - ------ - cder : ndarray - polynomial coefficients in Chebyshev form of the derivative - - Examples - -------- - - Fit exponential function: - >>> import pylab as pb - >>> ck = chebfit(pb.exp,7,0,2) - >>> x = pb.linspace(0,4) - >>> ck2 = chebder(ck,0,2); - >>> h=pb.plot(x,chebval(x,ck,0,2),'g',x,pb.exp(x),'r') - >>> pb.close() - - See also - -------- - chebint - chebfit - """ - - n = len(ck) - cder = np.zeros(n) - # n and n-1 are special cases. - # cder(n-1)=0; - cder[-2] = 2*(n-1)*ck[-1] - for j in range(n-3,-1,-1): - cder[j] = cder[j+2]+2*j*ck[j+1] - - return cder*2./(b-a) # Normalize to the interval b-a. - -def chebint(ck,a=-1,b=1): - """ - Integrate Chebyshef polynomial - - Parameters - ---------- - ck : array-like - polynomial coefficients in Chebyshev form of function to integrate. - a,b : real, scalars - limits for polynomial(Default -1,1) - - Return - ------ - cint : ndarray - polynomial coefficients in Chebyshev form of the integrated function - - Examples - -------- - Fit exponential function: - >>> import pylab as pb - >>> ck = chebfit(pb.exp,7,0,2) - >>> x = pb.linspace(0,4) - >>> ck2 = chebint(ck,0,2); - >>> h=pb.plot(x,chebval(x,ck,0,2),'g',x,pb.exp(x),'r.') - >>> pb.close() - - See also - -------- - chebder - chebfit - """ - - - n = len(ck) - - cint = np.zeros(n); - con = 0.25*(b-a); - - dif1 = np.diff(ck[::2]) - ix1= np.r_[1:n-1:2] - cint[ix1] = -(con*dif1)/(ix1-1) - if n>3: - dif2 = np.diff(ck[1::2]) - ix2=np.r_[2:n-1:2] - cint[ix2] = -(con*dif2)/(ix2-1) - - #% cint(n) is a special case - cint[-1] = (con*ck[n-2])/(n-1) - cint[0] = np.sum((-1)**np.r_[1:n]*cint[1::]) # Set constant of integration - - - return cint - -class Cheb1d(object): - coeffs = None - order = None - a = None - b = None - def __init__(self,ck,a=-1,b=1): - if isinstance(ck, poly1d): - for key in ck.__dict__.keys(): - self.__dict__[key] = ck.__dict__[key] - return - cki = np.trim_zeros(np.atleast_1d(ck),'b') - if len(cki.shape) > 1: - raise ValueError, "Polynomial must be 1d only." - self.__dict__['coeffs'] = cki - self.__dict__['order'] = len(cki) - 1 - self.__dict__['a'] = a - self.__dict__['b'] = b - - - def __call__(self,x): - return chebval(x,self.coeffs,self.a,self.b) - - def __array__(self, t=None): - if t: - return np.asarray(self.coeffs, t) - else: - return np.asarray(self.coeffs) - - def __repr__(self): - vals = repr(self.coeffs) - vals = vals[6:-1] - return "Cheb1d(%s)" % vals - - def __len__(self): - return self.order - - def __str__(self): - pass - def __neg__(self): - return Cheb1d(-self.coeffs,self.a,self.b) - - def __pos__(self): - return self - - - def __add__(self, other): - other = poly1d(other) - return poly1d(polyadd(self.coeffs, other.coeffs)) - - def __radd__(self, other): - other = poly1d(other) - return poly1d(polyadd(self.coeffs, other.coeffs)) - - def __sub__(self, other): - other = poly1d(other) - return poly1d(polysub(self.coeffs, other.coeffs)) - - def __rsub__(self, other): - other = poly1d(other) - return poly1d(polysub(other.coeffs, self.coeffs)) - - def __eq__(self, other): - return np.alltrue(self.coeffs == other.coeffs) - - def __ne__(self, other): - return np.any(self.coeffs != other.coeffs) or (self.a!=other.a) or (self.b !=other.b) - - def __setattr__(self, key, val): - raise ValueError, "Attributes cannot be changed this way." - - def __getattr__(self, key): - if key in ['c','coef','coefficients']: - return self.coeffs - elif key in ['o']: - return self.order - elif key in ['a']: - return self.a - elif key in ['b']: - return self.b - else: - try: - return self.__dict__[key] - except KeyError: - raise AttributeError("'%s' has no attribute '%s'" % (self.__class__, key)) - def __getitem__(self, val): - if val > self.order: - return 0 - if val < 0: - return 0 - return self.coeffs[val] - - def __setitem__(self, key, val): - ind = self.order - key - if key < 0: - raise ValueError, "Does not support negative powers." - if key > self.order: - zr = NX.zeros(key-self.order, self.coeffs.dtype) - self.__dict__['coeffs'] = NX.concatenate((self.coeffs,zr)) - self.__dict__['order'] = key - self.__dict__['coeffs'][key] = val - return - - def __iter__(self): - return iter(self.coeffs) - def integ(self, m=1): - """ - Return an antiderivative (indefinite integral) of this polynomial. - - Refer to `chebint` for full documentation. - - See Also - -------- - chebint : equivalent function - - """ - return Cheb1d(chebint(self.coeffs, self.a,self.b)) - - def deriv(self, m=1): - """ - Return a derivative of this polynomial. - - Refer to `chebder` for full documentation. - - See Also - -------- - chebder : equivalent function - - """ - return Cheb1d(chebder(self.coeffs,self.a,self.b)) -def test_doctstrings(): - if False: #True: # - x = np.arange(4) - dx = dct(x) - idx = idct(dx) - import pylab as plb - a = 0; - b = 2; - ck = chebfit(np.exp,6,a,b); - t = chebval(0,ck,a,b) - x=np.linspace(0,2,6); - plb.plot(x,np.exp(x),'r', x,chebval(x,ck,a,b),'g.') - #x1 = chebroot(9).'*(b-a)/2+(b+a)/2 ; - #ck1 =chebfit([x1 exp(x1)],9,a,b); - #plot(x,exp(x),'r'), hold on - #plot(x,chebval(x,ck1,a,b),'g'), hold off - - - t = poly2hstr([1,1,2]) - py = [1, 0] - px = polyshift(py,0,5); - t1=polyval(px,[0, 2.5, 5]) #% This is the same as the line below - t2=polyval(py,[-1, 0, 1 ]) - - px = [1, 0] - py = polyishift(px,0,5); - t1 = polyval(px,[0, 2.5, 5]) #% This is the same as the line below - t2 = polyval(py,[-1, 0, 1 ]) - print(t1,t2) - else: - import doctest - doctest.testmod() -if __name__== '__main__': - test_doctstrings() diff --git a/wafo/powerpoint.py b/wafo/powerpoint.py index d39b487..e012fce 100644 --- a/wafo/powerpoint.py +++ b/wafo/powerpoint.py @@ -3,12 +3,12 @@ Created on 15. des. 2009 @author: pab ''' -#import os -#import sys -#import win32com -#from win32com.client.selecttlb import EnumTlbs -#typelib_mso = None -#typelib_msppt = None +# import os +# import sys +# import win32com +# from win32com.client.selecttlb import EnumTlbs +# typelib_mso = None +# typelib_msppt = None # for typelib in EnumTlbs(): # d = typelib.desc.split(' ') # if d[0] == 'Microsoft' and d[1] == 'Office' and d[3] == 'Object' \ @@ -19,7 +19,7 @@ Created on 15. des. 2009 # typelib_msppt = typelib # if hasattr(sys, 'frozen'): # If we're an .exe file # win32com.__gen_path__ = os.path.dirname(sys.executable) -## win32com.__gen_path__ = os.environ['TEMP'] +# # win32com.__gen_path__ = os.environ['TEMP'] # if win32com.client.gencache.is_readonly: # win32com.client.gencache.is_readonly = False # win32com.client.gencache.Rebuild() @@ -49,7 +49,7 @@ class Powerpoint(object): def __init__(self, file_name=''): self.application = win32com.client.Dispatch("Powerpoint.Application") - #self.application.Visible = True + # self.application.Visible = True self._visible = self.application.Visible if file_name: self.presentation = self.application.Presentations.Open(file_name) @@ -129,7 +129,7 @@ class Powerpoint(object): titlerange.Font.Size = self.title_size if texts != '' and texts != ['']: - #textrange = slide.Shapes(textbox_id).TextFrame.TextRange + # textrange = slide.Shapes(textbox_id).TextFrame.TextRange self._add_text(slide, textbox_id, texts, maxlevel) if image_file != '' and image_file != ['']: @@ -213,7 +213,6 @@ def test_powerpoint(): # Make powerpoint ppt = Powerpoint() - # time. ppt.footer = 'This is the footer' ppt.add_title_slide('Title', 'Per A.') ppt.add_slide(title='alsfkasldk', texts='asdflaf', notes='asdfas') @@ -231,7 +230,7 @@ def make_ppt(): # title.TextFrame.TextRange.Text = 'Overskrift' title_id, textbox_id = 1, 2 slide1.Shapes(title_id).TextFrame.TextRange.Text = 'Overskrift' - #slide1.Shapes(title_id).TextFrame.Width = 190 + # slide1.Shapes(title_id).TextFrame.Width = 190 slide1.Shapes(textbox_id).TextFrame.TextRange.InsertAfter('Test') unused_tr = slide1.Shapes(textbox_id).TextFrame.TextRange.InsertAfter('\r') @@ -242,12 +241,12 @@ def make_ppt(): tr.IndentLevel = 2 tr1 = slide1.Shapes(textbox_id).TextFrame.TextRange.InsertAfter('test3') tr1.IndentLevel = 3 - #slide1.Shapes(textbox_id).TextFrame.TextRange.Text = 'Test \r test2' + # slide1.Shapes(textbox_id).TextFrame.TextRange.Text = 'Test \r test2' # textbox = slide1.Shapes.AddTextBox(Type=msoTextOrientationHorizontal, # Left=30, Top=100, Width=190, Height=400) # textbox.TextFrame.TextRange.Text = 'Test \r test2' - #picbox = slide1.Shapes(picb_id) + # picbox = slide1.Shapes(picb_id) filename = r'c:\temp\data1_report1_and_2_Tr120_1.png' slide1.Shapes.AddPicture(FileName=filename, LinkToFile=False, @@ -280,13 +279,8 @@ def make_ppt(): # application.Quit() def rename_ppt(): root = r'C:/pab/tsm_opeval/analysis_tsmps_aco_v2008b/plots' -# root = r'C:/pab/tsm_opeval/analysis_tsmps_mag_v2008b/plots' -# root = r'C:/pab/tsm_opeval/analysis_tsmps_mag_v2010a/plots' -# root = r'C:/pab/tsm_opeval/analysis_tsmps_aco_v2010a/plots' - #filename = r'mag_sweep_best_tsmps_ship_eff0-10.ppt' filenames = os.listdir(root) prefix = 'TSMPSv2008b_' - #prefix = 'TSMPSv2010a_' for filename in filenames: if filename.endswith('.ppt'): try: @@ -300,13 +294,8 @@ def rename_ppt(): def load_file_into_ppt(): root = r'C:/pab/tsm_opeval/analysis_tsmps_aco_v2008b/plots' -# root = r'C:/pab/tsm_opeval/analysis_tsmps_mag_v2008b/plots' -# root = r'C:/pab/tsm_opeval/analysis_tsmps_mag_v2010a/plots' -# root = r'C:/pab/tsm_opeval/analysis_tsmps_aco_v2010a/plots' - #filename = r'mag_sweep_best_tsmps_ship_eff0-10.ppt' filenames = os.listdir(root) prefix = 'TSMPSv2008b_' - #prefix = 'TSMPSv2010a_' for filename in filenames: if filename.startswith(prefix) and filename.endswith('.ppt'): try: diff --git a/wafo/setup.py b/wafo/setup.py index e0a51bb..a535100 100644 --- a/wafo/setup.py +++ b/wafo/setup.py @@ -5,6 +5,7 @@ Created on Sun Oct 25 14:55:34 2015 @author: dave """ + def configuration(parent_package='', top_path=None): from numpy.distutils.misc_util import Configuration config = Configuration('wafo', parent_package, top_path) diff --git a/wafo/tests/test_gaussian.py b/wafo/tests/test_gaussian.py index ae0fff8..1c7e795 100644 --- a/wafo/tests/test_gaussian.py +++ b/wafo/tests/test_gaussian.py @@ -3,12 +3,11 @@ Created on 17. juli 2010 @author: pab ''' -import numpy as np # @UnusedImport -from numpy import pi, inf # @UnusedImport +import numpy as np +from numpy import pi, inf from numpy.testing import assert_array_almost_equal from wafo.gaussian import (Rind, prbnormtndpc, prbnormndpc, prbnormnd, cdfnorm2d, prbnorm2d) -from numpy.ma.testutils import assert_array_almost_equal def test_rind(): @@ -31,7 +30,7 @@ def test_rind(): A = np.repeat(Blo, n) B = np.repeat(Bup, n) # Integration limits - E1, err1, terr1 = rind(np.triu(Sc), m, A, B) # same as E0 + E1, _err1, _terr1 = rind(np.triu(Sc), m, A, B) # same as E0 assert(np.abs(E1 - Et) < err0 + terr0) t = 'E1 = %2.5f' % E1 @@ -59,8 +58,10 @@ def test_rind(): Bup2 = np.inf indI2 = [-1, 1] rind2 = Rind(method=1) - g2 = lambda x: ( - x * (np.pi / 2 + np.arcsin(x)) + np.sqrt(1 - x**2)) / (2 * np.pi) + + def g2(x): + return (x * (np.pi / 2 + np.arcsin(x)) + + np.sqrt(1 - x**2)) / (2 * np.pi) assert_array_almost_equal(g2(rho2), 0.24137214191774381) # exact value E3, err3, terr3 = rind(Sc2, m2, Blo2, Bup2, indI2, nt=0) @@ -84,18 +85,21 @@ def test_prbnormtndpc(): rho2 = np.random.rand(2) a2 = np.zeros(2) b2 = np.repeat(np.inf, 2) - [val2, err2, ift2] = prbnormtndpc(rho2, a2, b2) - g2 = lambda x: 0.25 + np.arcsin(x[0] * x[1]) / (2 * pi) - E2 = g2(rho2) # % exact value + val2, err2, _ift2 = prbnormtndpc(rho2, a2, b2) + + def g2(x): + return 0.25 + np.arcsin(x[0] * x[1]) / (2 * pi) + E2 = g2(rho2) # exact value assert(np.abs(E2 - val2) < err2) rho3 = np.random.rand(3) a3 = np.zeros(3) b3 = np.repeat(inf, 3) - [val3, err3, ift3] = prbnormtndpc(rho3, a3, b3) - g3 = lambda x: 0.5 - \ - sum(np.sort( - np.arccos([x[0] * x[1], x[0] * x[2], x[1] * x[2]]))) / (4 * pi) + val3, err3, _ift3 = prbnormtndpc(rho3, a3, b3) + + def g3(x): + return 0.5 - sum(np.sort(np.arccos([x[0] * x[1], x[0] * x[2], + x[1] * x[2]]))) / (4 * pi) E3 = g3(rho3) # Exact value assert(np.abs(E3 - val3) < err3) @@ -105,17 +109,21 @@ def test_prbnormndpc(): rho2 = np.random.rand(2) a2 = np.zeros(2) b2 = np.repeat(np.inf, 2) - [val2, err2, ift2] = prbnormndpc(rho2, a2, b2) - g2 = lambda x: 0.25 + np.arcsin(x[0] * x[1]) / (2 * pi) - E2 = g2(rho2) # % exact value + val2, err2, _ift2 = prbnormndpc(rho2, a2, b2) + + def g2(x): + return 0.25 + np.arcsin(x[0] * x[1]) / (2 * pi) + E2 = g2(rho2) # exact value assert(np.abs(E2 - val2) < err2) rho3 = np.random.rand(3) a3 = np.zeros(3) b3 = np.repeat(inf, 3) - [val3, err3, ift3] = prbnormndpc(rho3, a3, b3) - g3 = lambda x: 0.5 - sum(np.sort(np.arccos([x[0] * x[1], x[0] * x[2], - x[1] * x[2]]))) / (4 * pi) + val3, err3, _ift3 = prbnormndpc(rho3, a3, b3) + + def g3(x): + return 0.5 - sum(np.sort(np.arccos([x[0] * x[1], x[0] * x[2], + x[1] * x[2]]))) / (4 * pi) E3 = g3(rho3) # Exact value assert(np.abs(E3 - val3) < err3) @@ -153,7 +161,7 @@ def test_prbnorm2d(): a = [-1, -2] b = [1, 1] r = 0.3 - assert_array_almost_equal(prbnorm2d(a,b,r), [ 0.56659121]) + assert_array_almost_equal(prbnorm2d(a, b, r), 0.56659121) if __name__ == '__main__': import doctest diff --git a/wafo/tests/test_misc.py b/wafo/tests/test_misc.py index adda210..8b7713d 100644 --- a/wafo/tests/test_misc.py +++ b/wafo/tests/test_misc.py @@ -1,6 +1,7 @@ from numpy.testing import (run_module_suite, assert_equal, assert_almost_equal, - assert_array_equal, assert_array_almost_equal) + assert_array_equal, assert_array_almost_equal, + TestCase, assert_, assert_raises,) import numpy as np from numpy import array, cos, exp, linspace, pi, sin, diff, arange, ones @@ -10,7 +11,10 @@ from wafo.misc import (JITImport, Bunch, detrendma, DotDict, findcross, ecross, findoutliers, common_shape, argsreduce, stirlerr, getshipchar, betaloge, hygfz, gravity, nextpow2, discretize, polar2cart, - cart2polar, tranproc) + cart2polar, tranproc, + rotation_matrix, rotate_2d, spaceline, + args_flat, sub2index, index2sub, piecewise, + parse_kwargs) def test_JITImport(): @@ -342,6 +346,17 @@ def test_stirlerr(): 0.02767793, 0.02079067])) +def test_parse_kwargs(): + opt = dict(arg1=1, arg2=3) + opt = parse_kwargs(opt, arg1=5) + assert(opt['arg1'] == 5) + assert(opt['arg2'] == 3) + opt2 = dict(arg3=15) + + opt = parse_kwargs(opt, **opt2) + assert('arg3' not in opt) + + def test_getshipchar(): sc = getshipchar(10, 'service_speed') true_sc = dict(beam=29, @@ -380,7 +395,7 @@ def test_nextpow2(): def test_discretize(): - x, y = discretize(np.cos, 0, np.pi, tol=0.0051) + x, y = discretize(np.cos, 0, np.pi, tol=0.05) assert_array_almost_equal( x, np.array( @@ -470,5 +485,241 @@ def test_tranproc(): 0.86643821, 0.83096482])) +class TestPiecewise(TestCase): + def test_condition_is_single_bool_list(self): + assert_raises(ValueError, piecewise, [0, 0], [True, False], [1]) + + def test_condition_is_list_of_single_bool_list(self): + x = piecewise([0, 0], [[True, False]], [1]) + assert_array_equal(x, [1, 0]) + + def test_conditions_is_list_of_single_bool_array(self): + x = piecewise([0, 0], [np.array([True, False])], [1]) + assert_array_equal(x, [1, 0]) + + def test_condition_is_single_int_array(self): + assert_raises(ValueError, piecewise, [0, 0], np.array([1, 0]), [1]) + + def test_condition_is_list_of_single_int_array(self): + x = piecewise([0, 0], [np.array([1, 0])], [1]) + assert_array_equal(x, [1, 0]) + + def test_simple(self): + x = piecewise([0, 0], [[False, True]], [lambda x:-1]) + assert_array_equal(x, [0, -1]) + + x = piecewise([1, 2], [[True, False], [False, True]], [3, 4]) + assert_array_equal(x, [3, 4]) + + def test_default(self): + # No value specified for x[1], should be 0 + x = piecewise([1, 2], [[True, False]], [2]) + assert_array_equal(x, [2, 0]) + + # Should set x[1] to 3 + x = piecewise([1, 2], [[True, False]], [2, 3]) + assert_array_equal(x, [2, 3]) + + def test_0d(self): + x = np.array(3) + y = piecewise(x, [x > 3], [4, 0]) + assert_(y.ndim == 0) + assert_(y == 0) + + x = 5 + y = piecewise(x, [[True], [False]], [1, 0]) + assert_(y.ndim == 0) + assert_(y == 1) + + def test_abs_function(self): + x = np.linspace(-2.5, 2.5, 6) + vals = piecewise((x,), [x < 0, x >= 0], [lambda x: -x, lambda x: x]) + assert_array_equal(vals, + [2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) + + def test_abs_function_with_scalar(self): + x = np.array(-2.5) + vals = piecewise((x,), [x < 0, x >= 0], [lambda x: -x, lambda x: x]) + assert_(vals == 2.5) + + def test_otherwise_condition(self): + x = np.linspace(-2.5, 2.5, 6) + vals = piecewise((x,), [x < 0, ], [lambda x: -x, lambda x: x]) + assert_array_equal(vals, [2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) + + def test_passing_further_args_to_fun(self): + def fun0(x, y, scale=1.): + return -x*y/scale + + def fun1(x, y, scale=1.): + return x*y/scale + x = np.linspace(-2.5, 2.5, 6) + vals = piecewise((x,), [x < 0, ], [fun0, fun1], args=(2.,), scale=2.) + assert_array_equal(vals, [2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) + + def test_step_function(self): + x = np.linspace(-2.5, 2.5, 6) + vals = piecewise(x, [x < 0, x >= 0], [-1, 1]) + assert_array_equal(vals, [-1., -1., -1., 1., 1., 1.]) + + def test_step_function_with_scalar(self): + x = 1 + vals = piecewise(x, [x < 0, x >= 0], [-1, 1]) + assert_(vals == 1) + + def test_function_with_two_args(self): + x = np.linspace(-2, 2, 5) + X, Y = np.meshgrid(x, x) + vals = piecewise( + (X, Y), [X * Y < 0, ], [lambda x, y: -x * y, lambda x, y: x * y]) + assert_array_equal(vals, [[4., 2., -0., 2., 4.], + [2., 1., -0., 1., 2.], + [-0., -0., 0., 0., 0.], + [2., 1., 0., 1., 2.], + [4., 2., 0., 2., 4.]]) + + def test_fill_value_and_function_with_two_args(self): + x = np.linspace(-2, 2, 5) + X, Y = np.meshgrid(x, x) + vals = piecewise((X, Y), [X * Y < -0.5, X * Y > 0.5], + [lambda x, y: -x * y, lambda x, y: x * y], + fill_value=np.nan) + nan = np.nan + assert_array_equal(vals, [[4., 2., nan, 2., 4.], + [2., 1., nan, 1., 2.], + [nan, nan, nan, nan, nan], + [2., 1., nan, 1., 2.], + [4., 2., nan, 2., 4.]]) + + def test_fill_value2_and_function_with_two_args(self): + x = np.linspace(-2, 2, 5) + X, Y = np.meshgrid(x, x) + vals = piecewise((X, Y), [X * Y < -0.5, X * Y > 0.5], + [lambda x, y: -x * y, lambda x, y: x * y, np.nan]) + nan = np.nan + assert_array_equal(vals, [[4., 2., nan, 2., 4.], + [2., 1., nan, 1., 2.], + [nan, nan, nan, nan, nan], + [2., 1., nan, 1., 2.], + [4., 2., nan, 2., 4.]]) + + +class TestRotationMatrix(TestCase): + + def test_h0_p0_r0(self): + vals = rotation_matrix(heading=0, pitch=0, roll=0).tolist() + truevals = [[1., 0., 0.], + [0., 1., 0.], + [0., 0., 1.]] + self.assertListEqual(vals, truevals) + + def test_h180_p0_r0(self): + vals = rotation_matrix(heading=180, pitch=0, roll=0).tolist() + truevals = [[-1.0, -1.2246467991473532e-16, 0.0], + [1.2246467991473532e-16, -1.0, 0.0], + [-0.0, 0.0, 1.0]] + self.assertListEqual(vals, truevals) + + def test_h0_p180_r0(self): + vals = rotation_matrix(heading=0, pitch=180, roll=0).tolist() + truevals = [[-1.0, 0.0, 1.2246467991473532e-16], + [-0.0, 1.0, 0.0], + [-1.2246467991473532e-16, -0.0, -1.0]] + self.assertListEqual(vals, truevals) + + def test_h0_p0_r180(self): + vals = rotation_matrix(heading=0, pitch=180, roll=0).tolist() + truevals = [[-1.0, 0.0, 1.2246467991473532e-16], + [-0.0, 1.0, 0.0], + [-1.2246467991473532e-16, -0.0, -1.0]] + self.assertListEqual(vals, truevals) + + +class TestRotate2d(TestCase): + + def test_rotate_0deg(self): + vals = list(rotate_2d(x=1, y=0, angle_deg=0)) + truevals = [1.0, 0.0] + self.assertListEqual(vals, truevals) + + def test_rotate_90deg(self): + vals = list(rotate_2d(x=1, y=0, angle_deg=90)) + truevals = [6.123233995736766e-17, 1.0] + self.assertListEqual(vals, truevals) + + def test_rotate_180deg(self): + vals = list(rotate_2d(x=1, y=0, angle_deg=180)) + truevals = [-1.0, 1.2246467991473532e-16] + self.assertListEqual(vals, truevals) + + def test_rotate_360deg(self): + vals = list(rotate_2d(x=1, y=0, angle_deg=360)) + truevals = [1.0, -2.4492935982947064e-16] + self.assertListEqual(vals, truevals) + + +class TestSpaceLine(TestCase): + + def test_space_line(self): + vals = spaceline((2, 0, 0), (3, 0, 0), num=5).tolist() + truevals = [[2., 0., 0.], + [2.25, 0., 0.], + [2.5, 0., 0.], + [2.75, 0., 0.], + [3., 0., 0.]] + self.assertListEqual(vals, truevals) + + +class TestArgsFlat(TestCase): + + def test_1_vector_and_2_scalar_args(self): + x = [1, 2, 3] + pos, c_shape = args_flat(x, 2, 3) + truepos = [[1, 2, 3], + [2, 2, 3], + [3, 2, 3]] + truec_shape = [3, ] + self.assertListEqual(pos.tolist(), truepos) + self.assertListEqual(list(c_shape), truec_shape) + + def test_1_vector_args(self): + pos1, c_shape1 = args_flat([1, 2, 3]) + truepos1 = [[1, 2, 3]] + truec_shape1 = None + self.assertListEqual(pos1.tolist(), truepos1) + self.assertIs(c_shape1, truec_shape1) + + def test_3_scalar_args(self): + pos1, c_shape1 = args_flat(1, 2, 3) + truepos1 = [[1, 2, 3]] + truec_shape1 = [] + self.assertListEqual(pos1.tolist(), truepos1) + self.assertListEqual(list(c_shape1), truec_shape1) + + def test_3_scalar_args_version2(self): + pos1, c_shape1 = args_flat([1], 2, 3) + truepos1 = [[1, 2, 3]] + truec_shape1 = [1, ] + self.assertListEqual(pos1.tolist(), truepos1) + self.assertListEqual(list(c_shape1), truec_shape1) + + +class TestSub2index2Sub(TestCase): + + def test_sub2index_and_index2sub(self): + shape = (3, 3, 4) + a = np.arange(np.prod(shape)).reshape(shape) + trueval = a[1, 2, 3] + order = 'C' + i = sub2index(shape, 1, 2, 3, order=order) + self.assertEquals(i, 23) + + val = a.ravel(order)[i] + self.assertEquals(val, trueval) + + sub = index2sub(shape, i, order=order) + for j, true_sub_j in enumerate([1, 2, 3]): + self.assertEquals(sub[j].tolist(), true_sub_j) + if __name__ == '__main__': run_module_suite() diff --git a/wafo/tests/test_numpy_utils.py b/wafo/tests/test_numpy_utils.py index c088c09..f7185a4 100644 --- a/wafo/tests/test_numpy_utils.py +++ b/wafo/tests/test_numpy_utils.py @@ -6,9 +6,9 @@ from numpy.testing import ( # ) import unittest as local_unittest +import numpy as np from wafo.numpy_utils import (rotation_matrix, rotate_2d, spaceline, args_flat, sub2index, index2sub, piecewise) -import numpy as np class TestPiecewise(TestCase):