Small updates

pywafo/src/Wafo.egg-info/SOURCES.txt

@@ -12,6 +12,7 @@ src/Wafo.egg-info/dependency_links.txt
 src/Wafo.egg-info/top_level.txt
 src/wafo/SpecData1D.mm
 src/wafo/__init__.py
+src/wafo/bitwise.py
 src/wafo/c_library.pyd
 src/wafo/dctpack.py
 src/wafo/definitions.py
@@ -56,6 +57,14 @@ src/wafo/data/wafoLogoNewWithoutBorder.png
 src/wafo/data/wafoLogoNewWithoutBorder.svg
 src/wafo/data/wafologoWithBorder.png
 src/wafo/data/yura87.dat
+src/wafo/doc/__init__.py
+src/wafo/doc/test_pyreport.py
+src/wafo/doc/tutorial_scripts/__init__.py
+src/wafo/doc/tutorial_scripts/chapter1.py
+src/wafo/doc/tutorial_scripts/chapter2.py
+src/wafo/doc/tutorial_scripts/chapter3.py
+src/wafo/doc/tutorial_scripts/chapter4.py
+src/wafo/doc/tutorial_scripts/chapter5.py
 src/wafo/source/c_codes/build_all.py
 src/wafo/source/c_codes/c_functions.c
 src/wafo/source/c_codes/c_library.pyd
@@ -264,8 +273,11 @@ src/wafo/stats/core.py
 src/wafo/stats/distributions.py
 src/wafo/stats/estimation.py
 src/wafo/stats/misc.py
+src/wafo/stats/tests/test_distributions.py
+src/wafo/stats/tests/test_estimation.py
 src/wafo/test/__init__.py
 src/wafo/test/test_gaussian.py
+src/wafo/test/test_kdetools.py
 src/wafo/test/test_misc.py
 src/wafo/transform/__init__.py
 src/wafo/transform/core.py

pywafo/src/wafo/covariance/core.py

@@ -121,7 +121,7 @@ class CovData1D(WafoData):
     --------
     >>> import numpy as np
     >>> import wafo.spectrum as sp
-    >>> Sj = sp.models.Jonswap(Hm0=3)
+    >>> Sj = sp.models.Jonswap(Hm0=3,Tp=7)
     >>> w = np.linspace(0,4,256)
     >>> S = sp.SpecData1D(Sj(w),w) #Make spectrum object from numerical values
 

pywafo/src/wafo/integrate.py

@@ -2,14 +2,40 @@ from __future__ import division
 import warnings
 import copy
 import numpy as np
-from numpy import pi, sqrt, ones, zeros #@UnresolvedImport
+from numpy import pi, sqrt, ones, zeros, exp #@UnresolvedImport
 from scipy import integrate as intg
 import scipy.special.orthogonal as ort
 from scipy import special as sp
 import pylab as plb
+from wafo.misc import meshgrid
 
 _POINTS_AND_WEIGHTS = {}
 
+def peaks(x=None, y=None, n=51):
+    '''
+    Return the "well" known MatLab (R) peaks function
+    evaluated in the [-3,3] x,y range
+    
+    Example
+    -------
+    >>> import pylab as plt
+    >>> x,y,z = peaks()
+    >>> plt.contourf(x,y,z)
+    
+    '''
+    if x is None:
+        x = np.linspace(-3, 3, n)
+    if y is None:
+        y = np.linspace(-3, 3, n)
+        
+    [x1, y1] = meshgrid(x, y)
+
+    z = (3 * (1 - x1) ** 2 * exp(-(x1 ** 2) - (y1 + 1) ** 2) 
+         - 10 * (x1 / 5 - x1 ** 3 - y1 ** 5) * exp(-x1 ** 2 - y1 ** 2) 
+         - 1. / 3 * exp(-(x1 + 1) ** 2 - y1 ** 2))
+
+    return x1, y1, z
+
 def humps(x=None):
     '''
     Computes a function that has three roots, and some humps.

pywafo/src/wafo/misc.py

@@ -4,7 +4,7 @@ Misc lsdkfalsdflasdfl
 from __future__ import division
 
 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,
@@ -384,8 +384,89 @@ def findextrema(x):
     '''
     xn = atleast_1d(x).ravel()
     return findcross(diff(xn), 0.0) + 1
