From e2793d38c7ebc399c64da612e5fab958e7d3f4c3 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Sat, 17 Jul 2010 11:44:32 +0000 Subject: [PATCH] Added more test to test/test_misc.py and test/test_gaussian.py --- pywafo/src/wafo/gaussian.py | 337 ++++++++++++++++++++++++- pywafo/src/wafo/misc.py | 77 +++--- pywafo/src/wafo/stats/core.py | 4 +- pywafo/src/wafo/stats/estimation.py | 1 + pywafo/src/wafo/test/test_gaussian.py | 111 ++++++++ pywafo/src/wafo/test/test_misc.py | 351 +++++++++++++++++++++++++- 6 files changed, 822 insertions(+), 59 deletions(-) create mode 100644 pywafo/src/wafo/test/test_gaussian.py diff --git a/pywafo/src/wafo/gaussian.py b/pywafo/src/wafo/gaussian.py index 44d6cab..62518f9 100644 --- a/pywafo/src/wafo/gaussian.py +++ b/pywafo/src/wafo/gaussian.py @@ -1,9 +1,12 @@ +import warnings import numpy as np from numpy import (r_, minimum, maximum, atleast_1d, atleast_2d, mod, zeros, #@UnresolvedImport - ones, floor, random, eye, nonzero, repeat, sqrt, inf, diag, triu) #@UnresolvedImport + ones, floor, random, eye, nonzero, repeat, sqrt, exp, inf, diag, triu) #@UnresolvedImport from scipy.special import ndtri as invnorm +from scipy.special import ndtr as cdfnorm import wafo.rindmod as rindmod - +import wafo.mvnprdmod as mvnprdmod +from wafo import mvn class Rind(object): ''' @@ -68,7 +71,7 @@ class Rind(object): >>> m = np.zeros(n); rho = 0.3; >>> Sc =(np.ones((n,n))-np.eye(n))*rho+np.eye(n) >>> rind = Rind() - >>> E0 = rind(Sc,m,Blo,Bup,indI) # exact prob. 0.001946 + >>> E0, err0, terr0 = rind(Sc,m,Blo,Bup,indI) # exact prob. 0.001946 >>> A = np.repeat(Blo,n); B = np.repeat(Bup,n) # Integration limits >>> E1 = rind(np.triu(Sc),m,A,B) #same as E0 @@ -358,6 +361,334 @@ def test_rind(): Blo[0, ind] = maximum(Blo[0, ind], -infinity * dev[indI[ind + 1]]) E3 = rind(Sc, m, Blo, Bup, indI, xc, nt=1) +def cdflomax(x,alpha,m0): + ''' + Return CDF for local maxima for a zero-mean Gaussian process + + Parameters + ---------- + x : array-like + evaluation points + alpha, m0 : real scalars + irregularity factor and zero-order spectral moment (variance of the + process), respectively. + + Returns + ------- + prb : ndarray + distribution function evaluated at x + + Notes + ----- + The cdf is calculated from an explicit expression involving the + standard-normal cdf. This relation is sometimes written as a convolution + + M = sqrt(m0)*( sqrt(1-a^2)*Z + a*R ) + + where M denotes local maximum, Z is a standard normal r.v., + R is a standard Rayleigh r.v., and "=" means equality in distribution. + + Note that all local maxima of the process are considered, not + only crests of waves. + + Example + ------- + >>> pylab + >>> import wafo.spectrum.models as wsm + >>> import wafo.objects as wo + >>> import wafo.stats as ws + >>> S = wsm.Jonswap(Hm0=10).tospecdata(); + >>> xs = S.sim(10000) + >>> ts = wo.mat2timeseries(xs) + >>> tp = ts.turning_points() + >>> mM = tp.cycle_pairs() + >>> m0 = S.moment(1)[0] + >>> alpha = S.characteristic('alpha')[0] + >>> x = linspace(-10,10,200); + >>> mcdf = ws.edf(mM.data) + >>> mcdf.plot(), pylab.plot(x,cdflomax(x,alpha,m0)) + + See also + -------- + spec2mom, spec2bw + ''' + c1 = 1.0/(sqrt(1-alpha**2))*x/sqrt(m0) + c2 = alpha*c1 + return cdfnorm(c1)-alpha*exp(-x**2/2/m0)*cdfnorm(c2) + +def prbnormtndpc(rho,a,b,D=None,df=0,abseps=1e-4,IERC=0,HNC=0.24): + ''' + Return Multivariate normal or T probability with product correlation structure. + + Parameters + ---------- + rho : array-like + vector of coefficients defining the correlation coefficient by: + correlation(I,J) = rho[i]*rho[j]) for J!=I + where -1 < rho[i] < 1 + a,b : array-like + vector of lower and upper integration limits, respectively. + Note: any values greater the 37 in magnitude, are considered as infinite values. + D : array-like + vector of means (default zeros(size(rho))) + df = Degrees of freedom, NDF<=0 gives normal probabilities (default) + abseps = absolute error tolerance. (default 1e-4) + IERC = 1 if strict error control based on fourth derivative + 0 if intuitive error control based on halving the intervals (default) + HNC = start interval width of simpson rule (default 0.24) + + Returns + ------- + value = estimated value for the integral + bound = bound on the error of the approximation + inform = INTEGER, termination status parameter: + 0, if normal completion with ERROR < EPS; + 1, if N > 1000 or N < 1. + 2, IF any abs(rho)>=1 + 4, if ANY(b(I)<=A(i)) + 5, if number of terms computed exceeds maximum number of evaluation points + 6, if fault accurs in normal subroutines + 7, if subintervals are too narrow or too many + 8, if bounds exceeds abseps + + PRBNORMTNDPC calculates multivariate normal or student T probability + with product correlation structure for rectangular regions. + The accuracy is as best around single precision, i.e., about 1e-7. + + Example: + -------- + >>> 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 + >>> np.abs(E2-val2)>> 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) + >>> E3 = g3(rho3) # Exact value + >>> np.abs(E3-val3)>> 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 + >>> np.abs(E2-val2)>> 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) + >>> E3 = g3(rho3) # Exact value + >>> np.abs(E3-val3)0: + warnings.warn('Abnormal termination ier = %d\n\n%s' % (ier,_ERRORMESSAGE[ier])) + return val, err, ier + +_ERRORMESSAGE = {} +_ERRORMESSAGE[0] = '' +_ERRORMESSAGE[1] =''' + Maximum number of subdivisions allowed has been achieved. one can allow + more subdivisions by increasing the value of limit (and taking the + according dimension adjustments into account). however, if this yields + no improvement it is advised to analyze the integrand in order to + determine the integration difficulties. if the position of a local + difficulty can be determined (i.e. singularity discontinuity within + the interval), it should be supplied to the routine as an element of + the vector points. If necessary an appropriate special-purpose integrator + must be used, which is designed for handling the type of difficulty involved. + ''' +_ERRORMESSAGE[2] =''' + the occurrence of roundoff error is detected, which prevents the requested + tolerance from being achieved. The error may be under-estimated.''' + +_ERRORMESSAGE[3] =''' + Extremely bad integrand behaviour occurs at some points of the integration interval.''' +_ERRORMESSAGE[4] =''' + The algorithm does not converge. Roundoff error is detected in the extrapolation table. + It is presumed that the requested tolerance cannot be achieved, and that + the returned result is the best which can be obtained.''' +_ERRORMESSAGE[5] =''' + The integral is probably divergent, or slowly convergent. + It must be noted that divergence can occur with any other value of ier>0.''' +_ERRORMESSAGE[6] ='''the input is invalid because: + 1) npts2 < 2 + 2) break points are specified outside the integration range + 3) (epsabs<=0 and epsrel EPS and MAXPTS + function vaules used; increase MAXPTS to + decrease ERROR; + if INFORM = 2, N > NMAX or N < 1. where NMAX depends on the + integration method + Example:% Compute the probability that X1<0,X2<0,X3<0,X4<0,X5<0, + % Xi are zero-mean Gaussian variables with variances one + % and correlations Cov(X(i),X(j))=0.3: + % indI=[0 5], and barriers B_lo=[-inf 0], B_lo=[0 inf] + % gives H_lo = [-inf -inf -inf -inf -inf] H_lo = [0 0 0 0 0] + + N = 5; rho=0.3; NIT=3; Nt=N; indI=[0 N]; + B_lo=-10; B_up=0; m=1.2*ones(N,1); + Sc=(ones(N)-eye(N))*rho+eye(N); + E = rind(Sc,m,B_lo,B_up,indI,[],Nt) % exact prob. 0.00195 + A = [-inf -inf -inf -inf -inf], + B = [0 0 0 0 0]-m' + [val,err,inform] = prbnormnd(Sc,A,B); + + See also + -------- + prbnormndpc, rind + ''' + + + m,n = correl.shape + Na = len(a) + Nb = len(b) + if (m!=n or m!=Na or m!=Nb): + raise ValueError('Size of input is inconsistent!') + + if maxpts is None: + maxpts = 1000*n + + maxpts = max(round(maxpts),10*n); + +# % array of correlation coefficients; the correlation +# % coefficient in row I column J of the correlation matrix +# % should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I. +# % The correlation matrix must be positive semidefinite. + + D = np.diag(correl) + if (any(D!=1)): + raise ValueError('This is not a correlation matrix') + + + # Make sure integration limits are finite + A = np.clip(a,-100,100) + B = np.clip(b,-100, 100) + + #L = correl((triu(ones(m),1)~=0)); % return only off diagonal elements + return mvn.mvnun(A,B, correl,maxpts, abseps, releps) + + #CALL the mexroutine +# t0 = clock; +# if ((method==0) && (n<=100)), +# %NMAX = 100 +# [value, err,inform] = mexmvnprb(L,A,B,abseps,releps,maxpts); +# elseif ( (method<0) || ((method<=0) && (n>100)) ), +# % NMAX = 500 +# [value, err,inform] = mexmvnprb2(L,A,B,abseps,releps,maxpts); +# else +# [value, err,inform] = mexGenzMvnPrb(L,A,B,abseps,releps,maxpts,method); +# end +# exTime = etime(clock,t0); +# if __name__ == '__main__': if False: #True: # diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index 3a70072..db5169f 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -6,45 +6,11 @@ from __future__ import division import sys import numpy as np -from numpy import abs -from numpy import amax -from numpy import any -from numpy import arange -from numpy import arctan2 -from numpy import array -from numpy import asarray -from numpy import atleast_1d -from numpy import broadcast_arrays -from numpy import ceil -from numpy import cos -from numpy import diff -from numpy import empty_like -from numpy import exp -from numpy import extract -from numpy import finfo -from numpy import floor -from numpy import frexp -from numpy import hstack -from numpy import hypot -from numpy import inf -from numpy import interp -from numpy import isnan -from numpy import isscalar -from numpy import linspace -from numpy import log -from numpy import logical_and -from numpy import mod -from numpy import nonzero -from numpy import ones -from numpy import pi -from numpy import r_ -from numpy import sign -from numpy import sin -from numpy import sqrt -from numpy import unique -from numpy import vstack -from numpy import where -from numpy import zeros +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, mod, diff, empty_like, + finfo, inf, pi, interp, isnan, isscalar, zeros, ones, + r_, sign, unique, hstack, vstack, nonzero, where, extract) from scipy.special import gammaln import types import warnings @@ -99,7 +65,12 @@ class DotDict(dict): class Bunch(object): ''' Implement keyword argument initialization of class - + + Example + ------- + >>> d = Bunch(test1=1,test2=3) + >>> d.test1 + 1 ''' def __init__(self, **kwargs): self.__dict__.update(kwargs) @@ -1293,7 +1264,7 @@ def betaloge(z, w): for large arguments Example - + >>> betaloge(3,2) See also @@ -1496,7 +1467,7 @@ def discretize2(fun, a,b,tol=0.005, n=5): return x, fx -def pol2cart(theta, rho): +def pol2cart(theta, rho, z=None): ''' Transform polar coordinates into 2D cartesian coordinates. @@ -1509,9 +1480,14 @@ def pol2cart(theta, rho): -------- cart2pol ''' - return rho * cos(theta), rho * sin(theta) + x, y = rho * cos(theta), rho * sin(theta) + if z is None: + return x, y + else: + return x,y,z + -def cart2pol(x, y): +def cart2pol(x, y, z=None): ''' Transform 2D cartesian coordinates into polar coordinates. Returns @@ -1525,7 +1501,12 @@ def cart2pol(x, y): -------- pol2cart ''' - return arctan2(y, x), hypot(x, y) + t, r = arctan2(y, x), hypot(x, y) + if z is None: + return t,r + else: + return t, r, z + def meshgrid(*xi, ** kwargs): """ @@ -1865,10 +1846,10 @@ def tranproc(x, f, x0, *xi): return y #y0,y1,y2,y3,y4 -def test_find_cross(): +def _test_find_cross(): t = findcross([0, 0, 1, -1, 1],0) -def test_common_shape(): +def _test_common_shape(): A = ones((4, 1)) B = 2 @@ -1883,7 +1864,7 @@ def test_common_shape(): common_shape(A, B, C, shape=(4, 5)) -def test_meshgrid(): +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) diff --git a/pywafo/src/wafo/stats/core.py b/pywafo/src/wafo/stats/core.py index be470cc..3718a65 100644 --- a/pywafo/src/wafo/stats/core.py +++ b/pywafo/src/wafo/stats/core.py @@ -206,7 +206,7 @@ def reslife(data, u=None, umin=None, umax=None, nu=None, nmin=3, alpha=0.05, plo res = WafoData(mrl, u, xlab='Threshold', ylab='Mean Excess', title=titleTxt) res.workspace = dict(numdata=num, umin=umin, umax=umax, nu=nu, nmin=nmin, alpha=alpha) res.children = [WafoData(vstack([mrll, mrlu]).T, u, xlab='Threshold', title=titleTxt)] - res.children_args = [':r'] + res.plot_args_children = [':r'] if plotflag: res.plot() return res @@ -371,7 +371,7 @@ def dispersion_idx(data, t=None, u=None, umin=None, umax=None, nu=None, nmin=10, #'caption',CItxt); res.workspace = dict(umin=umin, umax=umax, nu=nu, nmin=nmin, alpha=alpha) res.children = [WafoData(vstack([diLo * ones(nu), diUp * ones(nu)]).T, u, xlab='Threshold', title=CItxt)] - res.children_args = ['--r'] + res.plot_args_children = ['--r'] if plotflag: res.plot(di) return res, b_u, ok_u diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py index 81a000b..98f373f 100644 --- a/pywafo/src/wafo/stats/estimation.py +++ b/pywafo/src/wafo/stats/estimation.py @@ -36,6 +36,7 @@ all = alltrue def chi2isf(p, df): return special.chdtri(df, p) + def chi2sf(x, df): return special.chdtrc(df, x) diff --git a/pywafo/src/wafo/test/test_gaussian.py b/pywafo/src/wafo/test/test_gaussian.py new file mode 100644 index 0000000..d8ed75e --- /dev/null +++ b/pywafo/src/wafo/test/test_gaussian.py @@ -0,0 +1,111 @@ +''' +Created on 17. juli 2010 + +@author: pab +''' +import numpy as np +from numpy import pi, inf +from wafo.gaussian import Rind, prbnormtndpc, prbnormndpc + +def test_rind(): + ''' + >>> Et = 0.001946 # # exact prob. + >>> n = 5 + >>> Blo =-np.inf; Bup=-1.2; indI=[-1, n-1] # Barriers + >>> m = np.zeros(n); rho = 0.3; + >>> Sc =(np.ones((n,n))-np.eye(n))*rho+np.eye(n) + >>> rind = Rind() + >>> E0, err0, terr0 = rind(Sc,m,Blo,Bup,indI); + + >>> np.abs(E0-Et)< err0+terr0 + array([ True], dtype=bool) + >>> 'E0 = %2.6f' % E0 + 'E0 = 0.001946' + + >>> 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 + >>> np.abs(E1-Et)< err0+terr0 + array([ True], dtype=bool) + >>> 'E1 = %2.6f' % E1 + 'E1 = 0.001946' + + Compute expectation E( abs(X1*X2*...*X5) ) + >>> xc = np.zeros((0,1)) + >>> infinity = 37 + >>> dev = np.sqrt(np.diag(Sc)) # std + >>> ind = np.nonzero(indI[1:])[0] + >>> Bup, Blo = np.atleast_2d(Bup,Blo) + >>> Bup[0,ind] = np.minimum(Bup[0,ind] , infinity*dev[indI[ind+1]]) + >>> Blo[0,ind] = np.maximum(Blo[0,ind] ,-infinity*dev[indI[ind+1]]) + >>> rind(Sc,m,Blo,Bup,indI, xc, nt=0) + (array([ 0.05494076]), array([ 0.00083066]), array([ 1.00000000e-10])) + + Compute expectation E( X1^{+}*X2^{+} ) with random + correlation coefficient,Cov(X1,X2) = rho2. + >>> m2 = [0, 0]; + >>> rho2 = 0.3 #np.random.rand(1) + >>> Sc2 = [[1, rho2], [rho2 ,1]] + >>> Blo2 = 0; 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) + >>> E2 = g2(rho2); E2 # exact value + 0.24137214191774381 + + >>> E3, err3, terr3 = rind(Sc2,m2,Blo2,Bup2,indI2,nt=0); E3;err3;terr3 + array([ 0.24127499]) + array([ 0.00013838]) + array([ 1.00000000e-10]) + + >>> E4, err4, terr4 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0); E4;err4;terr4 + array([ 0.24127499]) + array([ 0.00013838]) + array([ 1.00000000e-10]) + + >>> E5, err5, terr5 = rind2(Sc2,m2,Blo2,Bup2,indI2,nt=0,abseps=1e-4); E5;err5;terr5 + array([ 0.24127499]) + array([ 0.00013838]) + array([ 1.00000000e-10]) + ''' +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 + >>> np.abs(E2-val2)>> 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) + >>> E3 = g3(rho3) # Exact value + >>> np.abs(E3-val3)>> 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 + >>> np.abs(E2-val2)>> 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) + >>> E3 = g3(rho3) # Exact value + >>> np.abs(E3-val3)>> np = JITImport('numpy') + >>> np.exp(0)==1.0 + True + ''' +def test_bunch(): + ''' + >>> d = Bunch(test1=1,test2=3) + >>> d.test1; d.test2 + 1 + 3 + ''' def test_dotdict(): ''' >>> d = DotDict(test1=1,test2=3) - >>> d.test1 + >>> d.test1; d.test2 1 + 3 ''' - pass def test_detrendma(): ''' @@ -123,9 +142,14 @@ def test_detrendma(): 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808])) ''' - pass + def test_findcross_and_ecross(): ''' + >>> findcross([0, 0, 1, -1, 1],0) + array([1, 2, 3]) + >>> findcross([0, 1, -1, 1],0) + array([0, 1, 2]) + >>> t = linspace(0,7*pi,250) >>> x = sin(t) >>> ind = findcross(x,0.75) @@ -136,6 +160,321 @@ def test_findcross_and_ecross(): array([ 0.84910514, 2.2933879 , 7.13205663, 8.57630119, 13.41484739, 14.85909194, 19.69776067, 21.14204343]) + ''' +def test_findextrema(): + ''' + >>> t = linspace(0,7*pi,250) + >>> x = sin(t) + >>> ind = findextrema(x) + >>> ind + array([ 18, 53, 89, 125, 160, 196, 231]) + ''' + +def test_findrfc(): + ''' + + >>> t = linspace(0,7*pi,250) + >>> x = sin(t)+0.1*sin(50*t) + >>> ind = findextrema(x) + >>> ind + array([ 1, 3, 4, 6, 7, 9, 11, 13, 14, 16, 18, 19, 21, + 23, 25, 26, 28, 29, 31, 33, 35, 36, 38, 39, 41, 43, + 45, 46, 48, 50, 51, 53, 55, 56, 58, 60, 61, 63, 65, + 67, 68, 70, 71, 73, 75, 77, 78, 80, 81, 83, 85, 87, + 88, 90, 92, 93, 95, 97, 99, 100, 102, 103, 105, 107, 109, + 110, 112, 113, 115, 117, 119, 120, 122, 124, 125, 127, 129, 131, + 132, 134, 135, 137, 139, 141, 142, 144, 145, 147, 149, 151, 152, + 154, 156, 157, 159, 161, 162, 164, 166, 167, 169, 171, 173, 174, + 176, 177, 179, 181, 183, 184, 186, 187, 189, 191, 193, 194, 196, + 198, 199, 201, 203, 205, 206, 208, 209, 211, 213, 215, 216, 218, + 219, 221, 223, 225, 226, 228, 230, 231, 233, 235, 237, 238, 240, + 241, 243, 245, 247, 248]) + >>> ti, tp = t[ind], x[ind] + >>> ind1 = findrfc(tp,0.3) + >>> ind1 + array([ 0, 9, 32, 53, 74, 95, 116, 137]) + >>> tp[ind1] + array([-0.00743352, 1.08753972, -1.07206545, 1.09550837, -1.07940458, + 1.07849396, -1.0995006 , 1.08094452]) + ''' + +def test_rfcfilter(): + ''' + # 1. Filtered signal y is the turning points of x. + >>> x = sea() + >>> y = rfcfilter(x[:,1], h=0, method=1) + >>> y[0:5] + array([-1.2004945 , 0.83950546, -0.09049454, -0.02049454, -0.09049454]) + + # 2. This removes all rainflow cycles with range less than 0.5. + >>> y1 = rfcfilter(x[:,1], h=0.5) + >>> y1[0:5] + array([-1.2004945 , 0.83950546, -0.43049454, 0.34950546, -0.51049454]) + + >>> t = linspace(0,7*pi,250) + >>> x = sin(t)+0.1*sin(50*t) + >>> ind = findextrema(x) + >>> ind + array([ 1, 3, 4, 6, 7, 9, 11, 13, 14, 16, 18, 19, 21, + 23, 25, 26, 28, 29, 31, 33, 35, 36, 38, 39, 41, 43, + 45, 46, 48, 50, 51, 53, 55, 56, 58, 60, 61, 63, 65, + 67, 68, 70, 71, 73, 75, 77, 78, 80, 81, 83, 85, 87, + 88, 90, 92, 93, 95, 97, 99, 100, 102, 103, 105, 107, 109, + 110, 112, 113, 115, 117, 119, 120, 122, 124, 125, 127, 129, 131, + 132, 134, 135, 137, 139, 141, 142, 144, 145, 147, 149, 151, 152, + 154, 156, 157, 159, 161, 162, 164, 166, 167, 169, 171, 173, 174, + 176, 177, 179, 181, 183, 184, 186, 187, 189, 191, 193, 194, 196, + 198, 199, 201, 203, 205, 206, 208, 209, 211, 213, 215, 216, 218, + 219, 221, 223, 225, 226, 228, 230, 231, 233, 235, 237, 238, 240, + 241, 243, 245, 247, 248]) + >>> ti, tp = t[ind], x[ind] + >>> tp03 = rfcfilter(tp,0.3) + >>> tp03 + array([-0.00743352, 1.08753972, -1.07206545, 1.09550837, -1.07940458, + 1.07849396, -1.0995006 , 1.08094452]) + ''' +def test_findtp(): + ''' + >>> x = sea() + >>> x1 = x[0:200,:] + >>> itp = findtp(x1[:,1],0,'Mw') + >>> itph = findtp(x1[:,1],0.3,'Mw') + >>> itp + array([ 11, 21, 22, 24, 26, 28, 31, 39, 43, 45, 47, 51, 56, + 64, 70, 78, 82, 84, 89, 94, 101, 108, 119, 131, 141, 148, + 149, 150, 159, 173, 184, 190, 199]) + >>> itph + array([ 11, 64, 28, 31, 47, 51, 39, 56, 70, 94, 78, 89, 101, + 108, 119, 148, 131, 141, 0, 159, 173, 184, 190]) + ''' +def test_findtc(): + ''' + >>> x = sea() + >>> x1 = x[0:200,:] + >>> itc, iv = findtc(x1[:,1],0,'dw') + >>> itc + array([ 28, 31, 39, 56, 64, 69, 78, 82, 83, 89, 94, 101, 108, + 119, 131, 140, 148, 159, 173, 184]) + >>> iv + array([ 19, 29, 34, 53, 60, 67, 76, 81, 82, 84, 90, 99, 103, + 112, 127, 137, 143, 154, 166, 180, 185]) + ''' + +def test_findoutliers(): + ''' + >>> xx = sea() + >>> dt = diff(xx[:2,0]) + >>> dcrit = 5*dt + >>> ddcrit = 9.81/2*dt*dt + >>> zcrit = 0 + >>> [inds, indg] = findoutliers(xx[:,1],zcrit,dcrit,ddcrit,verbose=True) + Found 0 spurious positive jumps of Dx + Found 0 spurious negative jumps of Dx + Found 37 spurious positive jumps of D^2x + Found 200 spurious negative jumps of D^2x + Found 244 consecutive equal values + Found the total of 1152 spurious points + >>> inds + array([ 6, 7, 8, ..., 9509, 9510, 9511]) + >>> indg + array([ 0, 1, 2, ..., 9521, 9522, 9523]) + ''' +def test_common_shape(): + ''' + >>> import numpy as np + >>> A = np.ones((4,1)) + >>> B = 2 + >>> C = np.ones((1,5))*5 + >>> common_shape(A,B,C) + (4, 5) + >>> common_shape(A,B,C,shape=(3,4,1)) + (3, 4, 5) + >>> A = np.ones((4,1)) + >>> B = 2 + >>> C = np.ones((1,5))*5 + >>> common_shape(A,B,C) + (4, 5) + >>> common_shape(A,B,C,shape=(3,4,1)) + (3, 4, 5) + ''' +def test_argsreduce(): + ''' + >>> import numpy as np + >>> rand = np.random.random_sample + >>> A = linspace(0,19,20).reshape((4,5)) + >>> B = 2 + >>> C = range(5) + >>> cond = np.ones(A.shape) + >>> [A1,B1,C1] = argsreduce(cond,A,B,C) + >>> B1.shape + (20,) + >>> cond[2,:] = 0 + >>> [A2,B2,C2] = argsreduce(cond,A,B,C) + >>> B2.shape + (15,) + >>> A2;B2;C2 + array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 15., + 16., 17., 18., 19.]) + array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]) + array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]) + ''' +def test_stirlerr(): + ''' + >>> stirlerr(range(5)) + array([ Inf, 0.08106147, 0.0413407 , 0.02767793, 0.02079067]) + ''' +def test_getshipchar(): + ''' + >>> sc = getshipchar(10,'service_speed') + >>> names = ['beam', 'beamSTD', 'draught', + ... 'draughtSTD', 'length', 'lengthSTD', + ... 'max_deadweight', 'max_deadweightSTD', 'propeller_diameter', + ... 'propeller_diameterSTD', 'service_speed', 'service_speedSTD'] + >>> for name in names: print( '%s : %g' % (name, sc[name])) + beam : 29 + beamSTD : 2.9 + draught : 9.6 + draughtSTD : 2.112 + length : 216 + lengthSTD : 2.01131 + max_deadweight : 30969 + max_deadweightSTD : 3096.9 + propeller_diameter : 6.76117 + propeller_diameterSTD : 0.20267 + service_speed : 10 + service_speedSTD : 0 + ''' +def test_betaloge(): + ''' + >>> betaloge(3, arange(4)) + array([ Inf, -1.09861229, -2.48490665, -3.40119738]) + ''' +def test_gravity(): + ''' + >>> phi = linspace(0,45,5) + >>> gravity(phi) + array([ 9.78049 , 9.78245014, 9.78803583, 9.79640552, 9.80629387]) + ''' +def test_nextpow2(): + ''' + >>> nextpow2(10) + 4 + >>> nextpow2(np.arange(5)) + 3 + ''' + +def test_discretize(): + ''' + >>> x, y = discretize(np.cos,0,np.pi) + >>> x; y + array([ 0. , 0.19634954, 0.39269908, 0.58904862, 0.78539816, + 0.9817477 , 1.17809725, 1.37444679, 1.57079633, 1.76714587, + 1.96349541, 2.15984495, 2.35619449, 2.55254403, 2.74889357, + 2.94524311, 3.14159265]) + array([ 1.00000000e+00, 9.80785280e-01, 9.23879533e-01, + 8.31469612e-01, 7.07106781e-01, 5.55570233e-01, + 3.82683432e-01, 1.95090322e-01, 6.12323400e-17, + -1.95090322e-01, -3.82683432e-01, -5.55570233e-01, + -7.07106781e-01, -8.31469612e-01, -9.23879533e-01, + -9.80785280e-01, -1.00000000e+00]) + ''' +def test_discretize2(): + ''' + >>> x, y = discretize2(np.cos,0,np.pi) + >>> x; y + array([ 0. , 0.19634954, 0.39269908, 0.58904862, 0.78539816, + 0.9817477 , 1.17809725, 1.37444679, 1.57079633, 1.76714587, + 1.96349541, 2.15984495, 2.35619449, 2.55254403, 2.74889357, + 2.94524311, 3.14159265]) + array([ 1.00000000e+00, 9.80785280e-01, 9.23879533e-01, + 8.31469612e-01, 7.07106781e-01, 5.55570233e-01, + 3.82683432e-01, 1.95090322e-01, 6.12323400e-17, + -1.95090322e-01, -3.82683432e-01, -5.55570233e-01, + -7.07106781e-01, -8.31469612e-01, -9.23879533e-01, + -9.80785280e-01, -1.00000000e+00]) + ''' +def test_pol2cart_n_cart2pol(): + ''' + >>> r = 5 + >>> t = linspace(0,pi,20) + >>> x, y = pol2cart(t,r) + >>> x; y + array([ 5. , 4.93180652, 4.72908621, 4.39736876, 3.94570255, + 3.38640786, 2.73474079, 2.00847712, 1.22742744, 0.41289673, + -0.41289673, -1.22742744, -2.00847712, -2.73474079, -3.38640786, + -3.94570255, -4.39736876, -4.72908621, -4.93180652, -5. ]) + array([ 0.00000000e+00, 8.22972951e-01, 1.62349735e+00, + 2.37973697e+00, 3.07106356e+00, 3.67861955e+00, + 4.18583239e+00, 4.57886663e+00, 4.84700133e+00, + 4.98292247e+00, 4.98292247e+00, 4.84700133e+00, + 4.57886663e+00, 4.18583239e+00, 3.67861955e+00, + 3.07106356e+00, 2.37973697e+00, 1.62349735e+00, + 8.22972951e-01, 6.12323400e-16]) + >>> ti, ri = cart2pol(x,y) + >>> ti;ri + array([ 0. , 0.16534698, 0.33069396, 0.49604095, 0.66138793, + 0.82673491, 0.99208189, 1.15742887, 1.32277585, 1.48812284, + 1.65346982, 1.8188168 , 1.98416378, 2.14951076, 2.31485774, + 2.48020473, 2.64555171, 2.81089869, 2.97624567, 3.14159265]) + array([ 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., + 5., 5., 5., 5., 5., 5., 5.]) + ''' +def test_meshgrid(): + ''' + >>> x = np.linspace(0,1,3) # coordinates along x axis + >>> y = np.linspace(0,1,2) # coordinates along y axis + >>> xv, yv = meshgrid(x,y) # extend x and y for a 2D xy grid + >>> xv + array([[ 0. , 0.5, 1. ], + [ 0. , 0.5, 1. ]]) + >>> yv + array([[ 0., 0., 0.], + [ 1., 1., 1.]]) + >>> xv, yv = meshgrid(x,y, sparse=True) # make sparse output arrays + >>> xv + array([[ 0. , 0.5, 1. ]]) + >>> yv + array([[ 0.], + [ 1.]]) + + >>> meshgrid(x,y,sparse=True,indexing='ij') # change to matrix indexing + [array([[ 0. ], + [ 0.5], + [ 1. ]]), array([[ 0., 1.]])] + >>> meshgrid(x,y,indexing='ij') + [array([[ 0. , 0. ], + [ 0.5, 0.5], + [ 1. , 1. ]]), array([[ 0., 1.], + [ 0., 1.], + [ 0., 1.]])] + + >>> meshgrid(0,1,5) # just a 3D point + [array([[[0]]]), array([[[1]]]), array([[[5]]])] + >>> map(np.squeeze,meshgrid(0,1,5)) # just a 3D point + [array(0), array(1), array(5)] + >>> meshgrid(3) + array([3]) + >>> meshgrid(y) # 1D grid y is just returned + array([ 0., 1.]) + + `meshgrid` is very useful to evaluate functions on a grid. + + >>> x = np.arange(-5, 5, 0.1) + >>> y = np.arange(-5, 5, 0.1) + >>> xx, yy = meshgrid(x, y, sparse=True) + >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2) + ''' +def test_tranproc(): + ''' + >>> import wafo.transform.models as wtm + >>> tr = wtm.TrHermite() + >>> x = linspace(-5,5,501) + >>> g = tr(x) + >>> y0, y1 = tranproc(x, g, range(5), ones(5)) + >>> y0;y1 + array([ 0.02659612, 1.00115284, 1.92872532, 2.81453257, 3.66292878]) + array([ 1.00005295, 0.9501118 , 0.90589954, 0.86643821, 0.83096482]) + ''' if __name__=='__main__': import doctest