Added more test to test/test_misc.py and test/test_gaussian.py

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)<err2
+    True
+    
+    >>> 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)<err2
+    True
+    
+      
+    See also
+    --------  
+    prbnormndpc, prbnormnd, rind
+      
+    Reference
+    --------- 
+    Charles Dunnett (1989)
+    "Multivariate normal probability integrals with product correlation
+    structure", Applied statistics, Vol 38,No3, (Algorithm AS 251)    
+    '''
+      
+    if D is None:
+        D = zeros(len(rho))
+    # Make sure integration limits are finite
+    A = np.clip(a-D,-100,100)
+    B = np.clip(b-D,-100,100)
+    return mvnprdmod.prbnormtndpc(rho,A,B,df,abseps,IERC,HNC)
+
+def prbnormndpc(rho,a,b,abserr=1e-4,relerr=1e-4,usesimpson=True, usebreakpoints=False):
+    '''
+    Return Multivariate Normal probabilities with product correlation
+     
+    Parameters
+    ----------
+      rho  = vector defining the correlation structure, i.e., 
+              corr(Xi,Xj) = rho(i)*rho(j) for i~=j
+                          = 1             for i==j  
+                 -1 <= rho <= 1  
+      a,b   = lower and upper integration limits respectively.  
+      tol   = requested absolute tolerance
+      
+    Returns
+    -------
+    value = value of integral
+    error = estimated absolute error
+    
+    PRBNORMNDPC calculates multivariate normal probability
+    with product correlation structure for rectangular regions.
+    The accuracy is up to almost double precision, i.e., about 1e-14.
+      
+    Example:
+    -------
+    >>> 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)<err2
+    True
+    
+    >>> 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)<err2
+    True
+    
+    See also
+    --------
+    prbnormtndpc, prbnormnd, rind
+      
+    Reference
+    ---------
+    P. A. Brodtkorb (2004), 
+    "Evaluating multinormal probabilites with product correlation structure."
+    In Lund university report series
+    and in the Dr.Ing thesis: 
+    "The probability of Occurrence of dangerous Wave Situations at Sea."
+    Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
+    Trondheim, Norway.
+      
+    '''  
+    # Call fortran implementation
+    val,err,ier = mvnprdmod.prbnormndpc(rho,a,b,abserr,relerr,usebreakpoints,usesimpson);
+    
+    if ier>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<max(50*rel.mach.acc.,0.5d-28))
+        4) limit < npts2.'''
+    
+
+def  prbnormnd(correl,a,b,abseps=1e-4,releps=1e-3,maxpts=None,method=0):
+    '''
+    
+    Multivariate Normal probability by Genz' algorithm.
+    
+    
+    Parameters
+    CORREL = Positive semidefinite correlation matrix
+    A         = vector of lower integration limits.
+    B         = vector of upper integration limits.
+    ABSEPS = absolute error tolerance.
+    RELEPS = relative error tolerance.
+    MAXPTS = maximum number of function values allowed. This 
+                  parameter can be used to limit the time. A sensible 
+                  strategy is to start with MAXPTS = 1000*N, and then
+                  increase MAXPTS if ERROR is too large.
+    METHOD = integer defining the integration method
+                 -1 KRBVRC randomized Korobov rules for the first 20
+                    variables, randomized Richtmeyer rules for the rest, 
+                    NMAX = 500 
+                  0 KRBVRC, NMAX = 100 (default)
+                  1 SADAPT Subregion Adaptive integration method, NMAX = 20 
+                  2 KROBOV Randomized KOROBOV rules,              NMAX = 100
+                  3 RCRUDE Crude Monte-Carlo Algorithm with simple
+                    antithetic variates and weighted results on restart 
+                  4 SPHMVN Monte-Carlo algorithm by Deak (1980),  NMAX = 100
+    Returns
+    -------
+    VALUE  REAL estimated value for the integral
+    ERROR  REAL estimated absolute error, with 99% confidence level.
+    INFORM INTEGER, termination status parameter:
+                if INFORM = 0, normal completion with ERROR < EPS;
+                if INFORM = 1, completion with ERROR > 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: #  

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 (abs, amax, any, logical_and, arange, linspace, atleast_1d, 
-from numpy import amax
+                   array, asarray, broadcast_arrays, ceil,  floor, frexp, hypot,
-from numpy import any
+                   sqrt, arctan2, sin, cos, exp, log, mod, diff, empty_like, 
-from numpy import arange
+                   finfo, inf, pi, interp, isnan, isscalar, zeros, ones, 
-from numpy import arctan2
+                   r_, sign, unique, hstack, vstack, nonzero, where, extract)
-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 scipy.special import gammaln
 import types
 import warnings
@@ -100,6 +66,11 @@ 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)

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

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)
 

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)<err2
+    True
+    
+    >>> 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)<err2
+    True
+    '''
+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
+    >>> np.abs(E2-val2)<err2
+    True
+    
+    >>> 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)<err2
+    True
+    '''
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()

pywafo/src/wafo/test/test_misc.py

@@ -1,14 +1,33 @@
-from pylab import cos, randn, exp, linspace, pi, sin
+import numpy as np
-from wafo.misc import detrendma, DotDict, findcross, ecross
+from numpy import cos, exp, linspace, pi, sin, diff, arange, ones
-
+from numpy.random import randn
+from wafo.data import sea
+from wafo.misc import (JITImport, Bunch, detrendma, DotDict, findcross, ecross, findextrema, 
+                       findrfc, rfcfilter, findtp, findtc, findoutliers, 
+                       common_shape, argsreduce, stirlerr, getshipchar, betaloge, 
+                       gravity, nextpow2, discretize, discretize2, pol2cart, 
+                       cart2pol, meshgrid, tranproc)
 
+def test_JITImport():
+    '''
+    >>> 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