+def findpeaks(data, n=2, min_h=None, min_p=0.0):
+    '''
+    Find peaks of vector or matrix possibly rainflow filtered
+     
+    Parameters
+    ----------
+    data      = matrix or vector
+    n     = The n highest peaks are found (if exist). (default 2)
+    min_h  = The threshold in the rainflowfilter (default 0.05*range(S(:))).
+                 A zero value will return all the peaks of S.
+    min_p  = 0..1, Only the peaks that are higher than 
+                       min_p*max(max(S))  min_p*(the largest peak in S)
+                       are returned (default  0).
+    Returns
+    ix = 
+        linear index to peaks of S
+        
+    Example: 
+    
+    Find highest 8 peaks that are not
+    less that 0.3*"global max" and have 
+    rainflow amplitude larger than 5.
+    >>> import numpy as np
+    >>> x = np.arange(0,10,0.01) 
+    >>> data = x**2+10*np.sin(3*x)+0.5*np.sin(50*x)
+    >>> findpeaks(data, n=8, min_h=5, min_p=0.3) 
+    array([908, 694, 481])
     
-def findrfc(tp, hmin=0.0):
+    See also  
+    --------
+    findtp
+    '''
+    S = np.atleast_1d(data)
+    smax = S.max()
+    if min_h is None:
+        smin = S.min()
+        min_h = 0.05*(smax-smin)
+    ndim = S.ndim
+    S = np.atleast_2d(S)
+    nrows, mcols = S.shape
+
+    # Finding turningpoints of the spectrum
+    # Returning only those with rainflowcycle heights greater than h_min 
+    indP = [] # indices to peaks
+    ind = []
+    for iy in range(nrows): # % find all peaks
+        TuP = findtp(S[iy], min_h)
+        if len(TuP):
+            ind = TuP[1::2] #; % extract indices to maxima only
+        else: # % did not find any , try maximum
+            ind = S[iy].argmax()
+      
+        if ndim>1:
+            if iy==0:
+                ind2 = np.flatnonzero(S[iy,ind]>S[iy+1,ind])
+            elif iy==nrows-1:
+                ind2 = np.flatnonzero(S[iy,ind]>S[iy-1,ind])
+            else:
+                ind2 = np.flatnonzero((S[iy,ind]>S[iy-1,ind]) & (S[iy,ind]>S[iy+1,ind]))
+        
+            if len(ind2):
+                indP.append((ind[ind2] + iy*mcols))
+      
+    if ndim>1:
+        ind = np.hstack(indP) if len(indP) else []
+    if len(ind)==0:
+        return []
+    
+    peaks = S.take(ind)
+    ind2 = peaks.argsort()[::-1]
+     
+    
+    # keeping only the Np most significant peak frequencies.
+    nmax = min(n,len(ind))
+    ind = ind[ind2[:nmax]]
+    if (min_p >0 ) :
+        # Keeping only peaks larger than min_p percent relative to the maximum peak 
+        ind = ind[(S.take(ind) > min_p*smax)]
+    
+    return ind
+    
+
+def findrfc(tp, hmin=0.0, method='clib'):
     '''
     Return indices to rainflow cycles of a sequence of TP.
 
@@ -410,7 +491,8 @@ def findrfc(tp, hmin=0.0):
     >>> ind = findextrema(x)
     >>> ti, tp = t[ind], x[ind]
     >>> a = pb.plot(t,x,'.',ti,tp,'r.')
-    >>> ind1 = findrfc(tp,0.3)
+    >>> ind1 = findrfc(tp,0.3); ind1
+    array([  0,   9,  32,  53,  74,  95, 116, 137])
     >>> a = pb.plot(ti[ind1],tp[ind1])
     >>> pb.close('all')
 
@@ -445,7 +527,7 @@ def findrfc(tp, hmin=0.0):
         warnings.warn('This is not a sequence of turningpoints, exit')
         return ind
 
-    if clib is None:
+    if clib is None or method!='clib':
         ind = zeros(n, dtype=np.int)
         NC = np.int(NC)
         for i in xrange(NC):
@@ -506,7 +588,7 @@ def findrfc(tp, hmin=0.0):
         #  /* for i */
     else:
         ind, ix = clib.findrfc(y, hmin)
-    return ind[:ix]
+    return np.sort(ind[:ix])
 
 def rfcfilter(x, h, method=0):
     """ 
@@ -653,8 +735,8 @@ def findtp(x, h=0.0, kind=None):
             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,
+    array([ 11,  28,  31,  39,  47,  51,  56,  64,  70,  78,  89,  94, 101,
-           108, 119, 148, 131, 141,   0, 159, 173, 184, 190])
+           108, 119, 131, 141, 148, 159, 173, 184, 190, 199])
 
     See also
     ---------
@@ -1911,6 +1993,40 @@ def histgrm(data, n=None, odd=False, scale=False, lintype='b-'):
     binwidth = d
     return binwidth
 
+def num2pistr(x, n=3):
+    ''' 
+    Convert a scalar to a text string in fractions of pi
+        if the numerator is less than 10 and not equal 0 
+               and if the denominator is less than 10.
+
+    Parameters
+    ----------
+    x   = a scalar
+    n   = maximum digits of precision. (default 3)
+    Returns
+    -------
+    xtxt = a text string in fractions of pi
+    
+    Example
+    >>> num2pistr(np.pi*3/4)
+    '3\\pi/4'
+    '''
+    
+    frac = fractions.Fraction.from_float(x/pi).limit_denominator(10000000)
+    num = frac.numerator
+    den = frac.denominator
+    if (den<10) and (num<10) and (num!=0):
+        dtxt = '' if abs(den)==1 else '/%d' % den
+        if abs(num)==1: # % numerator
+            ntxt='-' if num==-1 else ''
+        else: 
+            ntxt = '%d' % num 
+        xtxt= ntxt+r'\pi'+dtxt
+    else:
+        format = '%0.' +'%dg' % n
+        xtxt = format % x
+    return xtxt     
+
 def _test_find_cross():
     t = findcross([0, 0, 1, -1, 1], 0)
     

pywafo/src/wafo/spectrum/__init__.py

@@ -5,6 +5,6 @@ Spectrum package in WAFO Toolbox.
 
 """
 
-from core import SpecData1D
+from core import SpecData1D, SpecData2D, cltext
 import models
 import dispersion_relation

pywafo/src/wafo/spectrum/dispersion_relation.py

@@ -136,7 +136,7 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100):
         return zeros_like(wi)
 
     k = 1.0*sign(wi)*wi**2.0  #% deep water
-    if h == inf:
+    if h > 10. ** 25:
         k2 = k*sin(th)/gi[-1] #%size np x nf
         k1 = k*cos(th)/gi[0]
         return k1, k2

pywafo/src/wafo/spectrum/models.py

@@ -47,13 +47,13 @@ import scipy.integrate as integrate
 import scipy.special as sp
 from scipy.fftpack import fft
 #from scipy.misc import ppimport
-#import numpy as np
+import numpy as np
 from numpy import (inf, atleast_1d, newaxis, any, minimum, maximum, array, #@UnresolvedImport
     asarray, exp, log, sqrt, where, pi, arange, linspace, sin, cos, abs, sinh, #@UnresolvedImport
     isfinite, mod, expm1, tanh, cosh, finfo, ones, ones_like, isnan, #@UnresolvedImport
     zeros_like, flatnonzero, sinc, hstack, vstack, real, flipud, clip) #@UnresolvedImport
 from dispersion_relation import w2k
-from wafo.spectrum import SpecData1D
+from wafo.spectrum import SpecData1D, SpecData2D
 sech = lambda x: 1.0 / cosh(x)
 
 eps = finfo(float).eps
@@ -1492,9 +1492,9 @@ class Spreading(object):
                                     p=self.fourier2x, s=self.fourier2b,
                                     w=self.fourier2d)
 
-    def __call__(self, *args):
+    def __call__(self, *args, **kwds):
         spreadfun = self._spreadfun[self.type[0]]
-        return spreadfun(*args)
+        return spreadfun(*args, **kwds)
 
     def  chk_input(self, theta, w=1, wc=1): # [s_par,TH,phi0,Nt] =
         ''' CHK_INPUT
