updated .travis.yml

moved some funtions from numpy_utils -> misc.
pep8

.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:

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__':

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()
 
 

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
@@ -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)
@@ -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)
     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,14 +1331,20 @@ 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::]):
@@ -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,8 +1756,7 @@ 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.' %
+                print('Found %d consecutive values less than %g apart.' %
                       (indz.size, zcrit))
     indg = ones(xn.size, dtype=bool)
 
@@ -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,8 +1918,7 @@ def stirlerr(n):
     -----------
     Catherine Loader (2000).
     Fast and Accurate Computation of Binomial Probabilities
-    <http://www.herine.net/stat/software/dbinom.html>
-    <http://www.citeseer.ist.psu.edu/312695.html>
+    <http://lists.gnu.org/archive/html/octave-maintainers/2011-09/pdfK0uKOST642.pdf>
     '''
 
     S0 = 0.083333333333333333333  # /* 1/12 */
@@ -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))
@@ -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 "<pre>"
+    print("<pre>")
     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 "</pre>"
+    print("</pre>")
 
 
 main = profile_main1

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,

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.
 
@@ -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,20 +1990,21 @@ 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')
+        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
@@ -2015,11 +2017,13 @@ def chebvandernd(deg, *xi):
 
     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
@@ -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
     -------
@@ -2200,6 +2207,7 @@ def chebvalnd(c, *xi):
         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`.
@@ -2252,14 +2260,13 @@ def chebgridnd(c, *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]
 
@@ -2276,21 +2283,23 @@ 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)))
+    # 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])
+    _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)
     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 = chebvalnd(dcoeff, xgrid, ygrid)

wafo/powerpoint.py

@@ -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')
@@ -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:

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)

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,16 +109,20 @@ 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],
+    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

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()

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):