@@ -1828,10 +1828,6 @@ class Spreading(object):
         r = clip(r1, 0., 1.0)
         return where(r <= 0, inf, sqrt(-2.0 * log(r)))
 
-
-
-
-
     def  spread_par_s(self, wn):
         ''' Return spread parameter, S, of COS2S function
 
@@ -1906,7 +1902,66 @@ class Spreading(object):
         xk = pi / x
         return where(x < 100., xk / sinh(xk), -2. * xk / (exp(xk) * expm1(-2. * xk)))
 
+    def tospecdata2d(self, specdata=None, theta=None, wc=0.52, nt=51):
+        '''
+         MKDSPEC Make a directional spectrum
+                 frequency spectrum times spreading function
+        
+         CALL:  Snew=mkdspec(S,D,plotflag)
+          
+               Snew = directional spectrum (spectrum struct)
+               S    = frequency spectrum (spectrum struct)
+                          (default jonswap)  
+               D    = spreading function (special struct)
+                          (default spreading([],'cos2s'))
+               plotflag = 1: plot the spectrum, else: do not plot (default 0)   
+        
+         Creates a directional spectrum through multiplication of a frequency
+         spectrum and a spreading function: S(w,theta)=S(w)*D(w,theta)
           
+         The spreading structure must contain the following fields:
+           .S (size [np 1] or [np nf])  and  .theta (length np)  
+         optional fields: .w (length nf), .note (memo) .phi (rotation-azymuth)   
+          
+         NB! S.w and D.w (if any) must be identical.
+        
+         Example 
+         -------
+         >>> S = Jonswap().tospecdata()
+         >>> D = Spreading('cos2s')
+         >>> SD = D.tospecdata2d(S)
+         >>> SD.plot() 
+          
+         See also  spreading, rotspec, jonswap, torsethaugen  
+        '''  
+        
+        if specdata is None:
+            specdata = Jonswap().tospecdata()
+        if theta is None:
+            pi = np.pi
+            theta = np.linspace(-pi,pi,nt)
+        else:
+            L = abs(theta[-1]-theta[0])
+            if abs(L-pi)>eps: 
+                raise ValueError('theta must cover all angles -pi -> pi')   
+            nt = len(theta)
+
+        if nt<40: 
+            warnings.warn('Number of angles is less than 40. Spreading too sparsely sampled!')
+
+        w = specdata.args
+        S = specdata.data
+        D, phi0 = self(theta, w=w, wc=wc)
+            
+        SD = D * S[None,:]
+        
+        Snew    = SpecData2D(SD,(w,theta), type='dir', freqtype=specdata.freqtype)
+        Snew.tr = specdata.tr
+        Snew.h  = specdata.h
+        Snew.phi = phi0
+        Snew.norm = specdata.norm
+        #Snew.note = specdata.note + ', spreading: %s' % self.type
+        return Snew
 
 def test_some_spectra():
     S = Jonswap()

pywafo/src/wafo/stats/distributions.py

@@ -2075,7 +2075,7 @@ class rv_continuous(rv_generic):
                 i_tie = nonzero(tie)
                 tiedata = x[i_tie]
 
-                logD[(i_tie[0] + 1,)] = log(self._pdf(tiedata, *args)) + log(scale)
+                logD[i_tie + 1] = log(self._pdf(tiedata, *args)) + log(scale)
 
             finiteD = numpy.isfinite(logD)
             nonfiniteD = 1 - finiteD
@@ -3754,6 +3754,23 @@ class genextreme_gen(rv_continuous):
         k = arange(0, n + 1)
         vals = 1.0 / c ** n * sum(comb(n, k) * (-1) ** k * special.gamma(c * k + 1), axis=0)
         return where(c * n > -1, vals, inf)
+    def _fitstart(self, data):
+        d = arr(data)
+        #Probability weighted moments
+        log = np.log
+        n = len(d)
+        d.sort()
+        koeff1 = np.r_[0:n] / (n - 1)
+        koeff2 = koeff1 * (np.r_[0:n] - 1) / (n - 2)
+        b2 = np.dot(koeff2, d) / n
+        b1 = np.dot(koeff1, d) / n
+        b0 = d.mean()
+        z = (2 * b1 - b0) / (3 * b2 - b0) - log(2) / log(3)
+        shape = 7.8590 * z + 2.9554 * z ** 2
+        scale = (2 * b1 - b0) * shape / (exp(gamln(1 + shape)) * (1 - 2 ** (-shape)))
+        loc = b0 + scale * (expm1(gamln(1 + shape))) / shape
+        return shape, loc, scale 
+        
 genextreme = genextreme_gen(name='genextreme',
                             longname="A generalized extreme value",
                             shapes='c', extradoc="""