Fixed more bugs in distributions.py

pywafo/src/wafo/misc.py

@@ -7,12 +7,13 @@ import sys
 import fractions
 import numpy as np
 from numpy import (
-    abs, amax, any, logical_and, arange, linspace, atleast_1d,  # atleast_2d,
+                   meshgrid,
-    array, asarray, broadcast_arrays, ceil, floor, frexp, hypot,
+    abs, amax, any, logical_and, arange, linspace, atleast_1d,
+    array, asarray, ceil, floor, frexp, hypot,
     sqrt, arctan2, sin, cos, exp, log, mod, diff, empty_like,
     finfo, inf, pi, interp, isnan, isscalar, zeros, ones, linalg,
     r_, sign, unique, hstack, vstack, nonzero, where, extract)
-from scipy.special import gammaln
+from scipy.special import gammaln, gamma, psi
 from scipy.integrate import trapz, simps
 import warnings
 from plotbackend import plotbackend
@@ -24,7 +25,8 @@ try:
 except:
     clib = None
 floatinfo = finfo(float)
-
+_TINY = np.finfo(float).tiny
+_EPS = np.finfo(float).eps
 
 __all__ = [
     'is_numlike', 'JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_select',
@@ -1693,6 +1695,600 @@ def gravity(phi=45):
                       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 + 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 hyp2f1_taylor(a, b, c, z, tol=1e-13, itermax=500):
+    a, b, c, z = np.broadcast_arrays(*np.atleast_1d(a, b, c, z))
+    shape = a.shape
+    ak, bk, ck, zk = [d.ravel() for d in (a, b, c, z)]
+    ajm1 = np.ones(ak.shape)
+    bjm2 = 0.5 * np.ones(ak.shape)
+    bjm1 = np.ones(ak.shape)
+    hout = np.zeros(ak.shape)
+    k0 = np.arange(len(ak))
+    for j in range(0, itermax):
+        aj = ajm1 * (ak + j) * (bk + j) / (ck + j) * zk / (j + 1)
+        bj = bjm1 + aj
+        h, err = dea3(bjm2, bjm1, bj)
+        k = np.flatnonzero(err > tol * np.abs(h))
+        hout[k0] = h
+        if len(k) == 0:
+            break
+        k0 = k0[k]
+        ak, bk, ck, zk = ak[k], bk[k], ck[k], zk[k]
+        ajm1 = aj[k]
+        bjm2 = bjm1[k]
+        bjm1 = bj[k]
+    else:
+        warnings.warn(('Reached %d limit! \n' +
+                      '#%d values did not converge! Max error=%g') %
+                      (j, len(k), np.max(err)))
+    return hout.reshape(shape)
+
+
+def hyp2f1(a, b, c, z, rho=0.5):
+    e1 = gammaln(a)
+    e2 = gammaln(b)
+    e3 = gammaln(c)
+    e4 = gammaln(b - a)
+    e5 = gammaln(a - b)
+
+    e6 = gammaln(c - a)
+    e7 = gammaln(c - b)
+    e8 = gammaln(c - a - b)
+    e9 = gammaln(a + b - c)
+    _cmab = c-a-b
+    #~(np.round(cmab) == cmab & cmab <= 0)
+    if abs(z) <= rho:
+        h = hyp2f1_taylor(a, b, c, z, 1e-15)
+    elif abs(1 - z) <= rho:  # % Require that |arg(1-z)|<pi
+        h = exp(e3 + e8 - e6 - e7) * hyp2f1_taylor(a, b, a + b - c, 1 - z, 1e-15) \
+            + (1 - z) ** (c - a - b) * exp(e3 + e9 - e1 - e2) \
+            * hyp2f1_taylor(c - a, c - b, c - a - b + 1, 1 - z, 1e-15)
+    elif abs(z / (z - 1)) <= rho:
+        h = (1 - z) ** (-a) \
+            * hyp2f1_taylor(a, c - b, c, (z / (z - 1)), 1e-15)
+    elif abs(1 / z) <= rho:  # % Require that |arg(z)|<pi and |arg(1-z)|<pi
+        h = (-z + 0j) ** (-a) * exp(e3 + e4 - e2 - e6) \
+            * hyp2f1_taylor(a, a - c + 1, a - b + 1, 1. / z, 1e-15) \
+            + (-z + 0j) ** (-b) * exp(e3 + e5 - e1 - e7) \
+            * hyp2f1_taylor(b - c + 1, b, b - a + 1, (1. / z), 1e-15)
+    elif abs(1. / (1 - z)) <= rho:  # % Require that |arg(1-z)|<pi
+        h = (1 - z) ** (-a) * exp(e3 + e4 - e2 - e6) \
+            * hyp2f1_taylor(a, c - b, a - b + 1, (1. / (1 - z)), 1e-15)\
+            + (1 - z) ** (-b) * exp(e3 + e5 - e1 - e7) \
+            * hyp2f1_taylor(b, c - a, b - a + 1, (1. / (1 - z)), 1e-15)
+    elif abs(1 - 1 / z) < rho:  # % Require that |arg(z)|<pi and |arg(1-z)|<pi
+        h = z ** (-a) * exp(e3 + e8 - e6 - e7) \
+            * hyp2f1_taylor(a, a - c + 1, a + b - c + 1, (1 - 1 / z), 1e-15) \
+            + z ** (a - c) * (1 - z) ** (c - a - b) * exp(e3 + e9 - e1 - e2) \
+            * hyp2f1_taylor(c - a, 1 - a, c - a - b + 1, (1 - 1 / z), 1e-15)
+    else:
+        warnings.warn('Another method is needed')
+    return h
+
+
+def hyp2f1_wrong(a, b, c, z, tol=1e-13, itermax=500):
+    ajm1 = 0
+    bjm1 = 1
+    cjm1 = 1
+    xjm1 = np.ones(np.shape(c + a * b * z))
+    xjm2 = 2 * np.ones(xjm1.shape)
+
+    for j in range(1, itermax):
+        aj = (ajm1 + bjm1) * j * (c + j - 1)
+        bj = bjm1 * (a + j - 1) * (b + j - 1) * z
+        cj = cjm1 * j * (c + j - 1)
+        if np.any((aj == np.inf) | (bj == np.inf) | (cj == np.inf)):
+            break
+        xj = (aj + bj) / cj
+        h, err = dea3(xjm2, xjm1, xj)
+        if np.all(err <= tol * np.abs(h)) and j > 10:
+            break
+        xjm2 = xjm1
+        xjm1 = xj
+    else:
+        warnings.warn('Reached %d limit' % j)
+    return h
+
+
+def hygfz(A, B, C, Z):
+    ''' Return hypergeometric function for a complex argument, F(a,b,c,z)
+
+    Parameters
+    ----------
+     a, b, c:
+         parameters where c <> 0,-1,-2,...
+    z :--- Complex argument
+    '''
+    X = np.real(Z)
+    Y = np.imag(Z)
+    EPS = 1.0e-15
+    L0 = C == np.round(C) and C < 0.0e0
+    L1 = abs(1.0 - X) < EPS and Y == 0.0 and C - A - B <= 0.0
+    L2 = abs(Z + 1.0) < EPS and abs(C - A + B - 1.0) < EPS
+    L3 = A == np.round(A) and A < 0.0
+    L4 = B == np.round(B) and B < 0.0
+    L5 = C - A == np.round(C - A) and C - A <= 0.0
+    L6 = C - B == np.round(C - B) and C - B <= 0.0
+    AA = A
+    BB = B
+    A0 = abs(Z)
+    if (A0 > 0.95):
+        EPS = 1.0e-8
+    PI = 3.141592653589793
+    EL = .5772156649015329
+    if (L0 or L1):
+    # 'The hypergeometric series is divergent'
+        return np.inf
+
+    NM = 0
+    if (A0 == 0.0 or A == 0.0 or B == 0.0):
+        ZHF = 1.0
+    elif (Z == 1.0 and C - A - B > 0.0):
+        GC = gamma(C)
+        GCAB = gamma(C - A - B)
+        GCA = gamma(C - A)
+        GCB = gamma(C - B)
+        ZHF = GC * GCAB / (GCA * GCB)
+    elif L2:
+        G0 = sqrt(PI) * 2.0 ** (-A)
+        G1 = gamma(C)
+        G2 = gamma(1.0 + A / 2.0 - B)
+        G3 = gamma(0.5 + 0.5 * A)
+        ZHF = G0 * G1 / (G2 * G3)
+    elif L3 or L4:
+        if (L3):
+            NM = int(np.round(abs(A)))
+        if (L4):
+            NM = int(np.round(abs(B)))
+        ZHF = 1.0
+        ZR = 1.0
+        for K in range(NM):
+            ZR = ZR * (A + K) * (B + K) / ((K + 1.) * (C + K)) * Z
+            ZHF = ZHF + ZR
+    elif L5 or L6:
+        if (L5):
+            NM = np.round(abs(C - A))
+        if (L6):
+            NM = np.round(abs(C - B))
+        ZHF = 1.0 + 0j
+        ZR = 1.0 + 0j
+        for K in range(NM):
+            ZR *= (C - A + K) * (C - B + K) / ((K + 1.) * (C + K)) * Z
+            ZHF = ZHF + ZR
+        ZHF = (1.0 - Z) ** (C - A - B) * ZHF
+    elif (A0 <= 1.0):
+        if (X < 0.0):
+            Z1 = Z / (Z - 1.0)
+            if (C > A and B < A and B > 0.0):
+                A = BB
+                B = AA
+
+            ZC0 = 1.0 / ((1.0 - Z) ** A)
+            ZHF = 1.0 + 0j
+            ZR0 = 1.0 + 0j
+            ZW = 0
+            for K in range(500):
+                ZR0 *= (A + K) * (C - B + K) / ((K + 1.0) * (C + K)) * Z1
+                ZHF += ZR0
+                if (abs(ZHF - ZW) < abs(ZHF) * EPS):
+                    break
+                ZW = ZHF
+            ZHF = ZC0 * ZHF
+        elif (A0 >= 0.90):
+            ZW = 0.0
+            GM = 0.0
+            MCAB = np.round(C - A - B)
+            if (abs(C - A - B - MCAB) < EPS):
+                M = int(np.round(C - A - B))
+                GA = gamma(A)
+                GB = gamma(B)
+                GC = gamma(C)
+                GAM = gamma(A + M)
+                GBM = gamma(B + M)
+                PA = psi(A)
+                PB = psi(B)
+                if (M != 0):
+                    GM = 1.0
+                for j in range(1, abs(M)):
+                    GM *= j
+                RM = 1.0
+                for j in range(1, abs(M) + 1):  # DO 35 J=1,abs(M)
+                    RM *= j
+                ZF0 = 1.0
+                ZR0 = 1.0
+                ZR1 = 1.0
+                SP0 = 0.0
+                SP = 0.0
+                if (M >= 0):
+                    ZC0 = GM * GC / (GAM * GBM)
+                    ZC1 = -GC * (Z - 1.0) ** M / (GA * GB * RM)
+                    for K in range(1, M):
+                        ZR0 = ZR0 * \
+                            (A + K - 1.) * (B + K - 1.) / \
+                            (K * (K - M)) * (1. - Z)
+                        ZF0 = ZF0 + ZR0
+                    for K in range(M):
+                        SP0 = SP0 + 1.0 / \
+                            (A + K) + 1.0 / (B + K) - 1. / (K + 1.)
+                    ZF1 = PA + PB + SP0 + 2.0 * EL + np.log(1.0 - Z)
+                    for K in range(1, 501):
+                        SP = SP + \
+                            (1.0 - A) / (K * (A + K - 1.0)) + (
+                                1.0 - B) / (K * (B + K - 1.0))
+                        SM = 0.0
+                        for J in range(1, M):
+                            SM += (1.0 - A) / (
+                                (J + K) * (A + J + K - 1.0)) + 1.0 / (B + J + K - 1.0)
+
+                        ZP = PA + PB + 2.0 * EL + SP + SM + np.log(1.0 - Z)
+                        ZR1 = ZR1 * \
+                            (A + M + K - 1.0) * (B + M + K - 1.0) / (
+                                K * (M + K)) * (1.0 - Z)
+                        ZF1 = ZF1 + ZR1 * ZP
+                        if (abs(ZF1 - ZW) < abs(ZF1) * EPS):
+                            break
+                        ZW = ZF1
+                    ZHF = ZF0 * ZC0 + ZF1 * ZC1
+                elif (M < 0):
+                    M = -M
+                    ZC0 = GM * GC / (GA * GB * (1.0 - Z) ** M)
+                    ZC1 = -(-1) ** M * GC / (GAM * GBM * RM)
+                    for K in range(1, M):
+                        ZR0 = ZR0 * \
+                            (A - M + K - 1.0) * (B - M + K - 1.0) / (
+                                K * (K - M)) * (1.0 - Z)
+                        ZF0 = ZF0 + ZR0
+                    for K in range(1, M + 1):
+                        SP0 = SP0 + 1.0 / K
+                    ZF1 = PA + PB - SP0 + 2.0 * EL + np.log(1.0 - Z)
+                    for K in range(1, 501):
+                        SP = SP + \
+                            (1.0 - A) / (K * (A + K - 1.0)) + (
+                                1.0 - B) / (K * (B + K - 1.0))
+                        SM = 0.0
+                        for J in range(1, M + 1):
+                            SM = SM + 1.0 / (J + K)
+                        ZP = PA + PB + 2.0 * EL + SP - SM + np.log(1.0 - Z)
+                        ZR1 = ZR1 * \
+                            (A + K - 1.) * (B + K - 1.) / \
+                            (K * (M + K)) * (1. - Z)
+                        ZF1 = ZF1 + ZR1 * ZP
+                        if (abs(ZF1 - ZW) < abs(ZF1) * EPS):
+                            break
+                        ZW = ZF1
+                    ZHF = ZF0 * ZC0 + ZF1 * ZC1
+            else:
+                GA = gamma(A)
+                GB = gamma(B)
+                GC = gamma(C)
+                GCA = gamma(C - A)
+                GCB = gamma(C - B)
+                GCAB = gamma(C - A - B)
+                GABC = gamma(A + B - C)
+                ZC0 = GC * GCAB / (GCA * GCB)
+                ZC1 = GC * GABC / (GA * GB) * (1.0 - Z) ** (C - A - B)
+                ZHF = 0 + 0j
+                ZR0 = ZC0
+                ZR1 = ZC1
+                for K in range(1, 501):
+                    ZR0 = ZR0 * \
+                        (A + K - 1.) * (B + K - 1.) / \
+                        (K * (A + B - C + K)) * (1. - Z)
+                    ZR1 = ZR1 * \
+                        (C - A + K - 1.0) * (C - B + K - 1.0) / (
+                            K * (C - A - B + K)) * (1.0 - Z)
+                    ZHF = ZHF + ZR0 + ZR1
+                    if (abs(ZHF - ZW) < abs(ZHF) * EPS):
+                        break
+                    ZW = ZHF
+                ZHF = ZHF + ZC0 + ZC1
+        else:
+            ZW = 0.0
+            Z00 = 1.0 #+ 0j
+            if (C - A < A and C - B < B):
+                Z00 = (1.0 - Z) ** (C - A - B)
+                A = C - A
+                B = C - B
+            ZHF = 1.0
+            ZR = 1.0
+            for K in range(1, 501):
+                ZR = ZR * \
+                    (A + K - 1.0) * (B + K - 1.0) / (K * (C + K - 1.0)) * Z
+                ZHF = ZHF + ZR
+                if (abs(ZHF - ZW) <= abs(ZHF) * EPS):
+                    break
+                ZW = ZHF
+            ZHF = Z00 * ZHF
+    elif (A0 > 1.0):
+        MAB = np.round(A - B)
+        if (abs(A - B - MAB) < EPS and A0 <= 1.1):
+            B = B + EPS
+        if (abs(A - B - MAB) > EPS):
+            GA = gamma(A)
+            GB = gamma(B)
+            GC = gamma(C)
+            GAB = gamma(A - B)
+            GBA = gamma(B - A)
+            GCA = gamma(C - A)
+            GCB = gamma(C - B)
+            ZC0 = GC * GBA / (GCA * GB * (-Z) ** A)
+            ZC1 = GC * GAB / (GCB * GA * (-Z) ** B)
+            ZR0 = ZC0
+            ZR1 = ZC1
+            ZHF = 0.0 + 0j
+            for K in range(1, 501):
+                ZR0 = ZR0 * (A + K - 1.0) * (A - C + K) / ((A - B + K) * K * Z)
+                ZR1 = ZR1 * (B + K - 1.0) * (B - C + K) / ((B - A + K) * K * Z)
+                ZHF = ZHF + ZR0 + ZR1
+                if (abs((ZHF - ZW) / ZHF) <= EPS):
+                    break
+                ZW = ZHF
+            ZHF = ZHF + ZC0 + ZC1
+        else:
+            if (A - B < 0.0):
+                A = BB
+                B = AA
+            CA = C - A
+            CB = C - B
+            NCA = np.round(CA)
+            NCB = np.round(CB)
+            if (abs(CA - NCA) < EPS or abs(CB - NCB) < EPS):
+                C = C + EPS
+            GA = gamma(A)
+            GC = gamma(C)
+            GCB = gamma(C - B)
+            PA = psi(A)
+            PCA = psi(C - A)
+            PAC = psi(A - C)
+            MAB = np.round(A - B + EPS)
+            ZC0 = GC / (GA * (-Z) ** B)
+            GM = gamma(A - B)
+            ZF0 = GM / GCB * ZC0
+            ZR = ZC0
+            for K in range(1, MAB):
+                ZR = ZR * (B + K - 1.0) / (K * Z)
+                T0 = A - B - K
+                G0 = gamma(T0)
+                GCBK = gamma(C - B - K)
+                ZF0 = ZF0 + ZR * G0 / GCBK
+            if (MAB == 0):
+                ZF0 = 0.0 + 0j
+            ZC1 = GC / (GA * GCB * (-Z) ** A)
+            SP = -2.0 * EL - PA - PCA
+            for J in range(1, MAB + 1):
+                SP = SP + 1.0 / J
+            ZP0 = SP + np.log(-Z)
+            SQ = 1.0
+            for J in range(1, MAB + 1):
+                SQ = SQ * (B + J - 1.0) * (B - C + J) / J
+            ZF1 = (SQ * ZP0) * ZC1
+            ZR = ZC1
+            RK1 = 1.0
+            SJ1 = 0.0
+            W0 = 0.0
+            for K in range(1, 10001):
+                ZR = ZR / Z
+                RK1 = RK1 * (B + K - 1.0) * (B - C + K) / (K * K)
+                RK2 = RK1
+                for J in range(K + 1, K + MAB + 1):
+                    RK2 = RK2 * (B + J - 1.0) * (B - C + J) / J
+                SJ1 = SJ1 + \
+                    (A - 1.0) / (K * (A + K - 1.0)) + \
+                    (A - C - 1.0) / (K * (A - C + K - 1.0))
+                SJ2 = SJ1
+                for J in range(K + 1, K + MAB + 1):
+                    SJ2 = SJ2 + 1.0 / J
+                ZP = -2.0 * EL - PA - PAC + SJ2 - 1.0 / \
+                    (K + A - C) - PI / np.tan(PI * (K + A - C)) + np.log(-Z)
+                ZF1 = ZF1 + RK2 * ZR * ZP
+                WS = abs(ZF1)
+                if (abs((WS - W0) / WS) < EPS):
+                    break
+                W0 = WS
+            ZHF = ZF0 + ZF1
+    A = AA
+    B = BB
+    if (K > 150):
+        warnings.warn('Warning! You should check the accuracy')
+    return ZHF
+
+# def hypgf(a, b, c, x, abseps=0, releps=1e-13, kmax=10000):
+#     '''HYPGF  Hypergeometric function F(a,b,c,x)
+#
+#     CALL:   [y ,abserr] = hypgf(a,b,c,x,abseps,releps)
+#
+#     y = F(a,b,c,x)
+#     abserr = absolute error estimate
+#     a,b,c,x = input parameters
+#     abseps  = requested absolute error
+#     releps  = requested relative error
+#
+#     HYPGF calculates one solution to Gauss's hypergeometric differential
+#     equation:
+#
+#     x*(1-x)Y''(x)+[c-(a+b+1)*x]*Y'(x)-a*b*Y(x) = 0
+#     where
+#     F(a,b,c,x) = Y1(x) = 1 + a*b*x/c + a*(a+1)*b*(b+1)*x^2/(c*(c+1))+....
+#
+#
+#     Many elementary functions are special cases of F(a,b,c,x):
+#     1/(1-x) = F(1,1,1,x) = F(1,b,b,x) = F(a,1,a,x)
+#     (1+x)^n = F(-n,b,b,-x)
+#     atan(x) = x*F(.5,1,1.5,-x^2)
+#     asin(x) = x*F(.5,.5,1.5,x^2)
+#     log(x)  = x*F(1,1,2,-x)
+#     log(1+x)-log(1-x) = 2*x*F(.5,1,1.5,x^2)
+#
+#     NOTE: only real x, abs(x) < 1 and c~=0,-1,-2,... are allowed.
+#
+#     Examples:
+#     x = linspace(-.99,.99)';
+#     [Sn1,err1] = hypgf(1,1,1,x)
+#     plot(x,abs(Sn1-1./(1-x)),'b',x,err1,'r'),set(gca,'yscale','log')
+#     [Sn2,err2] = hypgf(.5,.5,1.5,x.^2);
+#     plot(x,abs(x.*Sn2-asin(x)),'b',x,abs(x.*err2),'r'),set(gca,'yscale','log')
+#
+#
+#     Reference:
+#     ---------
+#     Kreyszig, Erwin (1988)
+#     Advanced engineering mathematics
+#     John Wiley & Sons, sixth edition, pp 204.
+#     '''
+#     csize = common_shape(x, a, b, c)
+#     kmin = 2
+#     fsum = np.zeros(csize)
+#     delta = np.zeros(csize)
+#     err = np.zeros(csize)
+#
+#     ok = ~((np.round(c) == c & c <= 0) | np.abs(x) > 1)
+#     if np.any(~ok):
+#         warnings.warn('HYPGF', 'Illegal input: c = 0,-1,-2,... or abs(x)>1')
+#         fsum[~ok] = np.NaN
+#         err[~ok] = np.NaN
+#
+#     k0=find(ok & abs(x)==1);
+#     if any(k0)
+#         cmab = c(k0)-a(k0)-b(k0);
+#         fsum(k0) = exp(gammaln(c(k0))+gammaln(cmab)-...
+#                gammaln(c(k0)-a(k0))-gammaln(c(k0)-b(k0)));
+#         err(k0) = eps;
+#         k00 = find(real(cmab)<=0);
+#         if any(k00)
+#           err(k0(k00)) = nan;
+#           fsum(k0(k00)) = nan;
+#         end
+#     end
+#     k=find(ok & abs(x)<1);
+#       if any(k),
+#         delta(k) = ones(size(k));
+#         fsum(k)  = delta(k);
+#
+#         k1 = k;
+#         E = cell(1,3);
+#         E{3} = fsum(k);
+#         converge = 'n';
+#         for  ix=0:Kmax-1,
+#           delta(k1) = delta(k1).*((a(k1)+ix)./(ix+1)).*((b(k1)+ix)./(c(k1)+ ix)).*x(k1);
+#           fsum(k1) = fsum(k1)+delta(k1);
+#
+#           E(1:2) = E(2:3);
+#           E{3}   = fsum(k1);
+#
+#           if ix>Kmin
+#         if useDEA,
+#           [Sn, err(k1)] = dea3(E{:});
+#           k00 = find((abs(err(k1))) <= max(absEps,abs(relEps.*fsum(k1))));
+#           if any(k00)
+#             fsum(k1(k00)) = Sn(k00);
+#           end
+#           if (ix==Kmax-1)
+#             fsum(k1) = Sn;
+#           end
+#           k0 = (find((abs(err(k1))) > max(absEps,abs(relEps.*fsum(k1)))));
+#            if any(k0),% compute more terms
+#              %nk=length(k0);%# of values we have to compute again
+#              E{2} = E{2}(k0);
+#              E{3} = E{3}(k0);
+#            else
+#              converge='y';
+#              break;
+#            end
+#         else
+#           err(k1) = 10*abs(delta(k1));
+#            k0 = (find((abs(err(k1))) > max(absEps,abs(relEps.* ...
+#                                 fsum(k1)))));
+#            if any(k0),% compute more terms
+#              %nk=length(k0);%# of values we have to compute again
+#            else
+#              converge='y';
+#              break;
+#            end
+#         end
+#         k1 = k1(k0);
+#           end
+#         end
+#         if ~strncmpi(converge,'y',1)
+#           disp(sprintf('#%d values did not converge',length(k1)))
+#         end
+#       end
+#       %ix
+#       return
+
+
 def nextpow2(x):
     '''
     Return next higher power of 2
@@ -1761,8 +2357,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)
 
@@ -1777,7 +2371,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
 
 
@@ -1785,7 +2379,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)
@@ -1807,7 +2400,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()
 
@@ -1867,125 +2460,6 @@ def cart2polar(x, y, z=None):
         return t, r, z
 
 
-def meshgrid(*xi, **kwargs):
-    """
-    Return coordinate matrices from one or more coordinate vectors.
-
-    Make N-D coordinate arrays for vectorized evaluations of
-    N-D scalar/vector fields over N-D grids, given
-    one-dimensional coordinate arrays x1, x2,..., xn.
-
-    Parameters
-    ----------
-    x1, x2,..., xn : array_like
-        1-D arrays representing the coordinates of a grid.
-    indexing : 'xy' or 'ij' (optional)
-        cartesian ('xy', default) or matrix ('ij') indexing of output
-    sparse : True or False (default) (optional)
-         If True a sparse grid is returned in order to conserve memory.
-    copy : True (default) or False (optional)
-        If False a view into the original arrays are returned in order to
-        conserve memory
-
-    Returns
-    -------
-    X1, X2,..., XN : ndarray
-        For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` ,
-        return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij'
-        or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy'
-        with the elements of `xi` repeated to fill the matrix along
-        the first dimension for `x1`, the second for `x2` and so on.
-
-    See Also
-    --------
-    index_tricks.mgrid : Construct a multi-dimensional "meshgrid"
-                     using indexing notation.
-    index_tricks.ogrid : Construct an open multi-dimensional "meshgrid"
-                     using indexing notation.
-
-    Examples
-    --------
-    >>> 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)
-    """
-    copy_ = kwargs.get('copy', True)
-    args = atleast_1d(*xi)
-    if not isinstance(args, list):
-        if args.size > 0:
-            return args.copy() if copy_ else args
-        else:
-            raise TypeError('meshgrid() take 1 or more arguments (0 given)')
-
-    sparse = kwargs.get('sparse', False)
-    indexing = kwargs.get('indexing', 'xy')  # 'ij'
-
-    ndim = len(args)
-    s0 = (1,) * ndim
-    output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::])
-              for i, x in enumerate(args)]
-
-    shape = [x.size for x in output]
-
-    if indexing == 'xy':
-        # switch first and second axis
-        output[0].shape = (1, -1) + (1,) * (ndim - 2)
-        output[1].shape = (-1, 1) + (1,) * (ndim - 2)
-        shape[0], shape[1] = shape[1], shape[0]
-
-    if sparse:
-        if copy_:
-            return [x.copy() for x in output]
-        else:
-            return output
-    else:
-        # Return the full N-D matrix (not only the 1-D vector)
-        if copy_:
-            mult_fact = ones(shape, dtype=int)
-            return [x * mult_fact for x in output]
-        else:
-            return broadcast_arrays(*output)
-
-
 def ndgrid(*args, **kwargs):
     """
     Same as calling meshgrid with indexing='ij' (see meshgrid for
@@ -2059,8 +2533,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
@@ -2140,8 +2613,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):
@@ -2165,7 +2636,7 @@ def tranproc(x, f, x0, *xi):
     if N > 0:
         y = [y0]
         hn = xo[1] - xo[0]
-        if hn ** N < sqrt(eps):
+        if hn ** N < sqrt(_EPS):
             msg = ('Numerical problems may occur for the derivatives in ' +
             'tranproc.\nThe sampling of the transformation may be too small.')
             warnings.warn(msg)
@@ -2602,5 +3073,33 @@ def test_docstrings():
     import doctest
     doctest.testmod()
 
+
+def test_hyp2f1():
+    #     1/(1-x) = F(1,1,1,x) = F(1,b,b,x) = F(a,1,a,x)
+#     (1+x)^n = F(-n,b,b,-x)
+#     atan(x) = x*F(.5,1,1.5,-x^2)
+#     asin(x) = x*F(.5,.5,1.5,x^2)
+#     log(x)  = x*F(1,1,2,-x)
+#     log(1+x)-log(1-x) = 2*x*F(.5,1,1.5,x^2)
+    x = linspace(0., .7, 20)
+    y = hyp2f1_taylor(-1, -4, 1, .9)
+    y2 = hygfz(-1, -4, 1, .9)
+    y3 = hygfz(5, -300, 10, 0.5)
+    y4 = hyp2f1_taylor(5, -300, 10, 0.5)
+    #y = hyp2f1(0.1, 0.2, 0.3, 0.5)
+    #y = hyp2f1(1, 1.5, 3, -4 +3j)
+    #y = hyp2f1(5, 7.5, 2.5, 5)
+#     fun = lambda x : 1./(1-x)
+#     x = .99
+#     y = hyp2f1(1,1,1,x)
+#     print(y-fun(x))
+#
+    plt = plotbackend
+    plt.interactive(False)
+    plt.semilogy(x, np.abs(y- 1. / (1 - x)) + 1e-20, 'r')
+    plt.show()
+
+
 if __name__ == "__main__":
-    test_docstrings()
+    #test_docstrings()
+    test_hyp2f1()

pywafo/src/wafo/source/rind2007/krobovmod.mod

@@ -1,5 +1,5 @@
-GFORTRAN module version '4' created from intmodule.f on Sat May 05 23:15:41 2012
+GFORTRAN module version '4' created from intmodule.f on Fri Apr 05 14:43:34 2013
-MD5:eb0327a40d874f78d04c89aa93e323f2 -- If you edit this, you'll get what you deserve.
+MD5:99db0c86db329df2a1ee0bbf67b9ec99 -- If you edit this, you'll get what you deserve.
 
 (() () () () () () () () () () () () () () () () () () () () () () () ()
 () () ())
@@ -17,30 +17,30 @@ MD5:eb0327a40d874f78d04c89aa93e323f2 -- If you edit this, you'll get what you de
 (2 'krobov' 'krobovmod' 'krobov' 1 ((PROCEDURE UNKNOWN-INTENT
 MODULE-PROC DECL UNKNOWN 0 0 SUBROUTINE GENERIC) (UNKNOWN 0 0 0 UNKNOWN
 ()) 3 0 (4 5 6 7 8 9 10 11 12) () 0 () () () 0 0)
-7 'functn' '' 'functn' 3 ((PROCEDURE UNKNOWN-INTENT UNKNOWN-PROC BODY
+4 'ndim' '' 'ndim' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
-UNKNOWN 0 0 DUMMY FUNCTION ALWAYS_EXPLICIT) (REAL 8 0 0 REAL ()) 13 0 (
+DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
-14 15) () 7 () () () 0 0)
-5 'minvls' '' 'minvls' 3 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0
-0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 8 'abseps' '' 'abseps' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+6 'maxvls' '' 'maxvls' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
+9 'releps' '' 'releps' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+11 'finest' '' 'finest' 3 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0
+0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
 12 'inform' '' 'inform' 3 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0
 0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
+5 'minvls' '' 'minvls' 3 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0
+0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 10 'abserr' '' 'abserr' 3 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0
 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-11 'finest' '' 'finest' 3 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0
+7 'functn' '' 'functn' 3 ((PROCEDURE UNKNOWN-INTENT UNKNOWN-PROC BODY
-0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+UNKNOWN 0 0 DUMMY FUNCTION ALWAYS_EXPLICIT) (REAL 8 0 0 REAL ()) 13 0 (
-4 'ndim' '' 'ndim' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+14 15) () 7 () () () 0 0)
-DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
-9 'releps' '' 'releps' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
-DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-6 'maxvls' '' 'maxvls' 3 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
-DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
-14 'n' '' 'n' 13 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
-INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 15 'z' '' 'z' 13 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
 DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (1 ASSUMED_SHAPE (CONSTANT
 (INTEGER 4 0 0 INTEGER ()) 0 '1') ()) 0 () () () 0 0)
+14 'n' '' 'n' 13 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
+INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 )
 
 ('krobov' 0 2)

pywafo/src/wafo/source/rind2007/rcrudemod.mod

@@ -1,5 +1,5 @@
-GFORTRAN module version '4' created from intmodule.f on Sat May 05 23:15:40 2012
+GFORTRAN module version '4' created from intmodule.f on Fri Apr 05 14:43:34 2013
-MD5:f628260304c0d5215e1ef95941599430 -- If you edit this, you'll get what you deserve.
+MD5:c88c5a15c480306fb971bd1e5ced587e -- If you edit this, you'll get what you deserve.
 
 (() () () () () () () () () () () () () () () () () () () () () () () ()
 () () ())
@@ -17,6 +17,14 @@ MD5:f628260304c0d5215e1ef95941599430 -- If you edit this, you'll get what you de
 (2 'ranmc' 'rcrudemod' 'ranmc' 1 ((PROCEDURE UNKNOWN-INTENT MODULE-PROC
 DECL UNKNOWN 0 0 SUBROUTINE GENERIC) (UNKNOWN 0 0 0 UNKNOWN ()) 3 0 (4 5
 6 7 8 9 10 11) () 0 () () () 0 0)
+8 'releps' '' 'releps' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
+UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+9 'error' '' 'error' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
+UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+10 'value' '' 'value' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
+UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+11 'inform' '' 'inform' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
+UNKNOWN 0 0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 4 'n' '' 'n' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN UNKNOWN 0
 0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 5 'maxpts' '' 'maxpts' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
@@ -26,14 +34,6 @@ UNKNOWN 0 0 DUMMY FUNCTION ALWAYS_EXPLICIT) (REAL 8 0 0 REAL ()) 12 0 (
 13 14) () 6 () () () 0 0)
 7 'abseps' '' 'abseps' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
 UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-8 'releps' '' 'releps' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
-UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-9 'error' '' 'error' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
-UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-10 'value' '' 'value' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
-UNKNOWN 0 0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-11 'inform' '' 'inform' 3 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
-UNKNOWN 0 0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 13 'n' '' 'n' 12 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY) (
 INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 14 'z' '' 'z' 12 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0

pywafo/src/wafo/source/rind2007/rind71mod.mod

@@ -1,5 +1,5 @@
-GFORTRAN module version '4' created from rind71mod.f on Mon Feb 18 02:58:35 2013
+GFORTRAN module version '4' created from rind71mod.f on Fri Apr 05 14:43:37 2013
-MD5:520dd65f929350d1434842f22f38b888 -- If you edit this, you'll get what you deserve.
+MD5:c5460e9301460ce17aef8031cd82ad57 -- If you edit this, you'll get what you deserve.
 
 (() () () () () () () () () () () () () () () () () () () () () () ()
 () () () ())
@@ -27,16 +27,26 @@ UNKNOWN ()) 10 0 (11 12 13 14 15 16 17 18) () 0 () () () 0 0)
 4 'setdata' 'rind71mod' 'setdata' 1 ((PROCEDURE UNKNOWN-INTENT
 MODULE-PROC DECL UNKNOWN 0 0 SUBROUTINE GENERIC) (UNKNOWN 0 0 0 UNKNOWN
 ()) 19 0 (20 21 22 23 24 25 26 27 28) () 0 () () () 0 0)
-25 'dnit' '' 'dnit' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+27 'dnint' '' 'dnint' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
+28 'dxsplt' '' 'dxsplt' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
+0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+20 'method' '' 'method' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
+0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
+23 'dreps' '' 'dreps' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
 24 'deps2' '' 'deps2' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+25 'dnit' '' 'dnit' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 26 'dxc' '' 'dxc' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-27 'dnint' '' 'dnint' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+21 'scale' '' 'scale' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+22 'depss' '' 'depss' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
+DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+9 'speed' '' 'speed' 8 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
-28 'dxsplt' '' 'dxsplt' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
-0 DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
 7 'array' '' 'array' 6 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
 UNKNOWN 0 0 DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (2
 ASSUMED_SHAPE (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') () (CONSTANT (
@@ -68,16 +78,6 @@ DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (2 ASSUMED_SHAPE (CONSTANT
 DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (2 ASSUMED_SHAPE (CONSTANT
 (INTEGER 4 0 0 INTEGER ()) 0 '1') () (CONSTANT (INTEGER 4 0 0 INTEGER ())
 0 '1') ()) 0 () () () 0 0)
-9 'speed' '' 'speed' 8 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
-DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
-20 'method' '' 'method' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
-0 DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
-21 'scale' '' 'scale' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
-DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-22 'depss' '' 'depss' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
-DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-23 'dreps' '' 'dreps' 19 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
-DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
 )
 
 ('echo' 0 2 'initdata' 0 3 'rind71' 0 5 'setdata' 0 4)

pywafo/src/wafo/source/rind2007/rindmod.mod

@@ -1,5 +1,5 @@
-GFORTRAN module version '4' created from rindmod.f on Sat May 05 23:15:44 2012
+GFORTRAN module version '4' created from rindmod.f on Fri Apr 05 14:43:35 2013
-MD5:27b48943ab247880a4203cf14574fba3 -- If you edit this, you'll get what you deserve.
+MD5:dcdbb9dedca21469ecd6ba2a3e2bf880 -- If you edit this, you'll get what you deserve.
 
 (() () () () () () () () () () () () () () () () () () () () () () ()
 () () () ())
@@ -52,12 +52,6 @@ UNKNOWN ()) 15 0 (16 17 18 19 20 21 22 23 24 25 26) () 0 () () () 0 0)
 MODULE-PROC DECL UNKNOWN 0 0 SUBROUTINE GENERIC ALWAYS_EXPLICIT) (
 UNKNOWN 0 0 0 UNKNOWN ()) 27 0 (28 29 30 31 32 33 34 35 36 37) () 0 () ()
 () 0 0)
-28 'method' '' 'method' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
-0 OPTIONAL DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
-29 'xcscale' '' 'xcscale' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN
-0 0 OPTIONAL DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-30 'abseps' '' 'abseps' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
-0 OPTIONAL DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
 31 'releps' '' 'releps' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
 0 OPTIONAL DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
 32 'coveps' '' 'coveps' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
@@ -72,6 +66,12 @@ OPTIONAL DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 0 0 OPTIONAL DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
 37 'nc1c2' '' 'nc1c2' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
 OPTIONAL DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
+28 'method' '' 'method' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
+0 OPTIONAL DUMMY) (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
+29 'xcscale' '' 'xcscale' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN
+0 0 OPTIONAL DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+30 'abseps' '' 'abseps' 27 ((VARIABLE IN UNKNOWN-PROC UNKNOWN UNKNOWN 0
+0 OPTIONAL DUMMY) (REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
 16 'vals' '' 'vals' 15 ((VARIABLE OUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0
 DIMENSION DUMMY) (REAL 8 0 0 REAL ()) 0 0 () (1 ASSUMED_SHAPE (CONSTANT
 (INTEGER 4 0 0 INTEGER ()) 0 '1') ()) 0 () () () 0 0)

pywafo/src/wafo/source/rind2007/swapmod.mod

@@ -1,5 +1,5 @@
-GFORTRAN module version '4' created from swapmod.f on Sat May 05 23:15:42 2012
+GFORTRAN module version '4' created from swapmod.f on Fri Apr 05 14:43:34 2013
-MD5:52275e19413dc7ab9d6082dbb7b7af80 -- If you edit this, you'll get what you deserve.
+MD5:d3f134c81002cd5f6cec09ebff3e336f -- If you edit this, you'll get what you deserve.
 
 (() () () () () () () () () () () () () () () () () () () () () () () ()
 () () ())
@@ -25,6 +25,12 @@ DECL UNKNOWN 0 0 SUBROUTINE) (UNKNOWN 0 0 0 UNKNOWN ()) 11 0 (12 13) ()
 0 () () () 0 0)
 14 'swapmod' 'swapmod' 'swapmod' 1 ((MODULE UNKNOWN-INTENT UNKNOWN-PROC
 UNKNOWN UNKNOWN 0 0) (UNKNOWN 0 0 0 UNKNOWN ()) 0 0 () () 0 () () () 0 0)
+12 'a' '' 'a' 11 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
+(REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+13 'b' '' 'b' 11 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
+(REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
+9 'a' '' 'a' 8 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
+(INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 10 'b' '' 'b' 8 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
 (INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 6 'a' '' 'a' 5 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
@@ -33,12 +39,6 @@ UNKNOWN UNKNOWN 0 0) (UNKNOWN 0 0 0 UNKNOWN ()) 0 0 () () 0 () () () 0 0)
 7 'b' '' 'b' 5 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
 (CHARACTER 1 0 0 CHARACTER ((CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')))
 0 0 () () 0 () () () 0 0)
-12 'a' '' 'a' 11 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
-(REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-13 'b' '' 'b' 11 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
-(REAL 8 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-9 'a' '' 'a' 8 ((VARIABLE INOUT UNKNOWN-PROC UNKNOWN UNKNOWN 0 0 DUMMY)
-(INTEGER 4 0 0 INTEGER ()) 0 0 () () 0 () () () 0 0)
 )
 
 ('swap_c' 0 2 'swap_i' 0 3 'swap_r' 0 4 'swapmod' 0 14)

pywafo/src/wafo/stats/_binned_statistic.py

@@ -0,0 +1,402 @@
+from __future__ import division, print_function, absolute_import
+
+import warnings
+
+import numpy as np
+from scipy.lib.six import callable
+
+
+def binned_statistic(x, values, statistic='mean',
+                     bins=10, range=None):
+    """
+    Compute a binned statistic for a set of data.
+
+    This is a generalization of a histogram function.  A histogram divides
+    the space into bins, and returns the count of the number of points in
+    each bin.  This function allows the computation of the sum, mean, median,
+    or other statistic of the values within each bin.
+
+    .. versionadded:: 0.11.0
+
+    Parameters
+    ----------
+    x : array_like
+        A sequence of values to be binned.
+    values : array_like
+        The values on which the statistic will be computed.  This must be
+        the same shape as `x`.
+    statistic : string or callable, optional
+        The statistic to compute (default is 'mean').
+        The following statistics are available:
+
+          * 'mean' : compute the mean of values for points within each bin.
+            Empty bins will be represented by NaN.
+          * 'median' : compute the median of values for points within each
+            bin. Empty bins will be represented by NaN.
+          * 'count' : compute the count of points within each bin.  This is
+            identical to an unweighted histogram.  `values` array is not
+            referenced.
+          * 'sum' : compute the sum of values for points within each bin.
+            This is identical to a weighted histogram.
+          * function : a user-defined function which takes a 1D array of
+            values, and outputs a single numerical statistic. This function
+            will be called on the values in each bin.  Empty bins will be
+            represented by function([]), or NaN if this returns an error.
+
+    bins : int or sequence of scalars, optional
+        If `bins` is an int, it defines the number of equal-width
+        bins in the given range (10, by default). If `bins` is a sequence,
+        it defines the bin edges, including the rightmost edge, allowing
+        for non-uniform bin widths.
+    range : (float, float) or [(float, float)], optional
+        The lower and upper range of the bins.  If not provided, range
+        is simply ``(x.min(), x.max())``.  Values outside the range are
+        ignored.
+
+    Returns
+    -------
+    statistic : array
+        The values of the selected statistic in each bin.
+    bin_edges : array of dtype float
+        Return the bin edges ``(length(statistic)+1)``.
+    binnumber : 1-D ndarray of ints
+        This assigns to each observation an integer that represents the bin
+        in which this observation falls. Array has the same length as values.
+
+    See Also
+    --------
+    numpy.histogram, binned_statistic_2d, binned_statistic_dd
+
+    Notes
+    -----
+    All but the last (righthand-most) bin is half-open.  In other words, if
+    `bins` is::
+
+      [1, 2, 3, 4]
+
+    then the first bin is ``[1, 2)`` (including 1, but excluding 2) and the
+    second ``[2, 3)``.  The last bin, however, is ``[3, 4]``, which *includes*
+    4.
+
+    Examples
+    --------
+    >>> stats.binned_statistic([1, 2, 1, 2, 4], np.arange(5), statistic='mean',
+    ... bins=3)
+    (array([ 1.,  2.,  4.]), array([ 1.,  2.,  3.,  4.]), array([1, 2, 1, 2, 3]))
+
+    >>> stats.binned_statistic([1, 2, 1, 2, 4], np.arange(5), statistic='mean', bins=3)
+    (array([ 1.,  2.,  4.]), array([ 1.,  2.,  3.,  4.]), array([1, 2, 1, 2, 3]))
+
+    """
+    try:
+        N = len(bins)
+    except TypeError:
+        N = 1
+
+    if N != 1:
+        bins = [np.asarray(bins, float)]
+
+    if range is not None:
+        if len(range) == 2:
+            range = [range]
+
+    medians, edges, xy = binned_statistic_dd([x], values, statistic,
+                                             bins, range)
+
+    return medians, edges[0], xy
+
+
+def binned_statistic_2d(x, y, values, statistic='mean',
+                        bins=10, range=None):
+    """
+    Compute a bidimensional binned statistic for a set of data.
+
+    This is a generalization of a histogram2d function.  A histogram divides
+    the space into bins, and returns the count of the number of points in
+    each bin.  This function allows the computation of the sum, mean, median,
+    or other statistic of the values within each bin.
+
+    .. versionadded:: 0.11.0
+
+    Parameters
+    ----------
+    x : (N,) array_like
+        A sequence of values to be binned along the first dimension.
+    y : (M,) array_like
+        A sequence of values to be binned along the second dimension.
+    values : (N,) array_like
+        The values on which the statistic will be computed.  This must be
+        the same shape as `x`.
+    statistic : string or callable, optional
+        The statistic to compute (default is 'mean').
+        The following statistics are available:
+
+          * 'mean' : compute the mean of values for points within each bin.
+            Empty bins will be represented by NaN.
+          * 'median' : compute the median of values for points within each
+            bin. Empty bins will be represented by NaN.
+          * 'count' : compute the count of points within each bin.  This is
+            identical to an unweighted histogram.  `values` array is not
+            referenced.
+          * 'sum' : compute the sum of values for points within each bin.
+            This is identical to a weighted histogram.
+          * function : a user-defined function which takes a 1D array of
+            values, and outputs a single numerical statistic. This function
+            will be called on the values in each bin.  Empty bins will be
+            represented by function([]), or NaN if this returns an error.
+
+    bins : int or [int, int] or array-like or [array, array], optional
+        The bin specification:
+
+          * the number of bins for the two dimensions (nx=ny=bins),
+          * the number of bins in each dimension (nx, ny = bins),
+          * the bin edges for the two dimensions (x_edges = y_edges = bins),
+          * the bin edges in each dimension (x_edges, y_edges = bins).
+
+    range : (2,2) array_like, optional
+        The leftmost and rightmost edges of the bins along each dimension
+        (if not specified explicitly in the `bins` parameters):
+        [[xmin, xmax], [ymin, ymax]]. All values outside of this range will be
+        considered outliers and not tallied in the histogram.
+
+    Returns
+    -------
+    statistic : (nx, ny) ndarray
+        The values of the selected statistic in each two-dimensional bin
+    xedges : (nx + 1) ndarray
+        The bin edges along the first dimension.
+    yedges : (ny + 1) ndarray
+        The bin edges along the second dimension.
+    binnumber : 1-D ndarray of ints
+        This assigns to each observation an integer that represents the bin
+        in which this observation falls. Array has the same length as `values`.
+
+    See Also
+    --------
+    numpy.histogram2d, binned_statistic, binned_statistic_dd
+
+    """
+
+    # This code is based on np.histogram2d
+    try:
+        N = len(bins)
+    except TypeError:
+        N = 1
+
+    if N != 1 and N != 2:
+        xedges = yedges = np.asarray(bins, float)
+        bins = [xedges, yedges]
+
+    medians, edges, xy = binned_statistic_dd([x, y], values, statistic,
+                                             bins, range)
+
+    return medians, edges[0], edges[1], xy
+
+
+def binned_statistic_dd(sample, values, statistic='mean',
+                        bins=10, range=None):
+    """
+    Compute a multidimensional binned statistic for a set of data.
+
+    This is a generalization of a histogramdd function.  A histogram divides
+    the space into bins, and returns the count of the number of points in
+    each bin.  This function allows the computation of the sum, mean, median,
+    or other statistic of the values within each bin.
+
+    .. versionadded:: 0.11.0
+
+    Parameters
+    ----------
+    sample : array_like
+        Data to histogram passed as a sequence of D arrays of length N, or
+        as an (N,D) array.
+    values : array_like
+        The values on which the statistic will be computed.  This must be
+        the same shape as x.
+    statistic : string or callable, optional
+        The statistic to compute (default is 'mean').
+        The following statistics are available:
+
+          * 'mean' : compute the mean of values for points within each bin.
+            Empty bins will be represented by NaN.
+          * 'median' : compute the median of values for points within each
+            bin. Empty bins will be represented by NaN.
+          * 'count' : compute the count of points within each bin.  This is
+            identical to an unweighted histogram.  `values` array is not
+            referenced.
+          * 'sum' : compute the sum of values for points within each bin.
+            This is identical to a weighted histogram.
+          * function : a user-defined function which takes a 1D array of
+            values, and outputs a single numerical statistic. This function
+            will be called on the values in each bin.  Empty bins will be
+            represented by function([]), or NaN if this returns an error.
+
+    bins : sequence or int, optional
+        The bin specification:
+
+          * A sequence of arrays describing the bin edges along each dimension.
+          * The number of bins for each dimension (nx, ny, ... =bins)
+          * The number of bins for all dimensions (nx=ny=...=bins).
+
+    range : sequence, optional
+        A sequence of lower and upper bin edges to be used if the edges are
+        not given explicitely in `bins`. Defaults to the minimum and maximum
+        values along each dimension.
+
+    Returns
+    -------
+    statistic : ndarray, shape(nx1, nx2, nx3,...)
+        The values of the selected statistic in each two-dimensional bin
+    edges : list of ndarrays
+        A list of D arrays describing the (nxi + 1) bin edges for each
+        dimension
+    binnumber : 1-D ndarray of ints
+        This assigns to each observation an integer that represents the bin
+        in which this observation falls. Array has the same length as values.
+
+    See Also
+    --------
+    np.histogramdd, binned_statistic, binned_statistic_2d
+
+    """
+    if type(statistic) == str:
+        if statistic not in ['mean', 'median', 'count', 'sum', 'std']:
+            raise ValueError('unrecognized statistic "%s"' % statistic)
+    elif callable(statistic):
+        pass
+    else:
+        raise ValueError("statistic not understood")
+
+    # This code is based on np.histogramdd
+    try:
+        # Sample is an ND-array.
+        N, D = sample.shape
+    except (AttributeError, ValueError):
+        # Sample is a sequence of 1D arrays.
+        sample = np.atleast_2d(sample).T
+        N, D = sample.shape
+
+    nbin = np.empty(D, int)
+    edges = D * [None]
+    dedges = D * [None]
+
+    try:
+        M = len(bins)
+        if M != D:
+            raise AttributeError('The dimension of bins must be equal '
+                                 'to the dimension of the sample x.')
+    except TypeError:
+        bins = D * [bins]
+
+    # Select range for each dimension
+    # Used only if number of bins is given.
+    if range is None:
+        smin = np.atleast_1d(np.array(sample.min(0), float))
+        smax = np.atleast_1d(np.array(sample.max(0), float))
+    else:
+        smin = np.zeros(D)
+        smax = np.zeros(D)
+        for i in np.arange(D):
+            smin[i], smax[i] = range[i]
+
+    # Make sure the bins have a finite width.
+    for i in np.arange(len(smin)):
+        if smin[i] == smax[i]:
+            smin[i] = smin[i] - .5
+            smax[i] = smax[i] + .5
+
+    # Create edge arrays
+    for i in np.arange(D):
+        if np.isscalar(bins[i]):
+            nbin[i] = bins[i] + 2  # +2 for outlier bins
+            edges[i] = np.linspace(smin[i], smax[i], nbin[i] - 1)
+        else:
+            edges[i] = np.asarray(bins[i], float)
+            nbin[i] = len(edges[i]) + 1  # +1 for outlier bins
+        dedges[i] = np.diff(edges[i])
+
+    nbin = np.asarray(nbin)
+
+    # Compute the bin number each sample falls into.
+    Ncount = {}
+    for i in np.arange(D):
+        Ncount[i] = np.digitize(sample[:, i], edges[i])
+
+    # Using digitize, values that fall on an edge are put in the right bin.
+    # For the rightmost bin, we want values equal to the right
+    # edge to be counted in the last bin, and not as an outlier.
+    for i in np.arange(D):
+        # Rounding precision
+        decimal = int(-np.log10(dedges[i].min())) + 6
+        # Find which points are on the rightmost edge.
+        on_edge = np.where(np.around(sample[:, i], decimal)
+                           == np.around(edges[i][-1], decimal))[0]
+        # Shift these points one bin to the left.
+        Ncount[i][on_edge] -= 1
+
+    # Compute the sample indices in the flattened statistic matrix.
+    ni = nbin.argsort()
+    xy = np.zeros(N, int)
+    for i in np.arange(0, D - 1):
+        xy += Ncount[ni[i]] * nbin[ni[i + 1:]].prod()
+    xy += Ncount[ni[-1]]
+
+    result = np.empty(nbin.prod(), float)
+
+    if statistic == 'mean':
+        result.fill(np.nan)
+        flatcount = np.bincount(xy, None)
+        flatsum = np.bincount(xy, values)
+        a = flatcount.nonzero()
+        result[a] = flatsum[a] / flatcount[a]
+    elif statistic == 'std':
+        result.fill(0)
+        flatcount = np.bincount(xy, None)
+        flatsum = np.bincount(xy, values)
+        flatsum2 = np.bincount(xy, values ** 2)
+        a = flatcount.nonzero()
+        result[a] = np.sqrt(flatsum2[a] / flatcount[a]
+                            - (flatsum[a] / flatcount[a]) ** 2)
+    elif statistic == 'count':
+        result.fill(0)
+        flatcount = np.bincount(xy, None)
+        a = np.arange(len(flatcount))
+        result[a] = flatcount
+    elif statistic == 'sum':
+        result.fill(0)
+        flatsum = np.bincount(xy, values)
+        a = np.arange(len(flatsum))
+        result[a] = flatsum
+    elif statistic == 'median':
+        result.fill(np.nan)
+        for i in np.unique(xy):
+            result[i] = np.median(values[xy == i])
+    elif callable(statistic):
+        with warnings.catch_warnings():
+            # Numpy generates a warnings for mean/std/... with empty list
+            warnings.filterwarnings('ignore', category=RuntimeWarning)
+            old = np.seterr(invalid='ignore')
+            try:
+                null = statistic([])
+            except:
+                null = np.nan
+            np.seterr(**old)
+        result.fill(null)
+        for i in np.unique(xy):
+            result[i] = statistic(values[xy == i])
+
+    # Shape into a proper matrix
+    result = result.reshape(np.sort(nbin))
+    for i in np.arange(nbin.size):
+        j = ni.argsort()[i]
+        result = result.swapaxes(i, j)
+        ni[i], ni[j] = ni[j], ni[i]
+
+    # Remove outliers (indices 0 and -1 for each dimension).
+    core = D * [slice(1, -1)]
+    result = result[core]
+
+    if (result.shape != nbin - 2).any():
+        raise RuntimeError('Internal Shape Error')
+
+    return result, edges, xy

pywafo/src/wafo/stats/_constants.py

@@ -0,0 +1,24 @@
+"""
+Statistics-related constants.
+
+"""
+from __future__ import division, print_function, absolute_import
+
+import numpy as np
+
+
+# The smallest representable positive number such that 1.0 + _EPS != 1.0.
+_EPS = np.finfo(float).eps
+
+# The largest [in magnitude] usable floating value.
+_XMAX = np.finfo(float).machar.xmax
+
+# The smallest [in magnitude] usable floating value.
+_XMIN = np.finfo(float).machar.xmin
+
+# -special.psi(1)
+_EULER = 0.577215664901532860606512090082402431042
+
+# special.zeta(3, 1)  Apery's constant
+_ZETA3 = 1.202056903159594285399738161511449990765
+

pywafo/src/wafo/stats/_discrete_distns.py

@@ -19,10 +19,11 @@ __all__ = [
     'binom', 'bernoulli', 'nbinom', 'geom', 'hypergeom',
     'logser', 'poisson', 'planck', 'boltzmann', 'randint',
     'zipf', 'dlaplace', 'skellam'
-        ]
+]
 
 
 class binom_gen(rv_discrete):
+
     """A binomial discrete random variable.
 
     %(before_notes)s
@@ -40,6 +41,7 @@ class binom_gen(rv_discrete):
     %(example)s
 
     """
+
     def _rvs(self, n, p):
         return mtrand.binomial(n, p, self._size)
 
@@ -49,8 +51,8 @@ class binom_gen(rv_discrete):
 
     def _logpmf(self, x, n, p):
         k = floor(x)
-        combiln = (gamln(n+1) - (gamln(k+1) + gamln(n-k+1)))
+        combiln = (gamln(n + 1) - (gamln(k + 1) + gamln(n - k + 1)))
-        return combiln + special.xlogy(k, p) + special.xlog1py(n-k, -p)
+        return combiln + special.xlogy(k, p) + special.xlog1py(n - k, -p)
 
     def _pmf(self, x, n, p):
         return exp(self._logpmf(x, n, p))
@@ -66,16 +68,16 @@ class binom_gen(rv_discrete):
 
     def _ppf(self, q, n, p):
         vals = ceil(special.bdtrik(q, n, p))
-        vals1 = vals-1
+        vals1 = vals - 1
         temp = special.bdtr(vals1, n, p)
         return np.where(temp >= q, vals1, vals)
 
     def _stats(self, n, p):
-        q = 1.0-p
+        q = 1.0 - p
         mu = n * p
         var = n * p * q
-        g1 = (q-p) / sqrt(n*p*q)
+        g1 = (q - p) / sqrt(n * p * q)
-        g2 = (1.0-6*p*q)/(n*p*q)
+        g2 = (1.0 - 6 * p * q) / (n * p * q)
         return mu, var, g1, g2
 
     def _entropy(self, n, p):
@@ -87,6 +89,7 @@ binom = binom_gen(name='binom')
 
 
 class bernoulli_gen(binom_gen):
+
     """A Bernoulli discrete random variable.
 
     %(before_notes)s
@@ -105,6 +108,7 @@ class bernoulli_gen(binom_gen):
     %(example)s
 
     """
+
     def _rvs(self, p):
         return binom_gen._rvs(self, 1, p)
 
@@ -136,6 +140,7 @@ bernoulli = bernoulli_gen(b=1, name='bernoulli')
 
 
 class nbinom_gen(rv_discrete):
+
     """A negative binomial discrete random variable.
 
     %(before_notes)s
@@ -153,6 +158,7 @@ class nbinom_gen(rv_discrete):
     %(example)s
 
     """
+
     def _rvs(self, n, p):
         return mtrand.negative_binomial(n, p, self._size)
 
@@ -163,12 +169,12 @@ class nbinom_gen(rv_discrete):
         return exp(self._logpmf(x, n, p))
 
     def _logpmf(self, x, n, p):
-        coeff = gamln(n+x) - gamln(x+1) - gamln(n)
+        coeff = gamln(n + x) - gamln(x + 1) - gamln(n)
-        return coeff + n*log(p) + x*log1p(-p)
+        return coeff + n * log(p) + x * log1p(-p)
 
     def _cdf(self, x, n, p):
         k = floor(x)
-        return special.betainc(n, k+1, p)
+        return special.betainc(n, k + 1, p)
 
     def _sf_skip(self, x, n, p):
         # skip because special.nbdtrc doesn't work for 0<n<1
@@ -177,22 +183,23 @@ class nbinom_gen(rv_discrete):
 
     def _ppf(self, q, n, p):
         vals = ceil(special.nbdtrik(q, n, p))
-        vals1 = (vals-1).clip(0.0, np.inf)
+        vals1 = (vals - 1).clip(0.0, np.inf)
         temp = self._cdf(vals1, n, p)
         return np.where(temp >= q, vals1, vals)
 
     def _stats(self, n, p):
         Q = 1.0 / p
         P = Q - 1.0
-        mu = n*P
+        mu = n * P
-        var = n*P*Q
+        var = n * P * Q
-        g1 = (Q+P)/sqrt(n*P*Q)
+        g1 = (Q + P) / sqrt(n * P * Q)
-        g2 = (1.0 + 6*P*Q) / (n*P*Q)
+        g2 = (1.0 + 6 * P * Q) / (n * P * Q)
         return mu, var, g1, g2
 nbinom = nbinom_gen(name='nbinom')
 
 
 class geom_gen(rv_discrete):
+
     """A geometric discrete random variable.
 
     %(before_notes)s
@@ -210,6 +217,7 @@ class geom_gen(rv_discrete):
     %(example)s
 
     """
+
     def _rvs(self, p):
         return mtrand.geometric(p, size=self._size)
 
@@ -217,38 +225,39 @@ class geom_gen(rv_discrete):
         return (p <= 1) & (p >= 0)
 
     def _pmf(self, k, p):
-        return np.power(1-p, k-1) * p
+        return np.power(1 - p, k - 1) * p
 
     def _logpmf(self, k, p):
-        return (k-1)*log1p(-p) + log(p)
+        return (k - 1) * log1p(-p) + log(p)
 
     def _cdf(self, x, p):
         k = floor(x)
-        return -expm1(log1p(-p)*k)
+        return -expm1(log1p(-p) * k)
 
     def _sf(self, x, p):
         return np.exp(self._logsf(x, p))
 
     def _logsf(self, x, p):
         k = floor(x)
-        return k*log1p(-p)
+        return k * log1p(-p)
 
     def _ppf(self, q, p):
-        vals = ceil(log1p(-q)/log1p(-p))
+        vals = ceil(log1p(-q) / log1p(-p))
-        temp = self._cdf(vals-1, p)
+        temp = self._cdf(vals - 1, p)
-        return np.where((temp >= q) & (vals > 0), vals-1, vals)
+        return np.where((temp >= q) & (vals > 0), vals - 1, vals)
 
     def _stats(self, p):
-        mu = 1.0/p
+        mu = 1.0 / p
-        qr = 1.0-p
+        qr = 1.0 - p
         var = qr / p / p
-        g1 = (2.0-p) / sqrt(qr)
+        g1 = (2.0 - p) / sqrt(qr)
-        g2 = np.polyval([1, -6, 6], p)/(1.0-p)
+        g2 = np.polyval([1, -6, 6], p) / (1.0 - p)
         return mu, var, g1, g2
 geom = geom_gen(a=1, name='geom', longname="A geometric")
 
 
 class hypergeom_gen(rv_discrete):
+
     """A hypergeometric discrete random variable.
 
     The hypergeometric distribution models drawing objects from a bin.
@@ -298,22 +307,23 @@ class hypergeom_gen(rv_discrete):
     >>> R = hypergeom.rvs(M, n, N, size=10)
 
     """
+
     def _rvs(self, M, n, N):
-        return mtrand.hypergeometric(n, M-n, N, size=self._size)
+        return mtrand.hypergeometric(n, M - n, N, size=self._size)
 
     def _argcheck(self, M, n, N):
         cond = rv_discrete._argcheck(self, M, n, N)
         cond &= (n <= M) & (N <= M)
-        self.a = max(N-(M-n), 0)
+        self.a = max(N - (M - n), 0)
         self.b = min(n, N)
         return cond
 
     def _logpmf(self, k, M, n, N):
         tot, good = M, n
         bad = tot - good
-        return gamln(good+1) - gamln(good-k+1) - gamln(k+1) + gamln(bad+1) \
+        return gamln(good + 1) - gamln(good - k + 1) - gamln(k + 1) + \
-            - gamln(bad-N+k+1) - gamln(N-k+1) - gamln(tot+1) + gamln(tot-N+1) \
+            gamln(bad + 1) - gamln(bad - N + k + 1) - gamln(N - k + 1) - \
-            + gamln(N+1)
+            gamln(tot + 1) + gamln(tot - N + 1) + gamln(N + 1)
 
     def _pmf(self, k, M, n, N):
         # same as the following but numerically more precise
@@ -323,18 +333,19 @@ class hypergeom_gen(rv_discrete):
     def _stats(self, M, n, N):
         # tot, good, sample_size = M, n, N
         # "wikipedia".replace('N', 'M').replace('n', 'N').replace('K', 'n')
-        M, n, N = 1.*M, 1.*n, 1.*N
+        M, n, N = 1. * M, 1. * n, 1. * N
         m = M - n
-        p = n/M
+        p = n / M
-        mu = N*p
+        mu = N * p
 
-        var = m*n*N*(M - N)*1.0/(M*M*(M-1))
+        var = m * n * N * (M - N) * 1.0 / (M * M * (M - 1))
-        g1 = (m - n)*(M-2*N) / (M-2.0) * sqrt((M-1.0) / (m*n*N*(M-N)))
+        g1 = (m - n) * (M - 2 * N) / (M - 2.0) * \
+            sqrt((M - 1.0) / (m * n * N * (M - N)))
 
-        g2 = M*(M+1) - 6.*N*(M-N) - 6.*n*m
+        g2 = M * (M + 1) - 6. * N * (M - N) - 6. * n * m
-        g2 *= (M-1)*M*M
+        g2 *= (M - 1) * M * M
-        g2 += 6.*n*N*(M-N)*m*(5.*M-6)
+        g2 += 6. * n * N * (M - N) * m * (5. * M - 6)
-        g2 /= n * N * (M-N) * m * (M-2.) * (M-3.)
+        g2 /= n * N * (M - N) * m * (M - 2.) * (M - 3.)
         return mu, var, g1, g2
 
     def _entropy(self, M, n, N):
@@ -361,6 +372,7 @@ hypergeom = hypergeom_gen(name='hypergeom')
 
 # FIXME: Fails _cdfvec
 class logser_gen(rv_discrete):
+
     """A Logarithmic (Log-Series, Series) discrete random variable.
 
     %(before_notes)s
@@ -378,6 +390,7 @@ class logser_gen(rv_discrete):
     %(example)s
 
     """
+
     def _rvs(self, p):
         # looks wrong for p>0.5, too few k=1
         # trying to use generic is worse, no k=1 at all
@@ -392,21 +405,22 @@ class logser_gen(rv_discrete):
     def _stats(self, p):
         r = log1p(-p)
         mu = p / (p - 1.0) / r
-        mu2p = -p / r / (p - 1.0)**2
+        mu2p = -p / r / (p - 1.0) ** 2
-        var = mu2p - mu*mu
+        var = mu2p - mu * mu
-        mu3p = -p / r * (1.0+p) / (1.0 - p)**3
+        mu3p = -p / r * (1.0 + p) / (1.0 - p) ** 3
-        mu3 = mu3p - 3*mu*mu2p + 2*mu**3
+        mu3 = mu3p - 3 * mu * mu2p + 2 * mu ** 3
         g1 = mu3 / np.power(var, 1.5)
 
         mu4p = -p / r * (
-            1.0 / (p-1)**2 - 6*p / (p - 1)**3 + 6*p*p / (p-1)**4)
+            1.0 / (p - 1) ** 2 - 6 * p / (p - 1) ** 3 + 6 * p * p / (p - 1) ** 4)
-        mu4 = mu4p - 4*mu3p*mu + 6*mu2p*mu*mu - 3*mu**4
+        mu4 = mu4p - 4 * mu3p * mu + 6 * mu2p * mu * mu - 3 * mu ** 4
-        g2 = mu4 / var**2 - 3.0
+        g2 = mu4 / var ** 2 - 3.0
         return mu, var, g1, g2
 logser = logser_gen(a=1, name='logser', longname='A logarithmic')
 
 
 class poisson_gen(rv_discrete):
+
     """A Poisson discrete random variable.
 
     %(before_notes)s
@@ -424,11 +438,12 @@ class poisson_gen(rv_discrete):
     %(example)s
 
     """
+
     def _rvs(self, mu):
         return mtrand.poisson(mu, self._size)
 
     def _logpmf(self, k, mu):
-        Pk = k*log(mu)-gamln(k+1) - mu
+        Pk = k * log(mu) - gamln(k + 1) - mu
         return Pk
 
     def _pmf(self, k, mu):
@@ -458,6 +473,7 @@ poisson = poisson_gen(name="poisson", longname='A Poisson')
 
 
 class planck_gen(rv_discrete):
+
     """A Planck discrete exponential random variable.
 
     %(before_notes)s
@@ -475,6 +491,7 @@ class planck_gen(rv_discrete):
     %(example)s
 
     """
+
     def _argcheck(self, lambda_):
         if (lambda_ > 0):
             self.a = 0
@@ -496,16 +513,16 @@ class planck_gen(rv_discrete):
         return - expm1(-lambda_ * (k + 1))
 
     def _ppf(self, q, lambda_):
-        vals = ceil(-1.0/lambda_ * log1p(-q)-1)
+        vals = ceil(-1.0 / lambda_ * log1p(-q) - 1)
-        vals1 = (vals-1).clip(self.a, np.inf)
+        vals1 = (vals - 1).clip(self.a, np.inf)
         temp = self._cdf(vals1, lambda_)
         return np.where(temp >= q, vals1, vals)
 
     def _stats(self, lambda_):
-        mu = 1/(exp(lambda_)-1)
+        mu = 1 / (exp(lambda_) - 1)
-        var = exp(-lambda_)/(expm1(-lambda_))**2
+        var = exp(-lambda_) / (expm1(-lambda_)) ** 2
-        g1 = 2*cosh(lambda_/2.0)
+        g1 = 2 * cosh(lambda_ / 2.0)
-        g2 = 4+2*cosh(lambda_)
+        g2 = 4 + 2 * cosh(lambda_)
         return mu, var, g1, g2
 
     def _entropy(self, lambda_):
@@ -516,6 +533,7 @@ planck = planck_gen(name='planck', longname='A discrete exponential ')
 
 
 class boltzmann_gen(rv_discrete):
+
     """A Boltzmann (Truncated Discrete Exponential) random variable.
 
     %(before_notes)s
@@ -533,31 +551,33 @@ class boltzmann_gen(rv_discrete):
     %(example)s
 
     """
+
     def _pmf(self, k, lambda_, N):
-        fact = (expm1(-lambda_))/(expm1(-lambda_*N))
+        fact = (expm1(-lambda_)) / (expm1(-lambda_ * N))
-        return fact*exp(-lambda_*k)
+        return fact * exp(-lambda_ * k)
 
     def _cdf(self, x, lambda_, N):
         k = floor(x)
-        return (expm1(-lambda_*(k+1)))/(expm1(-lambda_*N))
+        return (expm1(-lambda_ * (k + 1))) / (expm1(-lambda_ * N))
 
     def _ppf(self, q, lambda_, N):
-        qnew = -q*(expm1(-lambda_*N))
+        qnew = -q * (expm1(-lambda_ * N))
-        vals = ceil(-1.0/lambda_ * log1p(-qnew)-1)
+        vals = ceil(-1.0 / lambda_ * log1p(-qnew) - 1)
-        vals1 = (vals-1).clip(0.0, np.inf)
+        vals1 = (vals - 1).clip(0.0, np.inf)
         temp = self._cdf(vals1, lambda_, N)
         return np.where(temp >= q, vals1, vals)
 
     def _stats(self, lambda_, N):
         z = exp(-lambda_)
-        zN = exp(-lambda_*N)
+        zN = exp(-lambda_ * N)
-        mu = z/(1.0-z)-N*zN/(1-zN)
+        mu = z / (1.0 - z) - N * zN / (1 - zN)
-        var = z/(1.0-z)**2 - N*N*zN/(1-zN)**2
+        var = z / (1.0 - z) ** 2 - N * N * zN / (1 - zN) ** 2
-        trm = (1-zN)/(1-z)
+        trm = (1 - zN) / (1 - z)
-        trm2 = (z*trm**2 - N*N*zN)
+        trm2 = (z * trm ** 2 - N * N * zN)
-        g1 = z*(1+z)*trm**3 - N**3*zN*(1+zN)
+        g1 = z * (1 + z) * trm ** 3 - N ** 3 * zN * (1 + zN)
-        g1 = g1 / trm2**(1.5)
+        g1 = g1 / trm2 ** (1.5)
-        g2 = z*(1+4*z+z*z)*trm**4 - N**4 * zN*(1+4*zN+zN*zN)
+        g2 = z * (1 + 4 * z + z * z) * \
+            trm ** 4 - N ** 4 * zN * (1 + 4 * zN + zN * zN)
         g2 = g2 / trm2 / trm2
         return mu, var, g1, g2
 boltzmann = boltzmann_gen(name='boltzmann',
@@ -565,6 +585,7 @@ boltzmann = boltzmann_gen(name='boltzmann',
 
 
 class randint_gen(rv_discrete):
+
     """A uniform discrete random variable.
 
     %(before_notes)s
@@ -585,6 +606,7 @@ class randint_gen(rv_discrete):
     %(example)s
 
     """
+
     def _argcheck(self, low, high):
         self.a = low
         self.b = high - 1
@@ -608,9 +630,9 @@ class randint_gen(rv_discrete):
         m2, m1 = np.asarray(high), np.asarray(low)
         mu = (m2 + m1 - 1.0) / 2
         d = m2 - m1
-        var = (d*d - 1) / 12.0
+        var = (d * d - 1) / 12.0
         g1 = 0.0
-        g2 = -6.0/5.0 * (d*d + 1.0) / (d*d - 1.0)
+        g2 = -6.0 / 5.0 * (d * d + 1.0) / (d * d - 1.0)
         return mu, var, g1, g2
 
     def _rvs(self, low, high=None):
@@ -628,6 +650,7 @@ randint = randint_gen(name='randint', longname='A discrete uniform '
 
 # FIXME: problems sampling.
 class zipf_gen(rv_discrete):
+
     """A Zipf discrete random variable.
 
     %(before_notes)s
@@ -645,6 +668,7 @@ class zipf_gen(rv_discrete):
     %(example)s
 
     """
+
     def _rvs(self, a):
         return mtrand.zipf(a, size=self._size)
 
@@ -652,7 +676,7 @@ class zipf_gen(rv_discrete):
         return a > 1
 
     def _pmf(self, k, a):
-        Pk = 1.0 / special.zeta(a, 1) / k**a
+        Pk = 1.0 / special.zeta(a, 1) / k ** a
         return Pk
 
     def _munp(self, n, a):
@@ -664,6 +688,7 @@ zipf = zipf_gen(a=1, name='zipf', longname='A Zipf')
 
 
 class dlaplace_gen(rv_discrete):
+
     """A  Laplacian discrete random variable.
 
     %(before_notes)s
@@ -681,35 +706,37 @@ class dlaplace_gen(rv_discrete):
     %(example)s
 
     """
+
     def _pmf(self, k, a):
-        return tanh(a/2.0) * exp(-a * abs(k))
+        return tanh(a / 2.0) * exp(-a * abs(k))
 
     def _cdf(self, x, a):
         k = floor(x)
         f = lambda k, a: 1.0 - exp(-a * k) / (exp(a) + 1)
-        f2 = lambda k, a: exp(a * (k+1)) / (exp(a) + 1)
+        f2 = lambda k, a: exp(a * (k + 1)) / (exp(a) + 1)
         return _lazywhere(k >= 0, (k, a), f=f, f2=f2)
 
     def _ppf(self, q, a):
         const = 1 + exp(a)
-        vals = ceil(np.where(q < 1.0 / (1 + exp(-a)), log(q*const) / a - 1,
+        vals = ceil(np.where(q < 1.0 / (1 + exp(-a)), log(q * const) / a - 1,
-                                                      -log((1-q) * const) / a))
+                             -log((1 - q) * const) / a))
         vals1 = vals - 1
         return np.where(self._cdf(vals1, a) >= q, vals1, vals)
 
     def _stats(self, a):
         ea = exp(a)
-        mu2 = 2.*ea/(ea-1.)**2
+        mu2 = 2. * ea / (ea - 1.) ** 2
-        mu4 = 2.*ea*(ea**2+10.*ea+1.) / (ea-1.)**4
+        mu4 = 2. * ea * (ea ** 2 + 10. * ea + 1.) / (ea - 1.) ** 4
-        return 0., mu2, 0., mu4/mu2**2 - 3.
+        return 0., mu2, 0., mu4 / mu2 ** 2 - 3.
 
     def _entropy(self, a):
-        return a / sinh(a) - log(tanh(a/2.0))
+        return a / sinh(a) - log(tanh(a / 2.0))
 dlaplace = dlaplace_gen(a=-np.inf,
                         name='dlaplace', longname='A discrete Laplacian')
 
 
 class skellam_gen(rv_discrete):
+
     """A  Skellam discrete random variable.
 
     %(before_notes)s
@@ -735,28 +762,29 @@ class skellam_gen(rv_discrete):
     %(example)s
 
     """
+
     def _rvs(self, mu1, mu2):
         n = self._size
         return mtrand.poisson(mu1, n) - mtrand.poisson(mu2, n)
 
     def _pmf(self, x, mu1, mu2):
         px = np.where(x < 0,
-                _ncx2_pdf(2*mu2, 2*(1-x), 2*mu1)*2,
+                      _ncx2_pdf(2 * mu2, 2 * (1 - x), 2 * mu1) * 2,
-                _ncx2_pdf(2*mu1, 2*(1+x), 2*mu2)*2)
+                      _ncx2_pdf(2 * mu1, 2 * (1 + x), 2 * mu2) * 2)
         # ncx2.pdf() returns nan's for extremely low probabilities
         return px
 
     def _cdf(self, x, mu1, mu2):
         x = floor(x)
         px = np.where(x < 0,
-                _ncx2_cdf(2*mu2, -2*x, 2*mu1),
+                      _ncx2_cdf(2 * mu2, -2 * x, 2 * mu1),
-                1-_ncx2_cdf(2*mu1, 2*(x+1), 2*mu2))
+                      1 - _ncx2_cdf(2 * mu1, 2 * (x + 1), 2 * mu2))
         return px
 
     def _stats(self, mu1, mu2):
         mean = mu1 - mu2
         var = mu1 + mu2
-        g1 = mean / sqrt((var)**3)
+        g1 = mean / sqrt((var) ** 3)
         g2 = 1 / var
         return mean, var, g1, g2
 skellam = skellam_gen(a=-np.inf, name="skellam", longname='A Skellam')

pywafo/src/wafo/stats/_distn_infrastructure.py

@@ -455,7 +455,7 @@ class rv_frozen(object):
     def __init__(self, dist, *args, **kwds):
         self.dist = dist
         args, loc, scale = dist._parse_args(*args, **kwds)
-        if len(args) == dist.numargs - 2:  # isinstance(dist, rv_continuous):
+        if isinstance(dist, rv_continuous):
             self.par = args + (loc, scale)
         else:  # rv_discrete
             self.par = args + (loc,)

pywafo/src/wafo/stats/kde_test.py

@@ -0,0 +1,15 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Dec 06 16:02:47 2011
+
+@author: pab
+"""
+import numpy as np
+import wafo.kdetools as wk
+n = 100
+x = np.sort(5*np.random.rand(1,n)-2.5, axis=-1).ravel()
+y = (np.cos(x)>2*np.random.rand(n, 1)-1).ravel()
+
+kreg = wk.KRegression(x,y)
+f = kreg(output='plotobj', title='Kernel regression', plotflag=1)
+f.plot()

pywafo/src/wafo/stats/morestats.py

@@ -8,23 +8,25 @@ import math
 import warnings
 
 import numpy as np
-from numpy import (isscalar, r_, log, sum, around, unique, asarray,
+from numpy import (isscalar, r_, log, sum, around, unique, asarray, zeros,
-     zeros, arange, sort, amin, amax, any, atleast_1d, sqrt, ceil,
+                   arange, sort, amin, amax, any, atleast_1d, sqrt, ceil,
-     floor, array, poly1d, compress, not_equal, pi, exp, ravel, angle)
+                   floor, array, poly1d, compress, not_equal, pi, exp, ravel,
+                   angle)
 from numpy.testing.decorators import setastest
 
 from scipy.lib.six import string_types
 from scipy import optimize
 from scipy import special
-from . import statlib
+from wafo.stats import statlib
-from . import stats
+from wafo.stats import stats
-from .stats import find_repeats
+from  wafo.stats.stats import find_repeats
-from . import distributions
+from wafo.stats import distributions
-from ._distn_infrastructure import rv_generic
+from wafo.stats._distn_infrastructure import rv_generic
 
 
 __all__ = ['mvsdist',
-           'bayes_mvs', 'kstat', 'kstatvar', 'probplot', 'ppcc_max', 'ppcc_plot',
+           'bayes_mvs', 'kstat', 'kstatvar', 'probplot', 'ppcc_max',
+           'ppcc_plot',
            'boxcox_llf', 'boxcox', 'boxcox_normmax', 'boxcox_normplot',
            'shapiro', 'anderson', 'ansari', 'bartlett', 'levene', 'binom_test',
            'fligner', 'mood', 'wilcoxon',
@@ -78,7 +80,8 @@ def bayes_mvs(data, alpha=0.90):
     """
     res = mvsdist(data)
     if alpha >= 1 or alpha <= 0:
-        raise ValueError("0 < alpha < 1 is required, but alpha=%s was given." % alpha)
+        raise ValueError(
+            "0 < alpha < 1 is required, but alpha=%s was given." % alpha)
     return tuple((x.mean(), x.interval(alpha)) for x in res)
 
 
@@ -134,20 +137,21 @@ def mvsdist(data):
     xbar = x.mean()
     C = x.var()
     if (n > 1000):  # gaussian approximations for large n
-        mdist = distributions.norm(loc=xbar, scale=math.sqrt(C/n))
+        mdist = distributions.norm(loc=xbar, scale=math.sqrt(C / n))
-        sdist = distributions.norm(loc=math.sqrt(C), scale=math.sqrt(C/(2.*n)))
+        sdist = distributions.norm(
-        vdist = distributions.norm(loc=C, scale=math.sqrt(2.0/n)*C)
+            loc=math.sqrt(C), scale=math.sqrt(C / (2. * n)))
+        vdist = distributions.norm(loc=C, scale=math.sqrt(2.0 / n) * C)
     else:
-        nm1 = n-1
+        nm1 = n - 1
-        fac = n*C/2.
+        fac = n * C / 2.
-        val = nm1/2.
+        val = nm1 / 2.
-        mdist = distributions.t(nm1,loc=xbar,scale=math.sqrt(C/nm1))
+        mdist = distributions.t(nm1, loc=xbar, scale=math.sqrt(C / nm1))
-        sdist = distributions.gengamma(val,-2,scale=math.sqrt(fac))
+        sdist = distributions.gengamma(val, -2, scale=math.sqrt(fac))
-        vdist = distributions.invgamma(val,scale=fac)
+        vdist = distributions.invgamma(val, scale=fac)
     return mdist, vdist, sdist
 
 
-def kstat(data,n=2):
+def kstat(data, n=2):
     """
     Return the nth k-statistic (1<=n<=4 so far).
 
@@ -197,26 +201,26 @@ def kstat(data,n=2):
     if n > 4 or n < 1:
         raise ValueError("k-statistics only supported for 1<=n<=4")
     n = int(n)
-    S = zeros(n+1,'d')
+    S = zeros(n + 1, 'd')
     data = ravel(data)
     N = len(data)
-    for k in range(1,n+1):
+    for k in range(1, n + 1):
-        S[k] = sum(data**k,axis=0)
+        S[k] = sum(data ** k, axis=0)
     if n == 1:
-        return S[1]*1.0/N
+        return S[1] * 1.0 / N
     elif n == 2:
-        return (N*S[2]-S[1]**2.0)/(N*(N-1.0))
+        return (N * S[2] - S[1] ** 2.0) / (N * (N - 1.0))
     elif n == 3:
-        return (2*S[1]**3 - 3*N*S[1]*S[2]+N*N*S[3]) / (N*(N-1.0)*(N-2.0))
+        return (2 * S[1] ** 3 - 3 * N * S[1] * S[2] + N * N * S[3]) / (N * (N - 1.0) * (N - 2.0))
     elif n == 4:
-        return (-6*S[1]**4 + 12*N*S[1]**2 * S[2] - 3*N*(N-1.0)*S[2]**2 -
+        return (-6 * S[1] ** 4 + 12 * N * S[1] ** 2 * S[2] - 3 * N * (N - 1.0) * S[2] ** 2 -
-                4*N*(N+1)*S[1]*S[3] + N*N*(N+1)*S[4]) / \
+                4 * N * (N + 1) * S[1] * S[3] + N * N * (N + 1) * S[4]) / \
-                (N*(N-1.0)*(N-2.0)*(N-3.0))
+            (N * (N - 1.0) * (N - 2.0) * (N - 3.0))
     else:
         raise ValueError("Should not be here.")
 
 
-def kstatvar(data,n=2):
+def kstatvar(data, n=2):
     """
     Returns an unbiased estimator of the variance of the k-statistic.
 
@@ -242,11 +246,11 @@ def kstatvar(data,n=2):
     data = ravel(data)
     N = len(data)
     if n == 1:
-        return kstat(data,n=2)*1.0/N
+        return kstat(data, n=2) * 1.0 / N
     elif n == 2:
-        k2 = kstat(data,n=2)
+        k2 = kstat(data, n=2)
-        k4 = kstat(data,n=4)
+        k4 = kstat(data, n=4)
-        return (2*k2*k2*N + (N-1)*k4)/(N*(N+1))
+        return (2 * k2 * k2 * N + (N - 1) * k4) / (N * (N + 1))
     else:
         raise ValueError("Only n=1 or n=2 supported.")
 
@@ -255,7 +259,7 @@ def _calc_uniform_order_statistic_medians(x):
     """See Notes section of `probplot` for details."""
     N = len(x)
     osm_uniform = np.zeros(N, dtype=np.float64)
-    osm_uniform[-1] = 0.5**(1.0 / N)
+    osm_uniform[-1] = 0.5 ** (1.0 / N)
     osm_uniform[0] = 1 - osm_uniform[-1]
     i = np.arange(2, N)
     osm_uniform[1:-1] = (i - 0.3175) / (N + 0.365)
@@ -418,10 +422,10 @@ def probplot(x, sparams=(), dist='norm', fit=True, plot=None):
     osr = sort(x)
     if fit or (plot is not None):
         # perform a linear fit.
-        slope, intercept, r, prob, sterrest = stats.linregress(osm, osr)
+        slope, intercept, r, _prob, _sterrest = stats.linregress(osm, osr)
 
     if plot is not None:
-        plot.plot(osm, osr, 'bo', osm, slope*osm + intercept, 'r-')
+        plot.plot(osm, osr, 'bo', osm, slope * osm + intercept, 'r-')
         try:
             if hasattr(plot, 'set_title'):
                 # Matplotlib Axes instance or something that looks like it
@@ -453,7 +457,7 @@ def probplot(x, sparams=(), dist='norm', fit=True, plot=None):
         return osm, osr
 
 
-def ppcc_max(x, brack=(0.0,1.0), dist='tukeylambda'):
+def ppcc_max(x, brack=(0.0, 1.0), dist='tukeylambda'):
     """Returns the shape parameter that maximizes the probability plot
     correlation coefficient for the given data to a one-parameter
     family of distributions.
@@ -470,13 +474,14 @@ def ppcc_max(x, brack=(0.0,1.0), dist='tukeylambda'):
     #  correlation
     def tempfunc(shape, mi, yvals, func):
         xvals = func(mi, shape)
-        r, prob = stats.pearsonr(xvals, yvals)
+        r, _prob = stats.pearsonr(xvals, yvals)
-        return 1-r
+        return 1 - r
 
-    return optimize.brent(tempfunc, brack=brack, args=(osm_uniform, osr, dist.ppf))
+    return optimize.brent(tempfunc, brack=brack,
+                          args=(osm_uniform, osr, dist.ppf))
 
 
-def ppcc_plot(x,a,b,dist='tukeylambda', plot=None, N=80):
+def ppcc_plot(x, a, b, dist='tukeylambda', plot=None, N=80):
     """Returns (shape, ppcc), and optionally plots shape vs. ppcc
     (probability plot correlation coefficient) as a function of shape
     parameter for a one-parameter family of distributions from shape
@@ -485,10 +490,10 @@ def ppcc_plot(x,a,b,dist='tukeylambda', plot=None, N=80):
     See also ppcc_max
     """
     svals = r_[a:b:complex(N)]
-    ppcc = svals*0.0
+    ppcc = svals * 0.0
     k = 0
     for sval in svals:
-        r1,r2 = probplot(x,sval,dist=dist,fit=1)
+        _r1, r2 = probplot(x, sval, dist=dist, fit=1)
         ppcc[k] = r2[-1]
         k += 1
     if plot is not None:
@@ -586,7 +591,7 @@ def boxcox_llf(lmb, data):
     y = boxcox(data, lmb)
     y_mean = np.mean(y, axis=0)
     llf = (lmb - 1) * np.sum(np.log(data), axis=0)
-    llf -= N / 2.0 * np.log(np.sum((y - y_mean)**2. / N, axis=0))
+    llf -= N / 2.0 * np.log(np.sum((y - y_mean) ** 2. / N, axis=0))
     return llf
 
 
@@ -719,7 +724,7 @@ def boxcox(x, lmbda=None, alpha=None):
         raise ValueError("Data must be positive.")
 
     if lmbda is not None:  # single transformation
-        return special.boxcox(x, lmbda)
+        return special.boxcox(x, lmbda)  # @UndefinedVariable
 
     # If lmbda=None, find the lmbda that maximizes the log-likelihood function.
     lmax = boxcox_normmax(x, method='mle')
@@ -810,7 +815,7 @@ def boxcox_normmax(x, brack=(-2.0, 2.0), method='pearsonr'):
             # correlation.
             y = boxcox(samps, lmbda)
             yvals = np.sort(y)
-            r, prob = stats.pearsonr(xvals, yvals)
+            r, _prob = stats.pearsonr(xvals, yvals)
             return 1 - r
 
         return optimize.brent(_eval_pearsonr, brack=brack, args=(xvals, x))
@@ -981,15 +986,15 @@ def shapiro(x, a=None, reta=False):
     if N < 3:
         raise ValueError("Data must be at least length 3.")
     if a is None:
-        a = zeros(N,'f')
+        a = zeros(N, 'f')
         init = 0
     else:
-        if len(a) != N//2:
+        if len(a) != N // 2:
             raise ValueError("len(a) must equal len(x)/2")
         init = 1
     y = sort(x)
-    a, w, pw, ifault = statlib.swilk(y, a[:N//2], init)
+    a, w, pw, ifault = statlib.swilk(y, a[:N // 2], init)
-    if not ifault in [0,2]:
+    if not ifault in [0, 2]:
         warnings.warn(str(ifault))
     if N > 5000:
         warnings.warn("p-value may not be accurate for N > 5000.")
@@ -1012,7 +1017,7 @@ _Avals_gumbel = array([0.474, 0.637, 0.757, 0.877, 1.038])
 _Avals_logistic = array([0.426, 0.563, 0.660, 0.769, 0.906, 1.010])
 
 
-def anderson(x,dist='norm'):
+def anderson(x, dist='norm'):
     """
     Anderson-Darling test for data coming from a particular distribution
 
@@ -1078,7 +1083,7 @@ def anderson(x,dist='norm'):
            pp. 591-595.
 
     """
-    if not dist in ['norm','expon','gumbel','extreme1','logistic']:
+    if not dist in ['norm', 'expon', 'gumbel', 'extreme1', 'logistic']:
         raise ValueError("Invalid distribution; dist must be 'norm', "
                          "'expon', 'gumbel', 'extreme1' or 'logistic'.")
     y = sort(x)
@@ -1086,52 +1091,52 @@ def anderson(x,dist='norm'):
     N = len(y)
     if dist == 'norm':
         s = np.std(x, ddof=1, axis=0)
-        w = (y-xbar)/s
+        w = (y - xbar) / s
         z = distributions.norm.cdf(w)
-        sig = array([15,10,5,2.5,1])
+        sig = array([15, 10, 5, 2.5, 1])
-        critical = around(_Avals_norm / (1.0 + 4.0/N - 25.0/N/N),3)
+        critical = around(_Avals_norm / (1.0 + 4.0 / N - 25.0 / N / N), 3)
     elif dist == 'expon':
         w = y / xbar
         z = distributions.expon.cdf(w)
-        sig = array([15,10,5,2.5,1])
+        sig = array([15, 10, 5, 2.5, 1])
-        critical = around(_Avals_expon / (1.0 + 0.6/N),3)
+        critical = around(_Avals_expon / (1.0 + 0.6 / N), 3)
     elif dist == 'logistic':
-        def rootfunc(ab,xj,N):
+        def rootfunc(ab, xj, N):
-            a,b = ab
+            a, b = ab
-            tmp = (xj-a)/b
+            tmp = (xj - a) / b
             tmp2 = exp(tmp)
-            val = [sum(1.0/(1+tmp2),axis=0)-0.5*N,
+            val = [sum(1.0 / (1 + tmp2), axis=0) - 0.5 * N,
-                   sum(tmp*(1.0-tmp2)/(1+tmp2),axis=0)+N]
+                   sum(tmp * (1.0 - tmp2) / (1 + tmp2), axis=0) + N]
             return array(val)
-        sol0 = array([xbar,np.std(x, ddof=1, axis=0)])
+        sol0 = array([xbar, np.std(x, ddof=1, axis=0)])
-        sol = optimize.fsolve(rootfunc,sol0,args=(x,N),xtol=1e-5)
+        sol = optimize.fsolve(rootfunc, sol0, args=(x, N), xtol=1e-5)
-        w = (y-sol[0])/sol[1]
+        w = (y - sol[0]) / sol[1]
         z = distributions.logistic.cdf(w)
-        sig = array([25,10,5,2.5,1,0.5])
+        sig = array([25, 10, 5, 2.5, 1, 0.5])
-        critical = around(_Avals_logistic / (1.0+0.25/N),3)
+        critical = around(_Avals_logistic / (1.0 + 0.25 / N), 3)
     else:  # (dist == 'gumbel') or (dist == 'extreme1'):
         # the following is incorrect, see ticket:1097
-##        def fixedsolve(th,xj,N):
+# def fixedsolve(th,xj,N):
 ##            val = stats.sum(xj)*1.0/N
 ##            tmp = exp(-xj/th)
 ##            term = sum(xj*tmp,axis=0)
 ##            term /= sum(tmp,axis=0)
-##            return val - term
+# return val - term
 ##        s = optimize.fixed_point(fixedsolve, 1.0, args=(x,N),xtol=1e-5)
 ##        xbar = -s*log(sum(exp(-x/s),axis=0)*1.0/N)
         xbar, s = distributions.gumbel_l.fit(x)
-        w = (y-xbar)/s
+        w = (y - xbar) / s
         z = distributions.gumbel_l.cdf(w)
-        sig = array([25,10,5,2.5,1])
+        sig = array([25, 10, 5, 2.5, 1])
-        critical = around(_Avals_gumbel / (1.0 + 0.2/sqrt(N)),3)
+        critical = around(_Avals_gumbel / (1.0 + 0.2 / sqrt(N)), 3)
 
-    i = arange(1,N+1)
+    i = arange(1, N + 1)
-    S = sum((2*i-1.0)/N*(log(z)+log(1-z[::-1])),axis=0)
+    S = sum((2 * i - 1.0) / N * (log(z) + log(1 - z[::-1])), axis=0)
-    A2 = -N-S
+    A2 = -N - S
     return A2, critical, sig
 
 
-def ansari(x,y):
+def ansari(x, y):
     """
     Perform the Ansari-Bradley test for equal scale parameters
 
@@ -1168,56 +1173,58 @@ def ansari(x,y):
            methods.  3rd ed. Chapman and Hall/CRC. 2001.  Section 5.8.2.
 
     """
-    x,y = asarray(x),asarray(y)
+    x, y = asarray(x), asarray(y)
     n = len(x)
     m = len(y)
     if m < 1:
         raise ValueError("Not enough other observations.")
     if n < 1:
         raise ValueError("Not enough test observations.")
-    N = m+n
+    N = m + n
-    xy = r_[x,y]  # combine
+    xy = r_[x, y]  # combine
     rank = stats.rankdata(xy)
-    symrank = amin(array((rank,N-rank+1)),0)
+    symrank = amin(array((rank, N - rank + 1)), 0)
-    AB = sum(symrank[:n],axis=0)
+    AB = sum(symrank[:n], axis=0)
     uxy = unique(xy)
     repeats = (len(uxy) != len(xy))
     exact = ((m < 55) and (n < 55) and not repeats)
     if repeats and ((m < 55) or (n < 55)):
         warnings.warn("Ties preclude use of exact statistic.")
     if exact:
-        astart, a1, ifault = statlib.gscale(n,m)
+        astart, a1, _ifault = statlib.gscale(n, m)
-        ind = AB-astart
+        ind = AB - astart
-        total = sum(a1,axis=0)
+        total = sum(a1, axis=0)
-        if ind < len(a1)/2.0:
+        if ind < len(a1) / 2.0:
             cind = int(ceil(ind))
             if (ind == cind):
-                pval = 2.0*sum(a1[:cind+1],axis=0)/total
+                pval = 2.0 * sum(a1[:cind + 1], axis=0) / total
             else:
-                pval = 2.0*sum(a1[:cind],axis=0)/total
+                pval = 2.0 * sum(a1[:cind], axis=0) / total
         else:
             find = int(floor(ind))
             if (ind == floor(ind)):
-                pval = 2.0*sum(a1[find:],axis=0)/total
+                pval = 2.0 * sum(a1[find:], axis=0) / total
             else:
-                pval = 2.0*sum(a1[find+1:],axis=0)/total
+                pval = 2.0 * sum(a1[find + 1:], axis=0) / total
-        return AB, min(1.0,pval)
+        return AB, min(1.0, pval)
 
     # otherwise compute normal approximation
     if N % 2:  # N odd
-        mnAB = n*(N+1.0)**2 / 4.0 / N
+        mnAB = n * (N + 1.0) ** 2 / 4.0 / N
-        varAB = n*m*(N+1.0)*(3+N**2)/(48.0*N**2)
+        varAB = n * m * (N + 1.0) * (3 + N ** 2) / (48.0 * N ** 2)
     else:
-        mnAB = n*(N+2.0)/4.0
+        mnAB = n * (N + 2.0) / 4.0
-        varAB = m*n*(N+2)*(N-2.0)/48/(N-1.0)
+        varAB = m * n * (N + 2) * (N - 2.0) / 48 / (N - 1.0)
     if repeats:   # adjust variance estimates
         # compute sum(tj * rj**2,axis=0)
-        fac = sum(symrank**2,axis=0)
+        fac = sum(symrank ** 2, axis=0)
         if N % 2:  # N odd
-            varAB = m*n*(16*N*fac-(N+1)**4)/(16.0 * N**2 * (N-1))
+            varAB = m * n * \
+                (16 * N * fac - (N + 1) ** 4) / (16.0 * N ** 2 * (N - 1))
         else:  # N even
-            varAB = m*n*(16*fac-N*(N+2)**2)/(16.0 * N * (N-1))
+            varAB = m * n * \
-    z = (AB - mnAB)/sqrt(varAB)
+                (16 * fac - N * (N + 2) ** 2) / (16.0 * N * (N - 1))
+    z = (AB - mnAB) / sqrt(varAB)
     pval = distributions.norm.sf(abs(z)) * 2.0
     return AB, pval
 
@@ -1255,20 +1262,21 @@ def bartlett(*args):
     if k < 2:
         raise ValueError("Must enter at least two input sample vectors.")
     Ni = zeros(k)
-    ssq = zeros(k,'d')
+    ssq = zeros(k, 'd')
     for j in range(k):
         Ni[j] = len(args[j])
         ssq[j] = np.var(args[j], ddof=1)
-    Ntot = sum(Ni,axis=0)
+    Ntot = sum(Ni, axis=0)
-    spsq = sum((Ni-1)*ssq,axis=0)/(1.0*(Ntot-k))
+    spsq = sum((Ni - 1) * ssq, axis=0) / (1.0 * (Ntot - k))
-    numer = (Ntot*1.0-k)*log(spsq) - sum((Ni-1.0)*log(ssq),axis=0)
+    numer = (Ntot * 1.0 - k) * log(spsq) - sum((Ni - 1.0) * log(ssq), axis=0)
-    denom = 1.0 + (1.0/(3*(k-1)))*((sum(1.0/(Ni-1.0),axis=0))-1.0/(Ntot-k))
+    denom = 1.0 + (1.0 / (3 * (k - 1))) * \
+        ((sum(1.0 / (Ni - 1.0), axis=0)) - 1.0 / (Ntot - k))
     T = numer / denom
-    pval = distributions.chi2.sf(T,k-1)  # 1 - cdf
+    pval = distributions.chi2.sf(T, k - 1)  # 1 - cdf
     return T, pval
 
 
-def levene(*args,**kwds):
+def levene(*args, **kwds):
     """
     Perform Levene test for equal variances.
 
@@ -1320,7 +1328,8 @@ def levene(*args,**kwds):
     proportiontocut = 0.05
     for kw, value in kwds.items():
         if kw not in ['center', 'proportiontocut']:
-            raise TypeError("levene() got an unexpected keyword argument '%s'" % kw)
+            raise TypeError(
+                "levene() got an unexpected keyword argument '%s'" % kw)
         if kw == 'center':
             center = value
         else:
@@ -1330,9 +1339,9 @@ def levene(*args,**kwds):
     if k < 2:
         raise ValueError("Must enter at least two input sample vectors.")
     Ni = zeros(k)
-    Yci = zeros(k,'d')
+    Yci = zeros(k, 'd')
 
-    if not center in ['mean','median','trimmed']:
+    if not center in ['mean', 'median', 'trimmed']:
         raise ValueError("Keyword argument <center> must be 'mean', 'median'"
                          + "or 'trimmed'.")
 
@@ -1341,42 +1350,43 @@ def levene(*args,**kwds):
     elif center == 'mean':
         func = lambda x: np.mean(x, axis=0)
     else:  # center == 'trimmed'
-        args = tuple(stats.trimboth(np.sort(arg), proportiontocut) for arg in args)
+        args = tuple(stats.trimboth(np.sort(arg), proportiontocut)
+                     for arg in args)
         func = lambda x: np.mean(x, axis=0)
 
     for j in range(k):
         Ni[j] = len(args[j])
         Yci[j] = func(args[j])
-    Ntot = sum(Ni,axis=0)
+    Ntot = sum(Ni, axis=0)
 
     # compute Zij's
-    Zij = [None]*k
+    Zij = [None] * k
     for i in range(k):
-        Zij[i] = abs(asarray(args[i])-Yci[i])
+        Zij[i] = abs(asarray(args[i]) - Yci[i])
     # compute Zbari
-    Zbari = zeros(k,'d')
+    Zbari = zeros(k, 'd')
     Zbar = 0.0
     for i in range(k):
         Zbari[i] = np.mean(Zij[i], axis=0)
-        Zbar += Zbari[i]*Ni[i]
+        Zbar += Zbari[i] * Ni[i]
     Zbar /= Ntot
 
-    numer = (Ntot-k)*sum(Ni*(Zbari-Zbar)**2,axis=0)
+    numer = (Ntot - k) * sum(Ni * (Zbari - Zbar) ** 2, axis=0)
 
     # compute denom_variance
     dvar = 0.0
     for i in range(k):
-        dvar += sum((Zij[i]-Zbari[i])**2,axis=0)
+        dvar += sum((Zij[i] - Zbari[i]) ** 2, axis=0)
 
-    denom = (k-1.0)*dvar
+    denom = (k - 1.0) * dvar
 
     W = numer / denom
-    pval = distributions.f.sf(W,k-1,Ntot-k)  # 1 - cdf
+    pval = distributions.f.sf(W, k - 1, Ntot - k)  # 1 - cdf
     return W, pval
 
 
 @setastest(False)
-def binom_test(x,n=None,p=0.5):
+def binom_test(x, n=None, p=0.5):
     """
     Perform a test that the probability of success is p.
 
@@ -1408,7 +1418,7 @@ def binom_test(x,n=None,p=0.5):
     """
     x = atleast_1d(x).astype(np.integer)
     if len(x) == 2:
-        n = x[1]+x[0]
+        n = x[1] + x[0]
         x = x[0]
     elif len(x) == 1:
         x = x[0]
@@ -1421,35 +1431,37 @@ def binom_test(x,n=None,p=0.5):
     if (p > 1.0) or (p < 0.0):
         raise ValueError("p must be in range [0,1]")
 
-    d = distributions.binom.pmf(x,n,p)
+    d = distributions.binom.pmf(x, n, p)
-    rerr = 1+1e-7
+    rerr = 1 + 1e-7
-    if (x == p*n):
+    if (x == p * n):
         # special case as shortcut, would also be handled by `else` below
         pval = 1.
-    elif (x < p*n):
+    elif (x < p * n):
-        i = np.arange(np.ceil(p*n),n+1)
+        i = np.arange(np.ceil(p * n), n + 1)
-        y = np.sum(distributions.binom.pmf(i,n,p) <= d*rerr,axis=0)
+        y = np.sum(distributions.binom.pmf(i, n, p) <= d * rerr, axis=0)
-        pval = distributions.binom.cdf(x,n,p) + distributions.binom.sf(n-y,n,p)
+        pval = distributions.binom.cdf(
+            x, n, p) + distributions.binom.sf(n - y, n, p)
     else:
-        i = np.arange(np.floor(p*n) + 1)
+        i = np.arange(np.floor(p * n) + 1)
-        y = np.sum(distributions.binom.pmf(i,n,p) <= d*rerr,axis=0)
+        y = np.sum(distributions.binom.pmf(i, n, p) <= d * rerr, axis=0)
-        pval = distributions.binom.cdf(y-1,n,p) + distributions.binom.sf(x-1,n,p)
+        pval = distributions.binom.cdf(
+            y - 1, n, p) + distributions.binom.sf(x - 1, n, p)
 
-    return min(1.0,pval)
+    return min(1.0, pval)
 
 
-def _apply_func(x,g,func):
+def _apply_func(x, g, func):
     # g is list of indices into x
     #  separating x into different groups
     #  func should be applied over the groups
-    g = unique(r_[0,g,len(x)])
+    g = unique(r_[0, g, len(x)])
     output = []
-    for k in range(len(g)-1):
+    for k in range(len(g) - 1):
-        output.append(func(x[g[k]:g[k+1]]))
+        output.append(func(x[g[k]:g[k + 1]]))
     return asarray(output)
 
 
-def fligner(*args,**kwds):
+def fligner(*args, **kwds):
     """
     Perform Fligner's test for equal variances.
 
@@ -1498,7 +1510,8 @@ def fligner(*args,**kwds):
     proportiontocut = 0.05
     for kw, value in kwds.items():
         if kw not in ['center', 'proportiontocut']:
-            raise TypeError("fligner() got an unexpected keyword argument '%s'" % kw)
+            raise TypeError(
+                "fligner() got an unexpected keyword argument '%s'" % kw)
         if kw == 'center':
             center = value
         else:
@@ -1508,7 +1521,7 @@ def fligner(*args,**kwds):
     if k < 2:
         raise ValueError("Must enter at least two input sample vectors.")
 
-    if not center in ['mean','median','trimmed']:
+    if not center in ['mean', 'median', 'trimmed']:
         raise ValueError("Keyword argument <center> must be 'mean', 'median'"
                          + "or 'trimmed'.")
 
@@ -1522,9 +1535,9 @@ def fligner(*args,**kwds):
 
     Ni = asarray([len(args[j]) for j in range(k)])
     Yci = asarray([func(args[j]) for j in range(k)])
-    Ntot = sum(Ni,axis=0)
+    Ntot = sum(Ni, axis=0)
     # compute Zij's
-    Zij = [abs(asarray(args[i])-Yci[i]) for i in range(k)]
+    Zij = [abs(asarray(args[i]) - Yci[i]) for i in range(k)]
     allZij = []
     g = [0]
     for i in range(k):
@@ -1532,14 +1545,14 @@ def fligner(*args,**kwds):
         g.append(len(allZij))
 
     ranks = stats.rankdata(allZij)
-    a = distributions.norm.ppf(ranks/(2*(Ntot+1.0)) + 0.5)
+    a = distributions.norm.ppf(ranks / (2 * (Ntot + 1.0)) + 0.5)
 
     # compute Aibar
-    Aibar = _apply_func(a,g,sum) / Ni
+    Aibar = _apply_func(a, g, sum) / Ni
     anbar = np.mean(a, axis=0)
-    varsq = np.var(a,axis=0, ddof=1)
+    varsq = np.var(a, axis=0, ddof=1)
-    Xsq = sum(Ni*(asarray(Aibar)-anbar)**2.0,axis=0)/varsq
+    Xsq = sum(Ni * (asarray(Aibar) - anbar) ** 2.0, axis=0) / varsq
-    pval = distributions.chi2.sf(Xsq,k-1)  # 1 - cdf
+    pval = distributions.chi2.sf(Xsq, k - 1)  # 1 - cdf
     return Xsq, pval
 
 
@@ -1618,7 +1631,8 @@ def mood(x, y, axis=0):
         axis = 0
 
     # Determine shape of the result arrays
-    res_shape = tuple([x.shape[ax] for ax in range(len(x.shape)) if ax != axis])
+    res_shape = tuple([x.shape[ax]
+                      for ax in range(len(x.shape)) if ax != axis])
     if not (res_shape == tuple([y.shape[ax] for ax in range(len(y.shape)) if
                                 ax != axis])):
         raise ValueError("Dimensions of x and y on all axes except `axis` "
@@ -1728,14 +1742,16 @@ def wilcoxon(x, y=None, zero_method="wilcox", correction=False):
         x, y = map(asarray, (x, y))
         if len(x) != len(y):
             raise ValueError('Unequal N in wilcoxon.  Aborting.')
-        d = x-y
+        d = x - y
 
     if zero_method == "wilcox":
-        d = compress(not_equal(d, 0), d, axis=-1)  # Keep all non-zero differences
+        # Keep all non-zero differences
+        d = compress(not_equal(d, 0), d, axis=-1)
 
     count = len(d)
     if (count < 10):
-        warnings.warn("Warning: sample size too small for normal approximation.")
+        warnings.warn(
+            "Warning: sample size too small for normal approximation.")
     r = stats.rankdata(abs(d))
     r_plus = sum((d > 0) * r, axis=0)
     r_minus = sum((d < 0) * r, axis=0)
@@ -1746,13 +1762,13 @@ def wilcoxon(x, y=None, zero_method="wilcox", correction=False):
         r_minus += r_zero / 2.
 
     T = min(r_plus, r_minus)
-    mn = count*(count + 1.) * 0.25
+    mn = count * (count + 1.) * 0.25
-    se = count*(count + 1.) * (2. * count + 1.)
+    se = count * (count + 1.) * (2. * count + 1.)
 
     if zero_method == "pratt":
         r = r[d != 0]
 
-    replist, repnum = find_repeats(r)
+    _replist, repnum = find_repeats(r)
     if repnum.size != 0:
         # Correction for repeated elements.
         se -= 0.5 * (repnum * (repnum * repnum - 1)).sum()
@@ -1770,35 +1786,35 @@ def _hermnorm(N):
     #  using the recursive relationship
     #  p_n+1 = p_n(x)' - x*p_n(x)
     #   and p_0(x) = 1
-    plist = [None]*N
+    plist = [None] * N
     plist[0] = poly1d(1)
-    for n in range(1,N):
+    for n in range(1, N):
-        plist[n] = plist[n-1].deriv() - poly1d([1,0])*plist[n-1]
+        plist[n] = plist[n - 1].deriv() - poly1d([1, 0]) * plist[n - 1]
     return plist
 
 
-def pdf_fromgamma(g1,g2,g3=0.0,g4=None):
+def pdf_fromgamma(g1, g2, g3=0.0, g4=None):
     if g4 is None:
-        g4 = 3*g2*g2
+        g4 = 3 * g2 * g2
-    sigsq = 1.0/g2
+    sigsq = 1.0 / g2
     sig = sqrt(sigsq)
-    mu = g1*sig**3.0
+    mu = g1 * sig ** 3.0
     p12 = _hermnorm(13)
     for k in range(13):
-        p12[k] = p12[k]/sig**k
+        p12[k] = p12[k] / sig ** k
 
     # Add all of the terms to polynomial
-    totp = p12[0] - (g1/6.0*p12[3]) + \
+    totp = p12[0] - (g1 / 6.0 * p12[3]) + \
-           (g2/24.0*p12[4] + g1*g1/72.0*p12[6]) - \
+        (g2 / 24.0 * p12[4] + g1 * g1 / 72.0 * p12[6]) - \
-           (g3/120.0*p12[5] + g1*g2/144.0*p12[7] + g1**3.0/1296.0*p12[9]) + \
+        (g3 / 120.0 * p12[5] + g1 * g2 / 144.0 * p12[7] + g1 ** 3.0 / 1296.0 * p12[9]) + \
-           (g4/720*p12[6] + (g2*g2/1152.0+g1*g3/720)*p12[8] +
+        (g4 / 720 * p12[6] + (g2 * g2 / 1152.0 + g1 * g3 / 720) * p12[8] +
-            g1*g1*g2/1728.0*p12[10] + g1**4.0/31104.0*p12[12])
+            g1 * g1 * g2 / 1728.0 * p12[10] + g1 ** 4.0 / 31104.0 * p12[12])
     # Final normalization
-    totp = totp / sqrt(2*pi)/sig
+    totp = totp / sqrt(2 * pi) / sig
 
     def thefunc(x):
-        xn = (x-mu)/sig
+        xn = (x - mu) / sig
-        return totp(xn)*exp(-xn*xn/2.0)
+        return totp(xn) * exp(-xn * xn / 2.0)
     return thefunc
 
 
@@ -1807,11 +1823,11 @@ def _circfuncs_common(samples, high, low):
     if samples.size == 0:
         return np.nan, np.nan
 
-    ang = (samples - low)*2*pi / (high-low)
+    ang = (samples - low) * 2 * pi / (high - low)
     return samples, ang
 
 
-def circmean(samples, high=2*pi, low=0, axis=None):
+def circmean(samples, high=2 * pi, low=0, axis=None):
     """
     Compute the circular mean for samples in a range.
 
@@ -1834,17 +1850,17 @@ def circmean(samples, high=2*pi, low=0, axis=None):
 
     """
     samples, ang = _circfuncs_common(samples, high, low)
-    res = angle(np.mean(exp(1j*ang), axis=axis))
+    res = angle(np.mean(exp(1j * ang), axis=axis))
     mask = res < 0
     if (mask.ndim > 0):
-        res[mask] += 2*pi
+        res[mask] += 2 * pi
     elif mask:
-        res = res + 2*pi
+        res = res + 2 * pi
 
-    return res*(high-low)/2.0/pi + low
+    return res * (high - low) / 2.0 / pi + low
 
 
-def circvar(samples, high=2*pi, low=0, axis=None):
+def circvar(samples, high=2 * pi, low=0, axis=None):
     """
     Compute the circular variance for samples assumed to be in a range
 
@@ -1872,12 +1888,12 @@ def circvar(samples, high=2*pi, low=0, axis=None):
 
     """
     samples, ang = _circfuncs_common(samples, high, low)
-    res = np.mean(exp(1j*ang), axis=axis)
+    res = np.mean(exp(1j * ang), axis=axis)
     R = abs(res)
-    return ((high-low)/2.0/pi)**2 * 2 * log(1/R)
+    return ((high - low) / 2.0 / pi) ** 2 * 2 * log(1 / R)
 
 
-def circstd(samples, high=2*pi, low=0, axis=None):
+def circstd(samples, high=2 * pi, low=0, axis=None):
     """
     Compute the circular standard deviation for samples assumed to be in the
     range [low to high].
@@ -1907,13 +1923,13 @@ def circstd(samples, high=2*pi, low=0, axis=None):
 
     """
     samples, ang = _circfuncs_common(samples, high, low)
-    res = np.mean(exp(1j*ang), axis=axis)
+    res = np.mean(exp(1j * ang), axis=axis)
     R = abs(res)
-    return ((high-low)/2.0/pi) * sqrt(-2*log(R))
+    return ((high - low) / 2.0 / pi) * sqrt(-2 * log(R))
 
 
 # Tests to include (from R) -- some of these already in stats.
-########
+#
 # X Ansari-Bradley
 # X Bartlett (and Levene)
 # X Binomial

pywafo/src/wafo/stats/stats.py

@@ -185,9 +185,9 @@ import numpy as np
 from . import futil
 from . import distributions
 try:
-    from ._rank import rankdata, tiecorrect
+    from scipy.stats._rank import rankdata, tiecorrect
 except:
-    rankdata=tiecorrect=None
+    rankdata = tiecorrect = None
 __all__ = ['find_repeats', 'gmean', 'hmean', 'mode',
            'tmean', 'tvar', 'tmin', 'tmax', 'tstd', 'tsem',
            'moment', 'variation', 'skew', 'kurtosis', 'describe',
@@ -257,12 +257,12 @@ def find_repeats(arr):
         (array([ 4., 5.]), array([2, 2], dtype=int32))
 
     """
-    v1,v2, n = futil.dfreps(arr)
+    v1, v2, n = futil.dfreps(arr)
-    return v1[:n],v2[:n]
+    return v1[:n], v2[:n]
 
-#######
+#
-### NAN friendly functions
+# NAN friendly functions
-########
+#
 
 
 def nanmean(x, axis=0):
@@ -387,8 +387,8 @@ def _nanmedian(arr1d):  # This only works on 1d arrays
     m : float
         The median.
     """
-    cond = 1-np.isnan(arr1d)
+    cond = 1 - np.isnan(arr1d)
-    x = np.sort(np.compress(cond,arr1d,axis=-1))
+    x = np.sort(np.compress(cond, arr1d, axis=-1))
     if x.size == 0:
         return np.nan
     return np.median(x)
@@ -451,9 +451,9 @@ def nanmedian(x, axis=0):
     return x
 
 
-#####################################
+#
-########  CENTRAL TENDENCY  ########
+# CENTRAL TENDENCY  ########
-#####################################
+#
 
 
 def gmean(a, axis=0, dtype=None):
@@ -498,10 +498,11 @@ def gmean(a, axis=0, dtype=None):
     arrays automatically mask any non-finite values.
 
     """
-    if not isinstance(a, np.ndarray):  # if not an ndarray object attempt to convert it
+    # if not an ndarray object attempt to convert it
+    if not isinstance(a, np.ndarray):
         log_a = np.log(np.array(a, dtype=dtype))
     elif dtype:  # Must change the default dtype allowing array type
-        if isinstance(a,np.ma.MaskedArray):
+        if isinstance(a, np.ma.MaskedArray):
             log_a = np.log(np.ma.asarray(a, dtype=dtype))
         else:
             log_a = np.log(np.asarray(a, dtype=dtype))
@@ -561,9 +562,10 @@ def hmean(a, axis=0, dtype=None):
                 size = a.shape[0]
             else:
                 size = a.shape[axis]
-        return size / np.sum(1.0/a, axis=axis, dtype=dtype)
+        return size / np.sum(1.0 / a, axis=axis, dtype=dtype)
     else:
-        raise ValueError("Harmonic mean only defined if all elements greater than zero")
+        raise ValueError(
+            "Harmonic mean only defined if all elements greater than zero")
 
 
 def mode(a, axis=0):
@@ -612,7 +614,7 @@ def mode(a, axis=0):
     oldcounts = np.zeros(testshape)
     for score in scores:
         template = (a == score)
-        counts = np.expand_dims(np.sum(template, axis),axis)
+        counts = np.expand_dims(np.sum(template, axis), axis)
         mostfrequent = np.where(counts > oldcounts, score, oldmostfreq)
         oldcounts = np.maximum(counts, oldcounts)
         oldmostfreq = mostfrequent
@@ -700,7 +702,7 @@ def tmean(a, limits=None, inclusive=(True, True)):
 
 def masked_var(am):
     m = am.mean()
-    s = ma.add.reduce((am - m)**2)
+    s = ma.add.reduce((am - m) ** 2)
     n = am.count() - 1.0
     return s / n
 
@@ -741,7 +743,7 @@ def tvar(a, limits=None, inclusive=(True, True)):
     a = a.astype(float).ravel()
     if limits is None:
         n = len(a)
-        return a.var()*(n/(n-1.))
+        return a.var() * (n / (n - 1.))
     am = mask_to_limits(a, limits, inclusive)
     return masked_var(am)
 
@@ -885,9 +887,9 @@ def tsem(a, limits=None, inclusive=(True, True)):
     return sd / np.sqrt(am.count())
 
 
-#####################################
+#
-############  MOMENTS  #############
+# MOMENTS  #############
-#####################################
+#
 
 def moment(a, moment=1, axis=0):
     """
@@ -926,8 +928,8 @@ def moment(a, moment=1, axis=0):
             # the input was 1D, so return a scalar instead of a rank-0 array
             return np.float64(0.0)
     else:
-        mn = np.expand_dims(np.mean(a,axis), axis)
+        mn = np.expand_dims(np.mean(a, axis), axis)
-        s = np.power((a-mn), moment)
+        s = np.power((a - mn), moment)
         return np.mean(s, axis)
 
 
@@ -951,7 +953,7 @@ def variation(a, axis=0):
 
     """
     a, axis = _chk_asarray(a, axis)
-    return a.std(axis)/a.mean(axis)
+    return a.std(axis) / a.mean(axis)
 
 
 def skew(a, axis=0, bias=True):
@@ -987,18 +989,18 @@ def skew(a, axis=0, bias=True):
        York. 2000.
 
     """
-    a, axis = _chk_asarray(a,axis)
+    a, axis = _chk_asarray(a, axis)
     n = a.shape[axis]
     m2 = moment(a, 2, axis)
     m3 = moment(a, 3, axis)
     zero = (m2 == 0)
-    vals = np.where(zero, 0, m3 / m2**1.5)
+    vals = np.where(zero, 0, m3 / m2 ** 1.5)
     if not bias:
         can_correct = (n > 2) & (m2 > 0)
         if can_correct.any():
             m2 = np.extract(can_correct, m2)
             m3 = np.extract(can_correct, m3)
-            nval = np.sqrt((n-1.0)*n)/(n-2.0)*m3/m2**1.5
+            nval = np.sqrt((n - 1.0) * n) / (n - 2.0) * m3 / m2 ** 1.5
             np.place(vals, can_correct, nval)
     if vals.ndim == 0:
         return vals.item()
@@ -1045,12 +1047,12 @@ def kurtosis(a, axis=0, fisher=True, bias=True):
     """
     a, axis = _chk_asarray(a, axis)
     n = a.shape[axis]
-    m2 = moment(a,2,axis)
+    m2 = moment(a, 2, axis)
-    m4 = moment(a,4,axis)
+    m4 = moment(a, 4, axis)
     zero = (m2 == 0)
     olderr = np.seterr(all='ignore')
     try:
-        vals = np.where(zero, 0, m4 / m2**2.0)
+        vals = np.where(zero, 0, m4 / m2 ** 2.0)
     finally:
         np.seterr(**olderr)
 
@@ -1059,8 +1061,10 @@ def kurtosis(a, axis=0, fisher=True, bias=True):
         if can_correct.any():
             m2 = np.extract(can_correct, m2)
             m4 = np.extract(can_correct, m4)
-            nval = 1.0/(n-2)/(n-3)*((n*n-1.0)*m4/m2**2.0-3*(n-1)**2.0)
+            nval = 1.0 / \
-            np.place(vals, can_correct, nval+3.0)
+                (n - 2) / (n - 3) * \
+                ((n * n - 1.0) * m4 / m2 ** 2.0 - 3 * (n - 1) ** 2.0)
+            np.place(vals, can_correct, nval + 3.0)
 
     if vals.ndim == 0:
         vals = vals.item()  # array scalar
@@ -1116,9 +1120,9 @@ def describe(a, axis=0):
     kurt = kurtosis(a, axis)
     return n, mm, m, v, sk, kurt
 
-#####################################
+#
-########  NORMALITY TESTS  ##########
+# NORMALITY TESTS  ##########
-#####################################
+#
 
 
 def skewtest(a, axis=0):
@@ -1206,17 +1210,20 @@ def kurtosistest(a, axis=0):
             "kurtosistest only valid for n>=20 ... continuing anyway, n=%i" %
             int(n))
     b2 = kurtosis(a, axis, fisher=False)
-    E = 3.0*(n-1) / (n+1)
+    E = 3.0 * (n - 1) / (n + 1)
-    varb2 = 24.0*n*(n-2)*(n-3) / ((n+1)*(n+1)*(n+3)*(n+5))
+    varb2 = 24.0 * n * \
-    x = (b2-E)/np.sqrt(varb2)
+        (n - 2) * (n - 3) / ((n + 1) * (n + 1) * (n + 3) * (n + 5))
-    sqrtbeta1 = 6.0*(n*n-5*n+2)/((n+7)*(n+9)) * np.sqrt((6.0*(n+3)*(n+5)) /
+    x = (b2 - E) / np.sqrt(varb2)
-                                                       (n*(n-2)*(n-3)))
+    sqrtbeta1 = 6.0 * (n * n - 5 * n + 2) / ((n + 7) * (n + 9)) * np.sqrt((6.0 * (n + 3) * (n + 5)) /
-    A = 6.0 + 8.0/sqrtbeta1 * (2.0/sqrtbeta1 + np.sqrt(1+4.0/(sqrtbeta1**2)))
+                                                                         (n * (n - 2) * (n - 3)))
-    term1 = 1 - 2/(9.0*A)
+    A = 6.0 + 8.0 / sqrtbeta1 * \
-    denom = 1 + x*np.sqrt(2/(A-4.0))
+        (2.0 / sqrtbeta1 + np.sqrt(1 + 4.0 / (sqrtbeta1 ** 2)))
+    term1 = 1 - 2 / (9.0 * A)
+    denom = 1 + x * np.sqrt(2 / (A - 4.0))
     denom = np.where(denom < 0, 99, denom)
-    term2 = np.where(denom < 0, term1, np.power((1-2.0/A)/denom,1/3.0))
+    term2 = np.where(
-    Z = (term1 - term2) / np.sqrt(2/(9.0*A))
+        denom < 0, term1, np.power((1 - 2.0 / A) / denom, 1 / 3.0))
+    Z = (term1 - term2) / np.sqrt(2 / (9.0 * A))
     Z = np.where(denom == 99, 0, Z)
     if Z.ndim == 0:
         Z = Z[()]
@@ -1261,10 +1268,10 @@ def normaltest(a, axis=0):
 
     """
     a, axis = _chk_asarray(a, axis)
-    s,p = skewtest(a,axis)
+    s, _p = skewtest(a, axis)
-    k,p = kurtosistest(a,axis)
+    k, _p = kurtosistest(a, axis)
-    k2 = s*s + k*k
+    k2 = s * s + k * k
-    return k2, chisqprob(k2,2)
+    return k2, chisqprob(k2, 2)
 
 
 def jarque_bera(x):
@@ -1315,17 +1322,18 @@ def jarque_bera(x):
 
     mu = x.mean()
     diffx = x - mu
-    skewness = (1 / n * np.sum(diffx**3)) / (1 / n * np.sum(diffx**2))**(3 / 2.)
+    skewness = (1 / n * np.sum(diffx ** 3)) / \
-    kurtosis = (1 / n * np.sum(diffx**4)) / (1 / n * np.sum(diffx**2))**2
+        (1 / n * np.sum(diffx ** 2)) ** (3 / 2.)
-    jb_value = n / 6 * (skewness**2 + (kurtosis - 3)**2 / 4)
+    kurtosis = (1 / n * np.sum(diffx ** 4)) / (1 / n * np.sum(diffx ** 2)) ** 2
+    jb_value = n / 6 * (skewness ** 2 + (kurtosis - 3) ** 2 / 4)
     p = 1 - distributions.chi2.cdf(jb_value, 2)
 
     return jb_value, p
 
 
-#####################################
+#
-######  FREQUENCY FUNCTIONS  #######
+# FREQUENCY FUNCTIONS  #######
-#####################################
+#
 
 def itemfreq(a):
     """
@@ -1597,10 +1605,11 @@ def histogram2(a, bins):
     # comment: probably obsoleted by numpy.histogram()
     n = np.searchsorted(np.sort(a), bins)
     n = np.concatenate([n, [len(a)]])
-    return n[1:]-n[:-1]
+    return n[1:] - n[:-1]
 
 
-def histogram(a, numbins=10, defaultlimits=None, weights=None, printextras=False):
+def histogram(a, numbins=10, defaultlimits=None, weights=None,
+              printextras=False):
     """
     Separates the range into several bins and returns the number of instances
     in each bin.
@@ -1714,9 +1723,9 @@ def cumfreq(a, numbins=10, defaultreallimits=None, weights=None):
     3
 
     """
-    h,l,b,e = histogram(a, numbins, defaultreallimits, weights=weights)
+    h, l, b, e = histogram(a, numbins, defaultreallimits, weights=weights)
-    cumhist = np.cumsum(h*1, axis=0)
+    cumhist = np.cumsum(h * 1, axis=0)
-    return cumhist,l,b,e
+    return cumhist, l, b, e
 
 
 def relfreq(a, numbins=10, defaultreallimits=None, weights=None):
@@ -1764,9 +1773,9 @@ def relfreq(a, numbins=10, defaultreallimits=None, weights=None):
     return h, l, b, e
 
 
-#####################################
+#
-######  VARIABILITY FUNCTIONS  #####
+# VARIABILITY FUNCTIONS  #####
-#####################################
+#
 
 def obrientransform(*args):
     """
@@ -1827,7 +1836,7 @@ def obrientransform(*args):
         a = np.asarray(arg)
         n = len(a)
         mu = np.mean(a)
-        sq = (a - mu)**2
+        sq = (a - mu) ** 2
         sumsq = sq.sum()
 
         # The O'Brien transform.
@@ -1879,7 +1888,7 @@ def signaltonoise(a, axis=0, ddof=0):
     a = np.asanyarray(a)
     m = a.mean(axis)
     sd = a.std(axis=axis, ddof=ddof)
-    return np.where(sd == 0, 0, m/sd)
+    return np.where(sd == 0, 0, m / sd)
 
 
 def sem(a, axis=0, ddof=1):
@@ -1927,7 +1936,7 @@ def sem(a, axis=0, ddof=1):
     """
     a, axis = _chk_asarray(a, axis)
     n = a.shape[axis]
-    s = np.std(a,axis=axis, ddof=ddof) / np.sqrt(n)  # JP check normalization
+    s = np.std(a, axis=axis, ddof=ddof) / np.sqrt(n)  # JP check normalization
     return s
 
 
@@ -1988,7 +1997,7 @@ def zscore(a, axis=0, ddof=0):
     sstd = a.std(axis=axis, ddof=ddof)
     if axis and mns.ndim < a.ndim:
         return ((a - np.expand_dims(mns, axis=axis)) /
-                 np.expand_dims(sstd,axis=axis))
+                np.expand_dims(sstd, axis=axis))
     else:
         return (a - mns) / sstd
 
@@ -2039,14 +2048,14 @@ def zmap(scores, compare, axis=0, ddof=0):
     sstd = compare.std(axis=axis, ddof=ddof)
     if axis and mns.ndim < compare.ndim:
         return ((scores - np.expand_dims(mns, axis=axis)) /
-                 np.expand_dims(sstd,axis=axis))
+                np.expand_dims(sstd, axis=axis))
     else:
         return (scores - mns) / sstd
 
 
-#####################################
+#
-#######  TRIMMING FUNCTIONS  #######
+# TRIMMING FUNCTIONS  #######
-#####################################
+#
 
 def threshold(a, threshmin=None, threshmax=None, newval=0):
     """
@@ -2150,10 +2159,10 @@ def sigmaclip(a, low=4., high=4.):
         c_std = c.std()
         c_mean = c.mean()
         size = c.size
-        critlower = c_mean - c_std*low
+        critlower = c_mean - c_std * low
-        critupper = c_mean + c_std*high
+        critupper = c_mean + c_std * high
         c = c[(c > critlower) & (c < critupper)]
-        delta = size-c.size
+        delta = size - c.size
     return c, critlower, critupper
 
 
@@ -2240,9 +2249,9 @@ def trim1(a, proportiontocut, tail='right'):
     a = asarray(a)
     if tail.lower() == 'right':
         lowercut = 0
-        uppercut = len(a) - int(proportiontocut*len(a))
+        uppercut = len(a) - int(proportiontocut * len(a))
     elif tail.lower() == 'left':
-        lowercut = int(proportiontocut*len(a))
+        lowercut = int(proportiontocut * len(a))
         uppercut = len(a)
 
     return a[lowercut:uppercut]
@@ -2421,7 +2430,7 @@ def pearsonr(x, y):
     n = len(x)
     mx = x.mean()
     my = y.mean()
-    xm, ym = x-mx, y-my
+    xm, ym = x - mx, y - my
     r_num = np.add.reduce(xm * ym)
     r_den = np.sqrt(ss(xm) * ss(ym))
     r = r_num / r_den
@@ -2429,12 +2438,12 @@ def pearsonr(x, y):
     # Presumably, if abs(r) > 1, then it is only some small artifact of floating
     # point arithmetic.
     r = max(min(r, 1.0), -1.0)
-    df = n-2
+    df = n - 2
     if abs(r) == 1.0:
         prob = 0.0
     else:
-        t_squared = r*r * (df / ((1.0 - r) * (1.0 + r)))
+        t_squared = r * r * (df / ((1.0 - r) * (1.0 + r)))
-        prob = betai(0.5*df, 0.5, df / (df + t_squared))
+        prob = betai(0.5 * df, 0.5, df / (df + t_squared))
     return r, prob
 
 
@@ -2497,7 +2506,8 @@ def fisher_exact(table, alternative='two-sided'):
 
     """
     hypergeom = distributions.hypergeom
-    c = np.asarray(table, dtype=np.int64)  # int32 is not enough for the algorithm
+    # int32 is not enough for the algorithm
+    c = np.asarray(table, dtype=np.int64)
     if not c.shape == (2, 2):
         raise ValueError("The input `table` must be of shape (2, 2).")
 
@@ -2509,14 +2519,14 @@ def fisher_exact(table, alternative='two-sided'):
         # the odds ratio is NaN.
         return np.nan, 1.0
 
-    if c[1,0] > 0 and c[0,1] > 0:
+    if c[1, 0] > 0 and c[0, 1] > 0:
-        oddsratio = c[0,0] * c[1,1] / float(c[1,0] * c[0,1])
+        oddsratio = c[0, 0] * c[1, 1] / float(c[1, 0] * c[0, 1])
     else:
         oddsratio = np.inf
 
-    n1 = c[0,0] + c[0,1]
+    n1 = c[0, 0] + c[0, 1]
-    n2 = c[1,0] + c[1,1]
+    n2 = c[1, 0] + c[1, 1]
-    n = c[0,0] + c[1,0]
+    n = c[0, 0] + c[1, 0]
 
     def binary_search(n, n1, n2, side):
         """Binary search for where to begin lower/upper halves in two-sided
@@ -2560,28 +2570,28 @@ def fisher_exact(table, alternative='two-sided'):
         return guess
 
     if alternative == 'less':
-        pvalue = hypergeom.cdf(c[0,0], n1 + n2, n1, n)
+        pvalue = hypergeom.cdf(c[0, 0], n1 + n2, n1, n)
     elif alternative == 'greater':
         # Same formula as the 'less' case, but with the second column.
-        pvalue = hypergeom.cdf(c[0,1], n1 + n2, n1, c[0,1] + c[1,1])
+        pvalue = hypergeom.cdf(c[0, 1], n1 + n2, n1, c[0, 1] + c[1, 1])
     elif alternative == 'two-sided':
         mode = int(float((n + 1) * (n1 + 1)) / (n1 + n2 + 2))
-        pexact = hypergeom.pmf(c[0,0], n1 + n2, n1, n)
+        pexact = hypergeom.pmf(c[0, 0], n1 + n2, n1, n)
         pmode = hypergeom.pmf(mode, n1 + n2, n1, n)
 
         epsilon = 1 - 1e-4
         if float(np.abs(pexact - pmode)) / np.abs(np.max(pexact, pmode)) <= 1 - epsilon:
             return oddsratio, 1.
 
-        elif c[0,0] < mode:
+        elif c[0, 0] < mode:
-            plower = hypergeom.cdf(c[0,0], n1 + n2, n1, n)
+            plower = hypergeom.cdf(c[0, 0], n1 + n2, n1, n)
             if hypergeom.pmf(n, n1 + n2, n1, n) > pexact / epsilon:
                 return oddsratio, plower
 
             guess = binary_search(n, n1, n2, "upper")
             pvalue = plower + hypergeom.sf(guess - 1, n1 + n2, n1, n)
         else:
-            pupper = hypergeom.sf(c[0,0] - 1, n1 + n2, n1, n)
+            pupper = hypergeom.sf(c[0, 0] - 1, n1 + n2, n1, n)
             if hypergeom.pmf(0, n1 + n2, n1, n) > pexact / epsilon:
                 return oddsratio, pupper
 
@@ -2691,24 +2701,24 @@ def spearmanr(a, b=None, axis=0):
 
     """
     a, axisout = _chk_asarray(a, axis)
-    ar = np.apply_along_axis(rankdata,axisout,a)
+    ar = np.apply_along_axis(rankdata, axisout, a)
 
     br = None
     if not b is None:
         b, axisout = _chk_asarray(b, axis)
-        br = np.apply_along_axis(rankdata,axisout,b)
+        br = np.apply_along_axis(rankdata, axisout, b)
     n = a.shape[axisout]
-    rs = np.corrcoef(ar,br,rowvar=axisout)
+    rs = np.corrcoef(ar, br, rowvar=axisout)
 
     olderr = np.seterr(divide='ignore')  # rs can have elements equal to 1
     try:
-        t = rs * np.sqrt((n-2) / ((rs+1.0)*(1.0-rs)))
+        t = rs * np.sqrt((n - 2) / ((rs + 1.0) * (1.0 - rs)))
     finally:
         np.seterr(**olderr)
-    prob = distributions.t.sf(np.abs(t),n-2)*2
+    prob = distributions.t.sf(np.abs(t), n - 2) * 2
 
-    if rs.shape == (2,2):
+    if rs.shape == (2, 2):
-        return rs[1,0], prob[1,0]
+        return rs[1, 0], prob[1, 0]
     else:
         return rs, prob
 
@@ -2770,13 +2780,13 @@ def pointbiserialr(x, y):
     y1m = y1.mean()
 
     # phat - phat**2 is more stable than phat*(1-phat)
-    rpb = (y1m - y0m) * np.sqrt(phat - phat**2) / y.std()
+    rpb = (y1m - y0m) * np.sqrt(phat - phat ** 2) / y.std()
 
-    df = n-2
+    df = n - 2
     # fixme: see comment about TINY in pearsonr()
     TINY = 1e-20
-    t = rpb*np.sqrt(df/((1.0-rpb+TINY)*(1.0+rpb+TINY)))
+    t = rpb * np.sqrt(df / ((1.0 - rpb + TINY) * (1.0 + rpb + TINY)))
-    prob = betai(0.5*df, 0.5, df/(df+t*t))
+    prob = betai(0.5 * df, 0.5, df / (df + t * t))
     return rpb, prob
 
 
@@ -2851,11 +2861,11 @@ def kendalltau(x, y, initial_lexsort=True):
         if length == 1:
             return 0
         if length == 2:
-            if y[perm[offs]] <= y[perm[offs+1]]:
+            if y[perm[offs]] <= y[perm[offs + 1]]:
                 return 0
             t = perm[offs]
-            perm[offs] = perm[offs+1]
+            perm[offs] = perm[offs + 1]
-            perm[offs+1] = t
+            perm[offs + 1] = t
             return 1
         length0 = length // 2
         length1 = length - length0
@@ -2879,7 +2889,7 @@ def kendalltau(x, y, initial_lexsort=True):
             if d > 0:
                 exchcnt += d
             i += 1
-        perm[offs:offs+length] = temp[0:length]
+        perm[offs:offs + length] = temp[0:length]
         return exchcnt
 
     # initial sort on values of x and, if tied, on values of y
@@ -2903,7 +2913,7 @@ def kendalltau(x, y, initial_lexsort=True):
     # compute ties in x
     first = 0
     u = 0
-    for i in xrange(1,n):
+    for i in xrange(1, n):
         if x[perm[first]] != x[perm[i]]:
             u += ((i - first) * (i - first - 1)) // 2
             first = i
@@ -2914,7 +2924,7 @@ def kendalltau(x, y, initial_lexsort=True):
     # compute ties in y after mergesort with counting
     first = 0
     v = 0
-    for i in xrange(1,n):
+    for i in xrange(1, n):
         if y[perm[first]] != y[perm[i]]:
             v += ((i - first) * (i - first - 1)) // 2
             first = i
@@ -2995,13 +3005,13 @@ def linregress(x, y=None):
         x = asarray(x)
         y = asarray(y)
     n = len(x)
-    xmean = np.mean(x,None)
+    xmean = np.mean(x, None)
-    ymean = np.mean(y,None)
+    ymean = np.mean(y, None)
 
     # average sum of squares:
     ssxm, ssxym, ssyxm, ssym = np.cov(x, y, bias=1).flat
     r_num = ssxym
-    r_den = np.sqrt(ssxm*ssym)
+    r_den = np.sqrt(ssxm * ssym)
     if r_den == 0.0:
         r = 0.0
     else:
@@ -3012,18 +3022,18 @@ def linregress(x, y=None):
         elif (r < -1.0):
             r = -1.0
 
-    df = n-2
+    df = n - 2
-    t = r*np.sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)))
+    t = r * np.sqrt(df / ((1.0 - r + TINY) * (1.0 + r + TINY)))
-    prob = distributions.t.sf(np.abs(t),df)*2
+    prob = distributions.t.sf(np.abs(t), df) * 2
     slope = r_num / ssxm
-    intercept = ymean - slope*xmean
+    intercept = ymean - slope * xmean
-    sterrest = np.sqrt((1-r*r)*ssym / ssxm / df)
+    sterrest = np.sqrt((1 - r * r) * ssym / ssxm / df)
     return slope, intercept, r, prob, sterrest
 
 
-#####################################
+#
-#####  INFERENTIAL STATISTICS  #####
+# INFERENTIAL STATISTICS  #####
-#####################################
+#
 
 def ttest_1samp(a, popmean, axis=0):
     """
@@ -3093,9 +3103,10 @@ def ttest_1samp(a, popmean, axis=0):
     return t, prob
 
 
-def _ttest_finish(df,t):
+def _ttest_finish(df, t):
     """Common code between all 3 t-test functions."""
-    prob = distributions.t.sf(np.abs(t), df) * 2  # use np.abs to get upper tail
+    # use np.abs to get upper tail
+    prob = distributions.t.sf(np.abs(t), df) * 2
     if t.ndim == 0:
         t = t[()]
 
@@ -3206,7 +3217,8 @@ def ttest_ind(a, b, axis=0, equal_var=True):
     else:
         vn1 = v1 / n1
         vn2 = v2 / n2
-        df = ((vn1 + vn2)**2) / ((vn1**2) / (n1 - 1) + (vn2**2) / (n2 - 1))
+        df = ((vn1 + vn2) ** 2) / \
+            ((vn1 ** 2) / (n1 - 1) + (vn2 ** 2) / (n2 - 1))
 
         # If df is undefined, variances are zero (assumes n1 > 0 & n2 > 0).
         # Hence it doesn't matter what df is as long as it's not NaN.
@@ -3411,8 +3423,8 @@ def kstest(rvs, cdf, args=(), N=20, alternative='two-sided', mode='approx'):
     if isinstance(cdf, string_types):
         cdf = getattr(distributions, cdf).cdf
     if callable(rvs):
-        kwds = {'size':N}
+        kwds = {'size': N}
-        vals = np.sort(rvs(*args,**kwds))
+        vals = np.sort(rvs(*args, **kwds))
     else:
         vals = np.sort(rvs)
         N = len(vals)
@@ -3423,25 +3435,25 @@ def kstest(rvs, cdf, args=(), N=20, alternative='two-sided', mode='approx'):
         alternative = 'two-sided'
 
     if alternative in ['two-sided', 'greater']:
-        Dplus = (np.arange(1.0, N+1)/N - cdfvals).max()
+        Dplus = (np.arange(1.0, N + 1) / N - cdfvals).max()
         if alternative == 'greater':
-            return Dplus, distributions.ksone.sf(Dplus,N)
+            return Dplus, distributions.ksone.sf(Dplus, N)
 
     if alternative in ['two-sided', 'less']:
-        Dmin = (cdfvals - np.arange(0.0, N)/N).max()
+        Dmin = (cdfvals - np.arange(0.0, N) / N).max()
         if alternative == 'less':
-            return Dmin, distributions.ksone.sf(Dmin,N)
+            return Dmin, distributions.ksone.sf(Dmin, N)
 
     if alternative == 'two-sided':
-        D = np.max([Dplus,Dmin])
+        D = np.max([Dplus, Dmin])
         if mode == 'asymp':
-            return D, distributions.kstwobign.sf(D*np.sqrt(N))
+            return D, distributions.kstwobign.sf(D * np.sqrt(N))
         if mode == 'approx':
-            pval_two = distributions.kstwobign.sf(D*np.sqrt(N))
+            pval_two = distributions.kstwobign.sf(D * np.sqrt(N))
-            if N > 2666 or pval_two > 0.80 - N*0.3/1000.0:
+            if N > 2666 or pval_two > 0.80 - N * 0.3 / 1000.0:
-                return D, distributions.kstwobign.sf(D*np.sqrt(N))
+                return D, distributions.kstwobign.sf(D * np.sqrt(N))
             else:
-                return D, distributions.ksone.sf(D,N)*2
+                return D, distributions.ksone.sf(D, N) * 2
 
 
 # Map from names to lambda_ values used in power_divergence().
@@ -3451,7 +3463,7 @@ _power_div_lambda_names = {
     "freeman-tukey": -0.5,
     "mod-log-likelihood": -1,
     "neyman": -2,
-    "cressie-read": 2/3,
+    "cressie-read": 2 / 3,
 }
 
 
@@ -3662,7 +3674,7 @@ def power_divergence(f_obs, f_exp=None, ddof=0, axis=0, lambda_=None):
     # cases of lambda_.
     if lambda_ == 1:
         # Pearson's chi-squared statistic
-        terms = (f_obs - f_exp)**2 / f_exp
+        terms = (f_obs - f_exp) ** 2 / f_exp
     elif lambda_ == 0:
         # Log-likelihood ratio (i.e. G-test)
         terms = 2.0 * special.xlogy(f_obs, f_obs / f_exp)
@@ -3671,7 +3683,7 @@ def power_divergence(f_obs, f_exp=None, ddof=0, axis=0, lambda_=None):
         terms = 2.0 * special.xlogy(f_exp, f_exp / f_obs)
     else:
         # General Cressie-Read power divergence.
-        terms = f_obs * ((f_obs / f_exp)**lambda_ - 1)
+        terms = f_obs * ((f_obs / f_exp) ** lambda_ - 1)
         terms /= 0.5 * lambda_ * (lambda_ + 1)
 
     stat = terms.sum(axis=axis)
@@ -3861,20 +3873,20 @@ def ks_2samp(data1, data2):
 
     """
     data1, data2 = map(asarray, (data1, data2))
-    n1 = data1.shape[0]
+    #n1 = data1.shape[0]
-    n2 = data2.shape[0]
+    #n2 = data2.shape[0]
     n1 = len(data1)
     n2 = len(data2)
     data1 = np.sort(data1)
     data2 = np.sort(data2)
-    data_all = np.concatenate([data1,data2])
+    data_all = np.concatenate([data1, data2])
-    cdf1 = np.searchsorted(data1,data_all,side='right')/(1.0*n1)
+    cdf1 = np.searchsorted(data1, data_all, side='right') / (1.0 * n1)
-    cdf2 = (np.searchsorted(data2,data_all,side='right'))/(1.0*n2)
+    cdf2 = (np.searchsorted(data2, data_all, side='right')) / (1.0 * n2)
-    d = np.max(np.absolute(cdf1-cdf2))
+    d = np.max(np.absolute(cdf1 - cdf2))
     # Note: d absolute not signed distance
-    en = np.sqrt(n1*n2/float(n1+n2))
+    en = np.sqrt(n1 * n2 / float(n1 + n2))
     try:
-        prob = ksprob((en+0.12+0.11/en)*d)
+        prob = ksprob((en + 0.12 + 0.11 / en) * d)
     except:
         prob = 1.0
     return d, prob
@@ -3915,22 +3927,24 @@ def mannwhitneyu(x, y, use_continuity=True):
     y = asarray(y)
     n1 = len(x)
     n2 = len(y)
-    ranked = rankdata(np.concatenate((x,y)))
+    ranked = rankdata(np.concatenate((x, y)))
     rankx = ranked[0:n1]       # get the x-ranks
-    u1 = n1*n2 + (n1*(n1+1))/2.0 - np.sum(rankx,axis=0)  # calc U for x
+    # calc U for x
-    u2 = n1*n2 - u1                            # remainder is U for y
+    u1 = n1 * n2 + (n1 * (n1 + 1)) / 2.0 - np.sum(rankx, axis=0)
-    bigu = max(u1,u2)
+    u2 = n1 * n2 - u1                            # remainder is U for y
-    smallu = min(u1,u2)
+    bigu = max(u1, u2)
+    smallu = min(u1, u2)
     T = tiecorrect(ranked)
     if T == 0:
         raise ValueError('All numbers are identical in amannwhitneyu')
-    sd = np.sqrt(T*n1*n2*(n1+n2+1)/12.0)
+    sd = np.sqrt(T * n1 * n2 * (n1 + n2 + 1) / 12.0)
 
     if use_continuity:
         # normal approximation for prob calc with continuity correction
-        z = abs((bigu-0.5-n1*n2/2.0) / sd)
+        z = abs((bigu - 0.5 - n1 * n2 / 2.0) / sd)
     else:
-        z = abs((bigu-n1*n2/2.0) / sd)  # normal approximation for prob calc
+        # normal approximation for prob calc
+        z = abs((bigu - n1 * n2 / 2.0) / sd)
     return smallu, distributions.norm.sf(z)  # (1.0 - zprob(z))
 
 
@@ -3966,16 +3980,16 @@ def ranksums(x, y):
     .. [1] http://en.wikipedia.org/wiki/Wilcoxon_rank-sum_test
 
     """
-    x,y = map(np.asarray, (x, y))
+    x, y = map(np.asarray, (x, y))
     n1 = len(x)
     n2 = len(y)
-    alldata = np.concatenate((x,y))
+    alldata = np.concatenate((x, y))
     ranked = rankdata(alldata)
     x = ranked[:n1]
     y = ranked[n1:]
-    s = np.sum(x,axis=0)
+    s = np.sum(x, axis=0)
-    expected = n1*(n1+n2+1) / 2.0
+    expected = n1 * (n1 + n2 + 1) / 2.0
-    z = (s - expected) / np.sqrt(n1*n2*(n1+n2+1)/12.0)
+    z = (s - expected) / np.sqrt(n1 * n2 * (n1 + n2 + 1) / 12.0)
     prob = 2 * distributions.norm.sf(abs(z))
     return z, prob
 
@@ -4033,7 +4047,7 @@ def kruskal(*args):
     j = np.insert(np.cumsum(n), 0, 0)
     ssbn = 0
     for i in range(na):
-        ssbn += square_of_sums(ranked[j[i]:j[i+1]]) / float(n[i])
+        ssbn += square_of_sums(ranked[j[i]:j[i + 1]]) / float(n[i])
 
     totaln = np.sum(n)
     h = 12.0 / (totaln * (totaln + 1)) * ssbn - 3 * (totaln + 1)
@@ -4080,7 +4094,8 @@ def friedmanchisquare(*args):
     """
     k = len(args)
     if k < 3:
-        raise ValueError('\nLess than 3 levels.  Friedman test not appropriate.\n')
+        raise ValueError(
+            '\nLess than 3 levels.  Friedman test not appropriate.\n')
 
     n = len(args[0])
     for i in range(1, k):
@@ -4096,19 +4111,19 @@ def friedmanchisquare(*args):
     # Handle ties
     ties = 0
     for i in range(len(data)):
-        replist, repnum = find_repeats(array(data[i]))
+        _replist, repnum = find_repeats(array(data[i]))
         for t in repnum:
-            ties += t*(t*t-1)
+            ties += t * (t * t - 1)
-    c = 1 - ties / float(k*(k*k-1)*n)
+    c = 1 - ties / float(k * (k * k - 1) * n)
 
-    ssbn = pysum(pysum(data)**2)
+    ssbn = pysum(pysum(data) ** 2)
-    chisq = (12.0 / (k*n*(k+1)) * ssbn - 3*n*(k+1)) / c
+    chisq = (12.0 / (k * n * (k + 1)) * ssbn - 3 * n * (k + 1)) / c
-    return chisq, chisqprob(chisq,k-1)
+    return chisq, chisqprob(chisq, k - 1)
 
 
-#####################################
+#
-####  PROBABILITY CALCULATIONS  ####
+# PROBABILITY CALCULATIONS  ####
-#####################################
+#
 
 zprob = special.ndtr
 
@@ -4132,7 +4147,7 @@ def chisqprob(chisq, df):
         distribution with degrees of freedom `df`.
 
     """
-    return special.chdtrc(df,chisq)
+    return special.chdtrc(df, chisq)
 
 ksprob = special.kolmogorov
 fprob = special.fdtrc
@@ -4169,9 +4184,9 @@ def betai(a, b, x):
     return special.betainc(a, b, x)
 
 
-#####################################
+#
-#######  ANOVA CALCULATIONS  #######
+# ANOVA CALCULATIONS  #######
-#####################################
+#
 
 def f_value_wilks_lambda(ER, EF, dfnum, dfden, a, b):
     """Calculation of Wilks lambda F-statistic for multivarite data, per
@@ -4182,12 +4197,13 @@ def f_value_wilks_lambda(ER, EF, dfnum, dfden, a, b):
     if isinstance(EF, (int, float)):
         EF = array([[EF]])
     lmbda = linalg.det(EF) / linalg.det(ER)
-    if (a-1)**2 + (b-1)**2 == 5:
+    if (a - 1) ** 2 + (b - 1) ** 2 == 5:
         q = 1
     else:
-        q = np.sqrt(((a-1)**2*(b-1)**2 - 2) / ((a-1)**2 + (b-1)**2 - 5))
+        q = np.sqrt(
-    n_um = (1 - lmbda**(1.0/q))*(a-1)*(b-1)
+            ((a - 1) ** 2 * (b - 1) ** 2 - 2) / ((a - 1) ** 2 + (b - 1) ** 2 - 5))
-    d_en = lmbda**(1.0/q) / (n_um*q - 0.5*(a-1)*(b-1) + 1)
+    n_um = (1 - lmbda ** (1.0 / q)) * (a - 1) * (b - 1)
+    d_en = lmbda ** (1.0 / q) / (n_um * q - 0.5 * (a - 1) * (b - 1) + 1)
     return n_um / d_en
 
 
@@ -4216,7 +4232,7 @@ def f_value(ER, EF, dfR, dfF):
     F-statistic : float
 
     """
-    return ((ER-EF)/float(dfR-dfF) / (EF/float(dfF)))
+    return ((ER - EF) / float(dfR - dfF) / (EF / float(dfF)))
 
 
 def f_value_multivariate(ER, EF, dfnum, dfden):
@@ -4251,9 +4267,9 @@ def f_value_multivariate(ER, EF, dfnum, dfden):
     return n_um / d_en
 
 
-#####################################
+#
-#######  SUPPORT FUNCTIONS  ########
+# SUPPORT FUNCTIONS  ########
-#####################################
+#
 
 def ss(a, axis=0):
     """
@@ -4291,7 +4307,7 @@ def ss(a, axis=0):
 
     """
     a, axis = _chk_asarray(a, axis)
-    return np.sum(a*a, axis)
+    return np.sum(a * a, axis)
 
 
 def square_of_sums(a, axis=0):
@@ -4326,11 +4342,11 @@ def square_of_sums(a, axis=0):
 
     """
     a, axis = _chk_asarray(a, axis)
-    s = np.sum(a,axis)
+    s = np.sum(a, axis)
     if not np.isscalar(s):
-        return s.astype(float)*s
+        return s.astype(float) * s
     else:
-        return float(s)*s
+        return float(s) * s
 
 
 def fastsort(a):

pywafo/src/wafo/stats/tests/common_tests.py

@@ -11,7 +11,7 @@ from scipy import stats
 
 
 #NUMPY_BELOW_1_7 = NumpyVersion(np.__version__) < '1.7.0'
-NUMPY_BELOW_1_7 =np.__version__ < '1.7.0'
+NUMPY_BELOW_1_7 = np.__version__ < '1.7.0'
 
 
 def check_normalization(distfn, args, distname):
@@ -60,14 +60,14 @@ def check_mean_expect(distfn, arg, m, msg):
 
 def check_var_expect(distfn, arg, m, v, msg):
     if np.isfinite(v):
-        m2 = distfn.expect(lambda x: x*x, arg)
+        m2 = distfn.expect(lambda x: x * x, arg)
-        npt.assert_almost_equal(m2, v + m*m, decimal=5, err_msg=msg +
+        npt.assert_almost_equal(m2, v + m * m, decimal=5, err_msg=msg +
                                 ' - 2st moment (expect)')
 
 
 def check_skew_expect(distfn, arg, m, v, s, msg):
     if np.isfinite(s):
-        m3e = distfn.expect(lambda x: np.power(x-m, 3), arg)
+        m3e = distfn.expect(lambda x: np.power(x - m, 3), arg)
         npt.assert_almost_equal(m3e, s * np.power(v, 1.5),
                                 decimal=5, err_msg=msg + ' - skew')
     else:
@@ -76,8 +76,9 @@ def check_skew_expect(distfn, arg, m, v, s, msg):
 
 def check_kurt_expect(distfn, arg, m, v, k, msg):
     if np.isfinite(k):
-        m4e = distfn.expect(lambda x: np.power(x-m, 4), arg)
+        m4e = distfn.expect(lambda x: np.power(x - m, 4), arg)
-        npt.assert_allclose(m4e, (k + 3.) * np.power(v, 2), atol=1e-5, rtol=1e-5,
+        npt.assert_allclose(
+            m4e, (k + 3.) * np.power(v, 2), atol=1e-5, rtol=1e-5,
             err_msg=msg + ' - kurtosis')
     else:
         npt.assert_(np.isnan(k))
@@ -115,7 +116,7 @@ def check_edge_support(distfn, args):
 
 
 def check_named_args(distfn, x, shape_args, defaults, meths):
-    ## Check calling w/ named arguments.
+    # Check calling w/ named arguments.
 
     # check consistency of shapes, numargs and _parse signature
     signature = inspect.getargspec(distfn._parse_args)
@@ -123,7 +124,8 @@ def check_named_args(distfn, x, shape_args, defaults, meths):
     npt.assert_(signature.keywords is None)
     npt.assert_(signature.defaults == defaults)
 
-    shape_argnames = signature.args[1:-len(defaults)]  # self, a, b, loc=0, scale=1
+    # self, a, b, loc=0, scale=1
+    shape_argnames = signature.args[1:-len(defaults)]
     if distfn.shapes:
         shapes_ = distfn.shapes.replace(',', ' ').split()
     else:

pywafo/src/wafo/stats/tests/test_binned_statistic.py

@@ -1,5 +1,4 @@
 from __future__ import division, print_function, absolute_import
-
 import numpy as np
 from numpy.testing import assert_array_almost_equal, run_module_suite
 from scipy.stats import \
@@ -235,4 +234,5 @@ class TestBinnedStatistic(object):
 
 
 if __name__ == "__main__":
+    #unittest.main()
     run_module_suite()

pywafo/src/wafo/stats/tests/test_continuous_basic.py

@@ -6,8 +6,9 @@ import numpy as np
 import numpy.testing as npt
 
 from scipy import integrate
-from scipy import stats
+from wafo import stats
-from common_tests import (check_normalization, check_moment, check_mean_expect,
+from wafo.stats.tests.common_tests import (check_normalization, check_moment,
+                           check_mean_expect,
         check_var_expect, check_skew_expect, check_kurt_expect,
         check_entropy, check_private_entropy, NUMPY_BELOW_1_7,
         check_edge_support, check_named_args)
@@ -158,7 +159,7 @@ distmissing = ['wald', 'gausshyper', 'genexpon', 'rv_continuous',
     'johnsonsb', 'truncexpon', 'rice', 'invgauss', 'invgamma',
     'powerlognorm']
 
-distmiss = [[dist,args] for dist,args in distcont if dist in distmissing]
+distmiss = [[dist, args] for dist, args in distcont if dist in distmissing]
 distslow = ['rdist', 'gausshyper', 'recipinvgauss', 'ksone', 'genexpon',
             'vonmises', 'vonmises_line', 'mielke', 'semicircular',
             'cosine', 'invweibull', 'powerlognorm', 'johnsonsu', 'kstwobign']
@@ -181,7 +182,8 @@ def _silence_fp_errors(func):
 def test_cont_basic():
     # this test skips slow distributions
     with warnings.catch_warnings():
-        warnings.filterwarnings('ignore', category=integrate.IntegrationWarning)
+#         warnings.filterwarnings('ignore',
+#                                 category=integrate.IntegrationWarning)
         for distname, arg in distcont[:]:
             if distname in distslow:
                 continue
@@ -233,14 +235,15 @@ def test_cont_basic():
 def test_cont_basic_slow():
     # same as above for slow distributions
     with warnings.catch_warnings():
-        warnings.filterwarnings('ignore', category=integrate.IntegrationWarning)
+#         warnings.filterwarnings('ignore',
+#                                 category=integrate.IntegrationWarning)
         for distname, arg in distcont[:]:
             if distname not in distslow:
                 continue
             distfn = getattr(stats, distname)
             np.random.seed(765456)
             sn = 500
-            rvs = distfn.rvs(size=sn,*arg)
+            rvs = distfn.rvs(size=sn, *arg)
             sm = rvs.mean()
             sv = rvs.var()
             m, v = distfn.stats(*arg)
@@ -284,7 +287,8 @@ def test_cont_basic_slow():
 @npt.dec.slow
 def test_moments():
     with warnings.catch_warnings():
-        warnings.filterwarnings('ignore', category=integrate.IntegrationWarning)
+#         warnings.filterwarnings('ignore',
+#                                 category=integrate.IntegrationWarning)
         knf = npt.dec.knownfailureif
         fail_normalization = set(['vonmises', 'ksone'])
         fail_higher = set(['vonmises', 'ksone', 'ncf'])
@@ -312,44 +316,45 @@ def check_sample_meanvar_(distfn, arg, m, v, sm, sv, sn, msg):
         check_sample_var(sv, sn, v)
 
 
-def check_sample_mean(sm,v,n, popmean):
+def check_sample_mean(sm, v, n, popmean):
     # from stats.stats.ttest_1samp(a, popmean):
     # Calculates the t-obtained for the independent samples T-test on ONE group
     # of scores a, given a population mean.
     #
     # Returns: t-value, two-tailed prob
-    df = n-1
+    df = n - 1
-    svar = ((n-1)*v) / float(df)    # looks redundant
+    svar = ((n - 1) * v) / float(df)    # looks redundant
-    t = (sm-popmean) / np.sqrt(svar*(1.0/n))
+    t = (sm - popmean) / np.sqrt(svar * (1.0 / n))
-    prob = stats.betai(0.5*df, 0.5, df/(df+t*t))
+    prob = stats.betai(0.5 * df, 0.5, df / (df + t * t))
 
     # return t,prob
     npt.assert_(prob > 0.01, 'mean fail, t,prob = %f, %f, m, sm=%f,%f' %
                 (t, prob, popmean, sm))
 
 
-def check_sample_var(sv,n, popvar):
+def check_sample_var(sv, n, popvar):
-    # two-sided chisquare test for sample variance equal to hypothesized variance
+    # two-sided chisquare test for sample variance equal to hypothesized
-    df = n-1
+    # variance
-    chi2 = (n-1)*popvar/float(popvar)
+    df = n - 1
-    pval = stats.chisqprob(chi2,df)*2
+    chi2 = (n - 1) * popvar / float(popvar)
+    pval = stats.chisqprob(chi2, df) * 2
     npt.assert_(pval > 0.01, 'var fail, t, pval = %f, %f, v, sv=%f, %f' %
-            (chi2,pval,popvar,sv))
+                (chi2, pval, popvar, sv))
 
 
-def check_cdf_ppf(distfn,arg,msg):
+def check_cdf_ppf(distfn, arg, msg):
     values = [0.001, 0.5, 0.999]
     npt.assert_almost_equal(distfn.cdf(distfn.ppf(values, *arg), *arg),
                             values, decimal=DECIMAL, err_msg=msg +
                             ' - cdf-ppf roundtrip')
 
 
-def check_sf_isf(distfn,arg,msg):
+def check_sf_isf(distfn, arg, msg):
-    npt.assert_almost_equal(distfn.sf(distfn.isf([0.1,0.5,0.9], *arg), *arg),
+    npt.assert_almost_equal(distfn.sf(distfn.isf([0.1, 0.5, 0.9], *arg), *arg),
-                            [0.1,0.5,0.9], decimal=DECIMAL, err_msg=msg +
+                            [0.1, 0.5, 0.9], decimal=DECIMAL, err_msg=msg +
                             ' - sf-isf roundtrip')
-    npt.assert_almost_equal(distfn.cdf([0.1,0.9], *arg),
+    npt.assert_almost_equal(distfn.cdf([0.1, 0.9], *arg),
-                            1.0-distfn.sf([0.1,0.9], *arg),
+                            1.0 - distfn.sf([0.1, 0.9], *arg),
                             decimal=DECIMAL, err_msg=msg +
                             ' - cdf-sf relationship')
 
@@ -360,15 +365,16 @@ def check_pdf(distfn, arg, msg):
     eps = 1e-6
     pdfv = distfn.pdf(median, *arg)
     if (pdfv < 1e-4) or (pdfv > 1e4):
-        # avoid checking a case where pdf is close to zero or huge (singularity)
+        # avoid checking a case where pdf is close to zero or huge
+        # (singularity)
         median = median + 0.1
         pdfv = distfn.pdf(median, *arg)
     cdfdiff = (distfn.cdf(median + eps, *arg) -
-               distfn.cdf(median - eps, *arg))/eps/2.0
+               distfn.cdf(median - eps, *arg)) / eps / 2.0
     # replace with better diff and better test (more points),
     # actually, this works pretty well
-    npt.assert_almost_equal(pdfv, cdfdiff,
+    npt.assert_almost_equal(pdfv, cdfdiff, decimal=DECIMAL,
-                decimal=DECIMAL, err_msg=msg + ' - cdf-pdf relationship')
+                            err_msg=msg + ' - cdf-pdf relationship')
 
 
 def check_pdf_logpdf(distfn, args, msg):
@@ -379,7 +385,8 @@ def check_pdf_logpdf(distfn, args, msg):
     logpdf = distfn.logpdf(vals, *args)
     pdf = pdf[pdf != 0]
     logpdf = logpdf[np.isfinite(logpdf)]
-    npt.assert_almost_equal(np.log(pdf), logpdf, decimal=7, err_msg=msg + " - logpdf-log(pdf) relationship")
+    npt.assert_almost_equal(np.log(pdf), logpdf, decimal=7,
+                            err_msg=msg + " - logpdf-log(pdf) relationship")
 
 
 def check_sf_logsf(distfn, args, msg):
@@ -390,7 +397,8 @@ def check_sf_logsf(distfn, args, msg):
     logsf = distfn.logsf(vals, *args)
     sf = sf[sf != 0]
     logsf = logsf[np.isfinite(logsf)]
-    npt.assert_almost_equal(np.log(sf), logsf, decimal=7, err_msg=msg + " - logsf-log(sf) relationship")
+    npt.assert_almost_equal(np.log(sf), logsf, decimal=7,
+                            err_msg=msg + " - logsf-log(sf) relationship")
 
 
 def check_cdf_logcdf(distfn, args, msg):
@@ -401,15 +409,16 @@ def check_cdf_logcdf(distfn, args, msg):
     logcdf = distfn.logcdf(vals, *args)
     cdf = cdf[cdf != 0]
     logcdf = logcdf[np.isfinite(logcdf)]
-    npt.assert_almost_equal(np.log(cdf), logcdf, decimal=7, err_msg=msg + " - logcdf-log(cdf) relationship")
+    npt.assert_almost_equal(np.log(cdf), logcdf, decimal=7,
+        err_msg=msg + " - logcdf-log(cdf) relationship")
 
 
 def check_distribution_rvs(dist, args, alpha, rvs):
     # test from scipy.stats.tests
     # this version reuses existing random variables
-    D,pval = stats.kstest(rvs, dist, args=args, N=1000)
+    D, pval = stats.kstest(rvs, dist, args=args, N=1000)
     if (pval < alpha):
-        D,pval = stats.kstest(dist,'',args=args, N=1000)
+        D, pval = stats.kstest(dist, '', args=args, N=1000)
         npt.assert_(pval > alpha, "D = " + str(D) + "; pval = " + str(pval) +
                     "; alpha = " + str(alpha) + "\nargs = " + str(args))
 
@@ -417,16 +426,17 @@ def check_distribution_rvs(dist, args, alpha, rvs):
 def check_vecentropy(distfn, args):
     npt.assert_equal(distfn.vecentropy(*args), distfn._entropy(*args))
 
+
 @npt.dec.skipif(NUMPY_BELOW_1_7)
 def check_loc_scale(distfn, arg, m, v, msg):
     loc, scale = 10.0, 10.0
     mt, vt = distfn.stats(loc=loc, scale=scale, *arg)
-    npt.assert_allclose(m*scale + loc, mt)
+    npt.assert_allclose(m * scale + loc, mt)
-    npt.assert_allclose(v*scale*scale, vt)
+    npt.assert_allclose(v * scale * scale, vt)
 
 
 def check_ppf_private(distfn, arg, msg):
-    #fails by design for truncnorm self.nb not defined
+    # fails by design for truncnorm self.nb not defined
     ppfs = distfn._ppf(np.array([0.1, 0.5, 0.9]), *arg)
     npt.assert_(not np.any(np.isnan(ppfs)), msg + 'ppf private is nan')
 

pywafo/src/wafo/stats/tests/test_discrete_basic.py

@@ -3,25 +3,26 @@ from __future__ import division, print_function, absolute_import
 import numpy.testing as npt
 import numpy as np
 try:
-    from scipy.lib.six import xrange
+    from wafo.stats.six import xrange
 except:
     pass
-from scipy import stats
+from wafo import stats
-from .common_tests import (check_normalization, check_moment, check_mean_expect,
+from wafo.stats.tests.common_tests import (check_normalization, check_moment,
+                                           check_mean_expect,
         check_var_expect, check_skew_expect, check_kurt_expect,
         check_entropy, check_private_entropy, check_edge_support,
         check_named_args)
 knf = npt.dec.knownfailureif
 
 distdiscrete = [
-    ['bernoulli',(0.3,)],
+    ['bernoulli', (0.3, )],
     ['binom', (5, 0.4)],
-    ['boltzmann',(1.4, 19)],
+    ['boltzmann', (1.4, 19)],
     ['dlaplace', (0.8,)],  # 0.5
     ['geom', (0.5,)],
-    ['hypergeom',(30, 12, 6)],
+    ['hypergeom', (30, 12, 6)],
-    ['hypergeom',(21,3,12)],  # numpy.random (3,18,12) numpy ticket:921
+    ['hypergeom', (21, 3, 12)],  # numpy.random (3,18,12) numpy ticket:921
-    ['hypergeom',(21,18,11)],  # numpy.random (18,3,11) numpy ticket:921
+    ['hypergeom', (21, 18, 11)],  # numpy.random (18,3,11) numpy ticket:921
     ['logser', (0.6,)],  # reenabled, numpy ticket:921
     ['nbinom', (5, 0.5)],
     ['nbinom', (0.4, 0.4)],  # from tickets: 583
@@ -39,7 +40,7 @@ def test_discrete_basic():
         np.random.seed(9765456)
         rvs = distfn.rvs(size=2000, *arg)
         supp = np.unique(rvs)
-        m, v = distfn.stats(*arg)
+        #_m, v = distfn.stats(*arg)
         yield check_cdf_ppf, distfn, arg, supp, distname + ' cdf_ppf'
 
         yield check_pmf_cdf, distfn, arg, distname
@@ -55,7 +56,7 @@ def test_discrete_basic():
         if distname in seen:
             continue
         seen.add(distname)
-        distfn = getattr(stats,distname)
+        distfn = getattr(stats, distname)
         locscale_defaults = (0,)
         meths = [distfn.pmf, distfn.logpmf, distfn.cdf, distfn.logcdf,
                  distfn.logsf]
@@ -73,7 +74,7 @@ def test_discrete_basic():
 
 def test_moments():
     for distname, arg in distdiscrete:
-        distfn = getattr(stats,distname)
+        distfn = getattr(stats, distname)
         m, v, s, k = distfn.stats(*arg, moments='mvsk')
         yield check_normalization, distfn, arg, distname
 
@@ -89,7 +90,7 @@ def test_moments():
 
         # frozen distr moments
         yield check_moment_frozen, distfn, arg, m, 1
-        yield check_moment_frozen, distfn, arg, v+m*m, 2
+        yield check_moment_frozen, distfn, arg, v + m * m, 2
 
 
 def check_cdf_ppf(distfn, arg, supp, msg):
@@ -107,7 +108,7 @@ def check_cdf_ppf(distfn, arg, supp, msg):
 def check_pmf_cdf(distfn, arg, distname):
     startind = np.int(distfn.ppf(0.01, *arg) - 1)
     index = list(range(startind, startind + 10))
-    cdfs, pmfs_cum = distfn.cdf(index,*arg), distfn.pmf(index, *arg).cumsum()
+    cdfs, pmfs_cum = distfn.cdf(index, *arg), distfn.pmf(index, *arg).cumsum()
 
     atol, rtol = 1e-10, 1e-10
     if distname == 'skellam':    # ncx2 accuracy
@@ -157,7 +158,7 @@ def check_discrete_chisquare(distfn, arg, rvs, alpha, msg):
     """
     n = len(rvs)
     nsupp = 20
-    wsupp = 1.0/nsupp
+    wsupp = 1.0 / nsupp
 
     # construct intervals with minimum mass 1/nsupp
     # intervals are left-half-open as in a cdf difference
@@ -166,30 +167,30 @@ def check_discrete_chisquare(distfn, arg, rvs, alpha, msg):
     distsupp = [max(distfn.a, -1000)]
     distmass = []
     for ii in distsupport:
-        current = distfn.cdf(ii,*arg)
+        current = distfn.cdf(ii, *arg)
-        if current - last >= wsupp-1e-14:
+        if current - last >= wsupp - 1e-14:
             distsupp.append(ii)
             distmass.append(current - last)
             last = current
-            if current > (1-wsupp):
+            if current > (1 - wsupp):
                 break
     if distsupp[-1] < distfn.b:
         distsupp.append(distfn.b)
-        distmass.append(1-last)
+        distmass.append(1 - last)
     distsupp = np.array(distsupp)
     distmass = np.array(distmass)
 
     # convert intervals to right-half-open as required by histogram
-    histsupp = distsupp+1e-8
+    histsupp = distsupp + 1e-8
     histsupp[0] = distfn.a
 
     # find sample frequencies and perform chisquare test
-    freq,hsupp = np.histogram(rvs,histsupp)
+    freq, _hsupp = np.histogram(rvs, histsupp)
-    cdfs = distfn.cdf(distsupp,*arg)
+    #cdfs = distfn.cdf(distsupp, *arg)
-    (chis,pval) = stats.chisquare(np.array(freq),n*distmass)
+    (_chis, pval) = stats.chisquare(np.array(freq), n * distmass)
 
     npt.assert_(pval > alpha, 'chisquare - test for %s'
-           ' at arg = %s with pval = %s' % (msg,str(arg),str(pval)))
+           ' at arg = %s with pval = %s' % (msg, str(arg), str(pval)))
 
 
 def check_scale_docstring(distfn):

pywafo/src/wafo/stats/tests/test_fit.py

@@ -7,7 +7,7 @@ from numpy.testing import dec
 
 from wafo import stats
 
-from .test_continuous_basic import distcont
+from wafo.stats.tests.test_continuous_basic import distcont
 
 # this is not a proper statistical test for convergence, but only
 # verifies that the estimate and true values don't differ by too much
@@ -45,14 +45,14 @@ skip_fit = [
 def test_cont_fit():
     # this tests the closeness of the estimated parameters to the true
     # parameters with fit method of continuous distributions
-    # Note: is slow, some distributions don't converge with sample size <= 10000
+    # Note: slow, some distributions don't converge with sample size <= 10000
 
     for distname, arg in distcont:
         if distname not in skip_fit:
-            yield check_cont_fit, distname,arg
+            yield check_cont_fit, distname, arg
 
 
-def check_cont_fit(distname,arg):
+def check_cont_fit(distname, arg):
     if distname in failing_fits:
         # Skip failing fits unless overridden
         xfail = True
@@ -62,14 +62,16 @@ def check_cont_fit(distname,arg):
             pass
         if xfail:
             msg = "Fitting %s doesn't work reliably yet" % distname
-            msg += " [Set environment variable SCIPY_XFAIL=1 to run this test nevertheless.]"
+            msg += " [Set environment variable SCIPY_XFAIL=1 to run this " + \
+                "test nevertheless.]"
             dec.knownfailureif(True, msg)(lambda: None)()
 
     distfn = getattr(stats, distname)
 
-    truearg = np.hstack([arg,[0.0,1.0]])
+    truearg = np.hstack([arg, [0.0, 1.0]])
-    diffthreshold = np.max(np.vstack([truearg*thresh_percent,
+    diffthreshold = np.max(np.vstack([
-                                      np.ones(distfn.numargs+2)*thresh_min]),0)
+        truearg * thresh_percent,
+        np.ones(distfn.numargs + 2) * thresh_min]), 0)
 
     for fit_size in fit_sizes:
         # Note that if a fit succeeds, the other fit_sizes are skipped
@@ -77,12 +79,16 @@ def check_cont_fit(distname,arg):
 
         with np.errstate(all='ignore'):
             rvs = distfn.rvs(size=fit_size, *arg)
-            est = distfn.fit(rvs)  # start with default values
+            #phat = distfn.fit2(rvs)
+            phat = distfn.fit2(rvs, method='mps')
+            est = phat.par
+            #est = distfn.fit(rvs)  # start with default values
 
         diff = est - truearg
 
         # threshold for location
-        diffthreshold[-2] = np.max([np.abs(rvs.mean())*thresh_percent,thresh_min])
+        diffthreshold[-2] = np.max([np.abs(rvs.mean()) * thresh_percent,
+                                    thresh_min])
 
         if np.any(np.isnan(est)):
             raise AssertionError('nan returned in fit')

pywafo/src/wafo/stats/tests/test_morestats.py

@@ -9,10 +9,11 @@ import warnings
 import numpy as np
 from numpy.random import RandomState
 from numpy.testing import (TestCase, run_module_suite, assert_array_equal,
-    assert_almost_equal, assert_array_less, assert_array_almost_equal,
+                           assert_almost_equal, assert_array_less,
-    assert_raises, assert_, assert_allclose, assert_equal, dec)
+                           assert_array_almost_equal, assert_raises, assert_,
+                           assert_allclose, assert_equal, dec)
 
-from scipy import stats
+from wafo import stats
 
 # Matplotlib is not a scipy dependency but is optionally used in probplot, so
 # check if it's available
@@ -36,19 +37,20 @@ g10 = [0.991, 0.995, 0.984, 0.994, 0.997, 0.997, 0.991, 0.998, 1.004, 0.997]
 
 
 class TestShapiro(TestCase):
+
     def test_basic(self):
-        x1 = [0.11,7.87,4.61,10.14,7.95,3.14,0.46,
+        x1 = [0.11, 7.87, 4.61, 10.14, 7.95, 3.14, 0.46,
-              4.43,0.21,4.75,0.71,1.52,3.24,
+              4.43, 0.21, 4.75, 0.71, 1.52, 3.24,
-              0.93,0.42,4.97,9.53,4.55,0.47,6.66]
+              0.93, 0.42, 4.97, 9.53, 4.55, 0.47, 6.66]
-        w,pw = stats.shapiro(x1)
+        w, pw = stats.shapiro(x1)
-        assert_almost_equal(w,0.90047299861907959,6)
+        assert_almost_equal(w, 0.90047299861907959, 6)
-        assert_almost_equal(pw,0.042089745402336121,6)
+        assert_almost_equal(pw, 0.042089745402336121, 6)
-        x2 = [1.36,1.14,2.92,2.55,1.46,1.06,5.27,-1.11,
+        x2 = [1.36, 1.14, 2.92, 2.55, 1.46, 1.06, 5.27, -1.11,
-              3.48,1.10,0.88,-0.51,1.46,0.52,6.20,1.69,
+              3.48, 1.10, 0.88, -0.51, 1.46, 0.52, 6.20, 1.69,
-              0.08,3.67,2.81,3.49]
+              0.08, 3.67, 2.81, 3.49]
-        w,pw = stats.shapiro(x2)
+        w, pw = stats.shapiro(x2)
-        assert_almost_equal(w,0.9590270,6)
+        assert_almost_equal(w, 0.9590270, 6)
-        assert_almost_equal(pw,0.52460,3)
+        assert_almost_equal(pw, 0.52460, 3)
 
     def test_bad_arg(self):
         # Length of x is less than 3.
@@ -57,24 +59,25 @@ class TestShapiro(TestCase):
 
 
 class TestAnderson(TestCase):
+
     def test_normal(self):
         rs = RandomState(1234567890)
         x1 = rs.standard_exponential(size=50)
         x2 = rs.standard_normal(size=50)
-        A,crit,sig = stats.anderson(x1)
+        A, crit, _sig = stats.anderson(x1)
         assert_array_less(crit[:-1], A)
-        A,crit,sig = stats.anderson(x2)
+        A, crit, _sig = stats.anderson(x2)
         assert_array_less(A, crit[-2:])
 
     def test_expon(self):
         rs = RandomState(1234567890)
         x1 = rs.standard_exponential(size=50)
         x2 = rs.standard_normal(size=50)
-        A,crit,sig = stats.anderson(x1,'expon')
+        A, crit, _sig = stats.anderson(x1, 'expon')
         assert_array_less(A, crit[-2:])
         olderr = np.seterr(all='ignore')
         try:
-            A,crit,sig = stats.anderson(x2,'expon')
+            A, crit, _sig = stats.anderson(x2, 'expon')
         finally:
             np.seterr(**olderr)
         assert_(A > crit[-1])
@@ -86,30 +89,31 @@ class TestAnderson(TestCase):
 class TestAnsari(TestCase):
 
     def test_small(self):
-        x = [1,2,3,3,4]
+        x = [1, 2, 3, 3, 4]
-        y = [3,2,6,1,6,1,4,1]
+        y = [3, 2, 6, 1, 6, 1, 4, 1]
-        W, pval = stats.ansari(x,y)
+        W, pval = stats.ansari(x, y)
-        assert_almost_equal(W,23.5,11)
+        assert_almost_equal(W, 23.5, 11)
-        assert_almost_equal(pval,0.13499256881897437,11)
+        assert_almost_equal(pval, 0.13499256881897437, 11)
 
     def test_approx(self):
         ramsay = np.array((111, 107, 100, 99, 102, 106, 109, 108, 104, 99,
                            101, 96, 97, 102, 107, 113, 116, 113, 110, 98))
-        parekh = np.array((107, 108, 106, 98, 105, 103, 110, 105, 104,
+        parekh = np.array((107, 108, 106, 98, 105, 103, 110, 105, 104, 100,
-                           100, 96, 108, 103, 104, 114, 114, 113, 108, 106, 99))
+                           96, 108, 103, 104, 114, 114, 113, 108, 106, 99))
 
         with warnings.catch_warnings():
             warnings.filterwarnings('ignore',
-                        message="Ties preclude use of exact statistic.")
+                                    message="Ties preclude use of exact " +
+                                    "statistic.")
             W, pval = stats.ansari(ramsay, parekh)
 
-        assert_almost_equal(W,185.5,11)
+        assert_almost_equal(W, 185.5, 11)
-        assert_almost_equal(pval,0.18145819972867083,11)
+        assert_almost_equal(pval, 0.18145819972867083, 11)
 
     def test_exact(self):
-        W,pval = stats.ansari([1,2,3,4],[15,5,20,8,10,12])
+        W, pval = stats.ansari([1, 2, 3, 4], [15, 5, 20, 8, 10, 12])
-        assert_almost_equal(W,10.0,11)
+        assert_almost_equal(W, 10.0, 11)
-        assert_almost_equal(pval,0.533333333333333333,7)
+        assert_almost_equal(pval, 0.533333333333333333, 7)
 
     def test_bad_arg(self):
         assert_raises(ValueError, stats.ansari, [], [1])
@@ -121,8 +125,8 @@ class TestBartlett(TestCase):
     def test_data(self):
         args = [g1, g2, g3, g4, g5, g6, g7, g8, g9, g10]
         T, pval = stats.bartlett(*args)
-        assert_almost_equal(T,20.78587342806484,7)
+        assert_almost_equal(T, 20.78587342806484, 7)
-        assert_almost_equal(pval,0.0136358632781,7)
+        assert_almost_equal(pval, 0.0136358632781, 7)
 
     def test_bad_arg(self):
         # Too few args raises ValueError.
@@ -134,14 +138,15 @@ class TestLevene(TestCase):
     def test_data(self):
         args = [g1, g2, g3, g4, g5, g6, g7, g8, g9, g10]
         W, pval = stats.levene(*args)
-        assert_almost_equal(W,1.7059176930008939,7)
+        assert_almost_equal(W, 1.7059176930008939, 7)
-        assert_almost_equal(pval,0.0990829755522,7)
+        assert_almost_equal(pval, 0.0990829755522, 7)
 
     def test_trimmed1(self):
         # Test that center='trimmed' gives the same result as center='mean'
         # when proportiontocut=0.
         W1, pval1 = stats.levene(g1, g2, g3, center='mean')
-        W2, pval2 = stats.levene(g1, g2, g3, center='trimmed', proportiontocut=0.0)
+        W2, pval2 = stats.levene(
+            g1, g2, g3, center='trimmed', proportiontocut=0.0)
         assert_almost_equal(W1, W2)
         assert_almost_equal(pval1, pval2)
 
@@ -152,8 +157,10 @@ class TestLevene(TestCase):
         x2 = np.random.permutation(x)
 
         # Use center='trimmed'
-        W0, pval0 = stats.levene(x, y, center='trimmed', proportiontocut=0.125)
+        W0, _pval0 = stats.levene(x, y, center='trimmed',
-        W1, pval1 = stats.levene(x2, y, center='trimmed', proportiontocut=0.125)
+                                  proportiontocut=0.125)
+        W1, pval1 = stats.levene(
+            x2, y, center='trimmed', proportiontocut=0.125)
         # Trim the data here, and use center='mean'
         W2, pval2 = stats.levene(x[1:-1], y[1:-1], center='mean')
         # Result should be the same.
@@ -162,21 +169,21 @@ class TestLevene(TestCase):
         assert_almost_equal(pval1, pval2)
 
     def test_equal_mean_median(self):
-        x = np.linspace(-1,1,21)
+        x = np.linspace(-1, 1, 21)
         np.random.seed(1234)
         x2 = np.random.permutation(x)
-        y = x**3
+        y = x ** 3
         W1, pval1 = stats.levene(x, y, center='mean')
         W2, pval2 = stats.levene(x2, y, center='median')
         assert_almost_equal(W1, W2)
         assert_almost_equal(pval1, pval2)
 
     def test_bad_keyword(self):
-        x = np.linspace(-1,1,21)
+        x = np.linspace(-1, 1, 21)
         assert_raises(TypeError, stats.levene, x, x, portiontocut=0.1)
 
     def test_bad_center_value(self):
-        x = np.linspace(-1,1,21)
+        x = np.linspace(-1, 1, 21)
         assert_raises(ValueError, stats.levene, x, x, center='trim')
 
     def test_too_few_args(self):
@@ -186,16 +193,16 @@ class TestLevene(TestCase):
 class TestBinomP(TestCase):
 
     def test_data(self):
-        pval = stats.binom_test(100,250)
+        pval = stats.binom_test(100, 250)
-        assert_almost_equal(pval,0.0018833009350757682,11)
+        assert_almost_equal(pval, 0.0018833009350757682, 11)
-        pval = stats.binom_test(201,405)
+        pval = stats.binom_test(201, 405)
-        assert_almost_equal(pval,0.92085205962670713,11)
+        assert_almost_equal(pval, 0.92085205962670713, 11)
-        pval = stats.binom_test([682,243],p=3.0/4)
+        pval = stats.binom_test([682, 243], p=3.0 / 4)
-        assert_almost_equal(pval,0.38249155957481695,11)
+        assert_almost_equal(pval, 0.38249155957481695, 11)
 
     def test_bad_len_x(self):
         # Length of x must be 1 or 2.
-        assert_raises(ValueError, stats.binom_test, [1,2,3])
+        assert_raises(ValueError, stats.binom_test, [1, 2, 3])
 
     def test_bad_n(self):
         # len(x) is 1, but n is invalid.
@@ -211,10 +218,10 @@ class TestBinomP(TestCase):
 class TestFindRepeats(TestCase):
 
     def test_basic(self):
-        a = [1,2,3,4,1,2,3,4,1,2,5]
+        a = [1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 5]
-        res,nums = stats.find_repeats(a)
+        res, nums = stats.find_repeats(a)
-        assert_array_equal(res,[1,2,3,4])
+        assert_array_equal(res, [1, 2, 3, 4])
-        assert_array_equal(nums,[3,3,2,2])
+        assert_array_equal(nums, [3, 3, 2, 2])
 
     def test_empty_result(self):
         # Check that empty arrays are returned when there are no repeats.
@@ -229,14 +236,16 @@ class TestFligner(TestCase):
     def test_data(self):
         # numbers from R: fligner.test in package stats
         x1 = np.arange(5)
-        assert_array_almost_equal(stats.fligner(x1,x1**2),
+        assert_array_almost_equal(stats.fligner(x1, x1 ** 2),
-                           (3.2282229927203536, 0.072379187848207877), 11)
+                                  (3.2282229927203536, 0.072379187848207877),
+                                  11)
 
     def test_trimmed1(self):
         # Test that center='trimmed' gives the same result as center='mean'
         # when proportiontocut=0.
         Xsq1, pval1 = stats.fligner(g1, g2, g3, center='mean')
-        Xsq2, pval2 = stats.fligner(g1, g2, g3, center='trimmed', proportiontocut=0.0)
+        Xsq2, pval2 = stats.fligner(
+            g1, g2, g3, center='trimmed', proportiontocut=0.0)
         assert_almost_equal(Xsq1, Xsq2)
         assert_almost_equal(pval1, pval2)
 
@@ -244,7 +253,8 @@ class TestFligner(TestCase):
         x = [1.2, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 100.0]
         y = [0.0, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 200.0]
         # Use center='trimmed'
-        Xsq1, pval1 = stats.fligner(x, y, center='trimmed', proportiontocut=0.125)
+        Xsq1, pval1 = stats.fligner(
+            x, y, center='trimmed', proportiontocut=0.125)
         # Trim the data here, and use center='mean'
         Xsq2, pval2 = stats.fligner(x[1:-1], y[1:-1], center='mean')
         # Result should be the same.
@@ -257,7 +267,7 @@ class TestFligner(TestCase):
     # errors) there are not.  This difference leads to differences in the
     # third significant digit of W.
     #
-    #def test_equal_mean_median(self):
+    # def test_equal_mean_median(self):
     #    x = np.linspace(-1,1,21)
     #    y = x**3
     #    W1, pval1 = stats.fligner(x, y, center='mean')
@@ -266,11 +276,11 @@ class TestFligner(TestCase):
     #    assert_almost_equal(pval1, pval2)
 
     def test_bad_keyword(self):
-        x = np.linspace(-1,1,21)
+        x = np.linspace(-1, 1, 21)
         assert_raises(TypeError, stats.fligner, x, x, portiontocut=0.1)
 
     def test_bad_center_value(self):
-        x = np.linspace(-1,1,21)
+        x = np.linspace(-1, 1, 21)
         assert_raises(ValueError, stats.fligner, x, x, center='trim')
 
     def test_bad_num_args(self):
@@ -279,11 +289,13 @@ class TestFligner(TestCase):
 
 
 class TestMood(TestCase):
+
     def test_mood(self):
         # numbers from R: mood.test in package stats
         x1 = np.arange(5)
-        assert_array_almost_equal(stats.mood(x1, x1**2),
+        assert_array_almost_equal(stats.mood(x1, x1 ** 2),
-                                  (-1.3830857299399906, 0.16663858066771478), 11)
+                                  (-1.3830857299399906, 0.16663858066771478),
+                                  11)
 
     def test_mood_order_of_args(self):
         # z should change sign when the order of arguments changes, pvalue
@@ -296,7 +308,7 @@ class TestMood(TestCase):
         assert_array_almost_equal([z1, p1], [-z2, p2])
 
     def test_mood_with_axis_none(self):
-        #Test with axis = None, compare with results from R
+        # Test with axis = None, compare with results from R
         x1 = [-0.626453810742332, 0.183643324222082, -0.835628612410047,
               1.59528080213779, 0.329507771815361, -0.820468384118015,
               0.487429052428485, 0.738324705129217, 0.575781351653492,
@@ -375,7 +387,8 @@ class TestMood(TestCase):
                                               stats.mood(slice1, slice2))
 
     def test_mood_bad_arg(self):
-        # Raise ValueError when the sum of the lengths of the args is less than 3
+        # Raise ValueError when the sum of the lengths of the args is less than
+        # 3
         assert_raises(ValueError, stats.mood, [1], [])
 
 
@@ -393,7 +406,7 @@ class TestProbplot(TestCase):
         assert_allclose(osr, np.sort(x))
         assert_allclose(osm, osm_expected)
 
-        res, res_fit = stats.probplot(x, fit=True)
+        _res, res_fit = stats.probplot(x, fit=True)
         res_fit_expected = [1.05361841, 0.31297795, 0.98741609]
         assert_allclose(res_fit, res_fit_expected)
 
@@ -410,7 +423,7 @@ class TestProbplot(TestCase):
         assert_allclose(osr1, osr2)
         assert_allclose(osr1, osr3)
         # Check giving (loc, scale) params for normal distribution
-        osm, osr = stats.probplot(x, sparams=(), fit=False)
+        _osm, _osr = stats.probplot(x, sparams=(), fit=False)
 
     def test_dist_keyword(self):
         np.random.seed(12345)
@@ -424,7 +437,9 @@ class TestProbplot(TestCase):
         assert_raises(AttributeError, stats.probplot, x, dist=[])
 
         class custom_dist(object):
+
             """Some class that looks just enough like a distribution."""
+
             def ppf(self, q):
                 return stats.norm.ppf(q, loc=2)
 
@@ -467,8 +482,8 @@ class TestProbplot(TestCase):
 def test_wilcoxon_bad_arg():
     # Raise ValueError when two args of different lengths are given or
     # zero_method is unknown.
-    assert_raises(ValueError, stats.wilcoxon, [1], [1,2])
+    assert_raises(ValueError, stats.wilcoxon, [1], [1, 2])
-    assert_raises(ValueError, stats.wilcoxon, [1,2], [1,2], "dummy")
+    assert_raises(ValueError, stats.wilcoxon, [1, 2], [1, 2], "dummy")
 
 
 def test_mvsdist_bad_arg():
@@ -504,7 +519,7 @@ class TestBoxcox_llf(TestCase):
         x = stats.norm.rvs(size=10000, loc=10)
         lmbda = 1
         llf = stats.boxcox_llf(lmbda, x)
-        llf_expected = -x.size / 2. * np.log(np.sum(x.std()**2))
+        llf_expected = -x.size / 2. * np.log(np.sum(x.std() ** 2))
         assert_allclose(llf, llf_expected)
 
     def test_array_like(self):
@@ -538,7 +553,7 @@ class TestBoxcox(TestCase):
         xt = stats.boxcox(x, lmbda=1)
         assert_allclose(xt, x - 1)
         xt = stats.boxcox(x, lmbda=-1)
-        assert_allclose(xt, 1 - 1/x)
+        assert_allclose(xt, 1 - 1 / x)
 
         xt = stats.boxcox(x, lmbda=0)
         assert_allclose(xt, np.log(x))
@@ -554,8 +569,8 @@ class TestBoxcox(TestCase):
         np.random.seed(1245)
         lmbda = 2.5
         x = stats.norm.rvs(loc=10, size=50000)
-        x_inv = (x * lmbda + 1)**(-lmbda)
+        x_inv = (x * lmbda + 1) ** (-lmbda)
-        xt, maxlog = stats.boxcox(x_inv)
+        _xt, maxlog = stats.boxcox(x_inv)
 
         assert_almost_equal(maxlog, -1 / lmbda, decimal=2)
 
@@ -586,6 +601,7 @@ class TestBoxcox(TestCase):
 
 
 class TestBoxcoxNormmax(TestCase):
+
     def setUp(self):
         np.random.seed(12345)
         self.x = stats.loggamma.rvs(5, size=50) + 5
@@ -608,6 +624,7 @@ class TestBoxcoxNormmax(TestCase):
 
 
 class TestBoxcoxNormplot(TestCase):
+
     def setUp(self):
         np.random.seed(7654321)
         self.x = stats.loggamma.rvs(5, size=500) + 5
@@ -638,15 +655,16 @@ class TestBoxcoxNormplot(TestCase):
         # `lb` has to be larger than `la`
         assert_raises(ValueError, stats.boxcox_normplot, self.x, 1, 0)
         # `x` can not contain negative values
-        assert_raises(ValueError, stats.boxcox_normplot, [-1, 1] , 0, 1)
+        assert_raises(ValueError, stats.boxcox_normplot, [-1, 1], 0, 1)
 
     def test_empty(self):
         assert_(stats.boxcox_normplot([], 0, 1).size == 0)
 
 
 class TestCircFuncs(TestCase):
+
     def test_circfuncs(self):
-        x = np.array([355,5,2,359,10,350])
+        x = np.array([355, 5, 2, 359, 10, 350])
         M = stats.circmean(x, high=360)
         Mval = 0.167690146
         assert_allclose(M, Mval, rtol=1e-7)
@@ -660,7 +678,7 @@ class TestCircFuncs(TestCase):
         assert_allclose(S, Sval, rtol=1e-7)
 
     def test_circfuncs_small(self):
-        x = np.array([20,21,22,18,19,20.5,19.2])
+        x = np.array([20, 21, 22, 18, 19, 20.5, 19.2])
         M1 = x.mean()
         M2 = stats.circmean(x, high=360)
         assert_allclose(M2, M1, rtol=1e-5)
@@ -674,9 +692,9 @@ class TestCircFuncs(TestCase):
         assert_allclose(S2, S1, rtol=1e-4)
 
     def test_circmean_axis(self):
-        x = np.array([[355,5,2,359,10,350],
+        x = np.array([[355, 5, 2, 359, 10, 350],
-                      [351,7,4,352,9,349],
+                      [351, 7, 4, 352, 9, 349],
-                      [357,9,8,358,4,356]])
+                      [357, 9, 8, 358, 4, 356]])
         M1 = stats.circmean(x, high=360)
         M2 = stats.circmean(x.ravel(), high=360)
         assert_allclose(M1, M2, rtol=1e-14)
@@ -686,13 +704,13 @@ class TestCircFuncs(TestCase):
         assert_allclose(M1, M2, rtol=1e-14)
 
         M1 = stats.circmean(x, high=360, axis=0)
-        M2 = [stats.circmean(x[:,i], high=360) for i in range(x.shape[1])]
+        M2 = [stats.circmean(x[:, i], high=360) for i in range(x.shape[1])]
         assert_allclose(M1, M2, rtol=1e-14)
 
     def test_circvar_axis(self):
-        x = np.array([[355,5,2,359,10,350],
+        x = np.array([[355, 5, 2, 359, 10, 350],
-                      [351,7,4,352,9,349],
+                      [351, 7, 4, 352, 9, 349],
-                  [357,9,8,358,4,356]])
+                      [357, 9, 8, 358, 4, 356]])
 
         V1 = stats.circvar(x, high=360)
         V2 = stats.circvar(x.ravel(), high=360)
@@ -703,13 +721,13 @@ class TestCircFuncs(TestCase):
         assert_allclose(V1, V2, rtol=1e-11)
 
         V1 = stats.circvar(x, high=360, axis=0)
-        V2 = [stats.circvar(x[:,i], high=360) for i in range(x.shape[1])]
+        V2 = [stats.circvar(x[:, i], high=360) for i in range(x.shape[1])]
         assert_allclose(V1, V2, rtol=1e-11)
 
     def test_circstd_axis(self):
-        x = np.array([[355,5,2,359,10,350],
+        x = np.array([[355, 5, 2, 359, 10, 350],
-                      [351,7,4,352,9,349],
+                      [351, 7, 4, 352, 9, 349],
-                      [357,9,8,358,4,356]])
+                      [357, 9, 8, 358, 4, 356]])
 
         S1 = stats.circstd(x, high=360)
         S2 = stats.circstd(x.ravel(), high=360)
@@ -720,11 +738,11 @@ class TestCircFuncs(TestCase):
         assert_allclose(S1, S2, rtol=1e-11)
 
         S1 = stats.circstd(x, high=360, axis=0)
-        S2 = [stats.circstd(x[:,i], high=360) for i in range(x.shape[1])]
+        S2 = [stats.circstd(x[:, i], high=360) for i in range(x.shape[1])]
         assert_allclose(S1, S2, rtol=1e-11)
 
     def test_circfuncs_array_like(self):
-        x = [355,5,2,359,10,350]
+        x = [355, 5, 2, 359, 10, 350]
         assert_allclose(stats.circmean(x, high=360), 0.167690146, rtol=1e-7)
         assert_allclose(stats.circvar(x, high=360), 42.51955609, rtol=1e-7)
         assert_allclose(stats.circstd(x, high=360), 6.520702116, rtol=1e-7)

pywafo/src/wafo/stats/tests/test_multivariate.py

@@ -11,11 +11,11 @@ import numpy
 import numpy as np
 
 import scipy.linalg
-import scipy.stats._multivariate
+#import wafo.stats._multivariate
-from scipy.stats import multivariate_normal
+from wafo.stats import multivariate_normal
-from scipy.stats import norm
+from wafo.stats import norm
 
-from scipy.stats._multivariate import _psd_pinv_decomposed_log_pdet
+from wafo.stats._multivariate import _psd_pinv_decomposed_log_pdet
 
 from scipy.integrate import romb
 
@@ -70,7 +70,7 @@ def test_large_pseudo_determinant():
     #assert_allclose(np.linalg.slogdet(cov[:npos, :npos]), (1, large_total_log))
 
     # Check the pseudo-determinant.
-    U, log_pdet = scipy.stats._multivariate._psd_pinv_decomposed_log_pdet(cov)
+    U, log_pdet = _psd_pinv_decomposed_log_pdet(cov)
     assert_allclose(log_pdet, large_total_log)
 
 

pywafo/src/wafo/stats/tests/test_rank.py

@@ -4,7 +4,7 @@ import numpy as np
 from numpy.testing import TestCase, run_module_suite, assert_equal, \
     assert_array_equal
 
-from scipy.stats import rankdata, tiecorrect
+from wafo.stats import rankdata, tiecorrect
 
 
 class TestTieCorrect(TestCase):

pywafo/src/wafo/stats/tests/test_tukeylambda_stats.py

@@ -18,7 +18,7 @@ def test_tukeylambda_stats_known_exact():
 
     # lambda = 0
     var = tukeylambda_variance(0)
-    assert_allclose(var, np.pi**2 / 3, atol=1e-12)
+    assert_allclose(var, np.pi ** 2 / 3, atol=1e-12)
     kurt = tukeylambda_kurtosis(0)
     assert_allclose(kurt, 1.2, atol=1e-10)
 
@@ -26,7 +26,7 @@ def test_tukeylambda_stats_known_exact():
     var = tukeylambda_variance(0.5)
     assert_allclose(var, 4 - np.pi, atol=1e-12)
     kurt = tukeylambda_kurtosis(0.5)
-    desired = (5./3 - np.pi/2) / (np.pi/4 - 1)**2 - 3
+    desired = (5. / 3 - np.pi / 2) / (np.pi / 4 - 1) ** 2 - 3
     assert_allclose(kurt, desired, atol=1e-10)
 
     # lambda = 1

pywafo/src/wafo/test/test_misc.py

@@ -1,51 +1,44 @@
-import numpy as np  # @UnusedImport
+
-#@UnusedImport
+from numpy.testing import (run_module_suite, assert_equal, assert_almost_equal,
-from numpy import cos, exp, linspace, pi, sin, diff, arange, ones
+                           assert_array_equal, assert_array_almost_equal)
-from numpy.random import randn  # @UnusedImport
+
-from wafo.data import sea  # @UnusedImport
+import numpy as np
-from wafo.misc import (JITImport, Bunch, detrendma, DotDict, findcross, ecross, findextrema,  # @UnusedImport
+from numpy import array, cos, exp, linspace, pi, sin, diff, arange, ones
-                       #@UnusedImport
+from wafo.data import sea
-                       findrfc, rfcfilter, findtp, findtc, findoutliers,
+from wafo.misc import (JITImport, Bunch, detrendma, DotDict, findcross, ecross,
-                       common_shape, argsreduce, stirlerr, getshipchar, betaloge,
+                       findextrema, findrfc, rfcfilter, findtp, findtc,
-                       #@UnusedImport
+                       findoutliers, common_shape, argsreduce, stirlerr,
-                       #@UnusedImport
+                       getshipchar, betaloge, hygfz,
                        gravity, nextpow2, discretize, polar2cart,
-                       cart2polar, meshgrid, tranproc)  # @UnusedImport
+                       cart2polar, tranproc)
 
 
 def test_JITImport():
-    '''
+    np = JITImport('numpy')
-    >>> np = JITImport('numpy')
+    assert_equal(1.0, np.exp(0))
-    >>> np.exp(0)==1.0
-    True
-    '''
 
 
 def test_bunch():
-    '''
+    d = Bunch(test1=1, test2=3)
-    >>> d = Bunch(test1=1,test2=3)
+    assert_equal(1, d.test1)
-    >>> d.test1; d.test2
+    assert_equal(3, d.test2)
-    1
-    3
-    '''
 
 
 def test_dotdict():
-    '''
+    d = DotDict(test1=1, test2=3)
-    >>> d = DotDict(test1=1,test2=3)
+    assert_equal(1, d.test1)
-    >>> d.test1; d.test2
+    assert_equal(3, d.test2)
-    1
-    3
-    '''
 
 
 def test_detrendma():
-    '''
+    x = linspace(0, 1, 200)
-    >>> x = linspace(0,1,200)
+    y = exp(x) + 0.1 * cos(20 * 2 * pi * x)
-    >>> y = exp(x)+0.1*cos(20*2*pi*x)
+    y0 = detrendma(y, 20)
-    >>> y0 = detrendma(y,20); tr = y-y0
+    tr = y - y0
-    >>> y0,tr
+    assert_array_almost_equal(
-    (array([ -1.05815186e-02,  -2.48280355e-02,  -7.01800760e-02,
+        y0,
+        array(
+            [-1.05815186e-02, -2.48280355e-02, -7.01800760e-02,
              -1.27193089e-01, -1.71915213e-01, -1.85125121e-01,
              -1.59745361e-01, -1.03571981e-01, -3.62676515e-02,
              1.82219951e-02, 4.09039083e-02, 2.50630186e-02,
@@ -111,7 +104,9 @@ def test_detrendma():
              2.43802139e-01, 2.39414013e-01, 2.03257341e-01,
              1.54325635e-01, 1.16564992e-01, 1.09638547e-01,
              1.41342814e-01, 2.04600808e-01, 2.80191671e-01,
-             3.44164010e-01,   3.77073744e-01]), array([
+             3.44164010e-01, 3.77073744e-01
+             ]))
+    assert_array_almost_equal(tr, array([
             1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152,
             1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152,
             1.11058152, 1.11058152, 1.11058152, 1.11058152, 1.11058152,
@@ -123,74 +118,66 @@ def test_detrendma():
             1.23051218, 1.23633498, 1.24209697, 1.24791509, 1.25389641,
             1.26009689, 1.26649987, 1.27302256, 1.27954802, 1.28597031,
             1.29223546, 1.29836228, 1.30443522, 1.31057183, 1.31687751,
-            1.32340488,  1.3301336 ,  1.33697825,  1.34382132,  1.35055864,
+            1.32340488, 1.3301336, 1.33697825, 1.34382132, 1.35055864,
             1.35713958, 1.36358668, 1.36998697, 1.37645853, 1.38310497,
             1.38997553, 1.39704621, 1.40422902, 1.41140604, 1.41847493,
-            1.4253885 ,  1.43217295,  1.43891784,  1.44574164,  1.45274607,
+            1.4253885, 1.43217295, 1.43891784, 1.44574164, 1.45274607,
             1.45997696, 1.46740665, 1.47494469, 1.48247285, 1.48989073,
             1.49715462, 1.50429437, 1.51140198, 1.51859618, 1.52597672,
-            1.53358594,  1.54139257,  1.5493038 ,  1.55720119,  1.56498641,
+            1.53358594, 1.54139257, 1.5493038, 1.55720119, 1.56498641,
-            1.57261924,  1.58013316,  1.58762252,  1.5952062 ,  1.60298187,
+            1.57261924, 1.58013316, 1.58762252, 1.5952062, 1.60298187,
-            1.61098836,  1.6191908 ,  1.62749412,  1.63577979,  1.64395163,
+            1.61098836, 1.6191908, 1.62749412, 1.63577979, 1.64395163,
             1.65197298, 1.65988092, 1.66777202, 1.67576523, 1.68395602,
             1.69237968, 1.70099778, 1.70971307, 1.71840707, 1.72698583,
             1.73541631, 1.74373911, 1.75205298, 1.76047677, 1.76910369,
             1.77796544, 1.78702008, 1.79616827, 1.80529169, 1.81429875,
-            1.82316   ,  1.83191959,  1.84067831,  1.84955481,  1.85863994,
+            1.82316, 1.83191959, 1.84067831, 1.84955481, 1.85863994,
-            1.86796178,  1.87747491,  1.88707803,  1.89665308,  1.9061109 ,
+            1.86796178, 1.87747491, 1.88707803, 1.89665308, 1.9061109,
-            1.91542572,  1.92464514,  1.9338719 ,  1.94322436,  1.9527909 ,
+            1.91542572, 1.92464514, 1.9338719, 1.94322436, 1.9527909,
-            1.96259596,  1.97259069,  1.9826719 ,  1.99272195,  2.00265419,
+            1.96259596, 1.97259069, 1.9826719, 1.99272195, 2.00265419,
-            2.01244653,  2.02215   ,  2.0318692 ,  2.04172204,  2.05179437,
+            2.01244653, 2.02215, 2.0318692, 2.04172204, 2.05179437,
             2.06210696, 2.07260759, 2.08319129, 2.09374092, 2.10417247,
-            2.11446752,  2.12468051,  2.13491776,  2.14529665,  2.1559004 ,
+            2.11446752, 2.12468051, 2.13491776, 2.14529665, 2.1559004,
             2.16674609, 2.17777817, 2.18889002, 2.19996511, 2.21092214,
             2.22174641, 2.23249567, 2.24327791, 2.25420982, 2.26537192,
-            2.2767776 ,  2.28836802,  2.30003501,  2.3116628 ,  2.32317284,
+            2.2767776, 2.28836802, 2.30003501, 2.3116628, 2.32317284,
             2.33455419, 2.34586786, 2.35722337, 2.36873665, 2.38048542,
-            2.39247934,  2.4046564 ,  2.41690694,  2.42911606,  2.44120808,
+            2.39247934, 2.4046564, 2.41690694, 2.42911606, 2.44120808,
             2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808,
             2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808,
             2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808,
             2.44120808, 2.44120808, 2.44120808, 2.44120808, 2.44120808]))
-    '''
 
 
 def test_findcross_and_ecross():
-    '''
+    assert_array_equal(findcross([0, 0, 1, -1, 1], 0), np.array([1, 2, 3]))
-    >>> findcross([0, 0, 1, -1, 1],0)
+    assert_array_equal(findcross([0, 1, -1, 1], 0), np.array([0, 1, 2]))
-    array([1, 2, 3])
+
-    >>> findcross([0, 1, -1, 1],0)
+    t = linspace(0, 7 * pi, 250)
-    array([0, 1, 2])
+    x = sin(t)
-
+    ind = findcross(x, 0.75)
-    >>> t = linspace(0,7*pi,250)
+    assert_array_equal(ind, np.array([9, 25, 80, 97, 151, 168, 223, 239]))
-    >>> x = sin(t)
+    t0 = ecross(t, x, ind, 0.75)
-    >>> ind = findcross(x,0.75)
+    assert_array_almost_equal(t0, np.array([0.84910514, 2.2933879, 7.13205663,
-    >>> ind
+                                     8.57630119, 13.41484739, 14.85909194,
-    array([  9,  25,  80,  97, 151, 168, 223, 239])
+                                     19.69776067, 21.14204343]))
-    >>> t0 = ecross(t,x,ind,0.75)
-    >>> t0
-    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)
-    >>> t = linspace(0,7*pi,250)
+    x = sin(t)
-    >>> x = sin(t)
+    ind = findextrema(x)
-    >>> ind = findextrema(x)
+    assert_array_almost_equal(ind, np.array([18, 53, 89, 125, 160, 196, 231]))
-    >>> ind
-    array([ 18,  53,  89, 125, 160, 196, 231])
-    '''
 
 
 def test_findrfc():
-    '''
+    t = linspace(0, 7 * pi, 250)
-    >>> t = linspace(0,7*pi,250)
+    x = sin(t) + 0.1 * sin(50 * t)
-    >>> x = sin(t)+0.1*sin(50*t)
+    ind = findextrema(x)
-    >>> ind = findextrema(x)
+    assert_array_almost_equal(
-    >>> ind
+        ind,
-    array([  1,   3,   4,   6,   7,   9,  11,  13,  14,  16,  18,  19,  21,
+        np.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,
@@ -201,35 +188,42 @@ def test_findrfc():
              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])
+             241, 243, 245, 247, 248]))
-    >>> ti, tp = t[ind], x[ind]
+    _ti, tp = t[ind], x[ind]
-    >>> ind1 = findrfc(tp,0.3)
+    ind1 = findrfc(tp, 0.3)
-    >>> ind1
+    assert_array_almost_equal(
-    array([  0,   9,  32,  53,  74,  95, 116, 137])
+        ind1,
-    >>> tp[ind1]
+        np.array([0,  9, 32, 53, 74, 95, 116, 137]))
-    array([-0.00743352,  1.08753972, -1.07206545,  1.09550837, -1.07940458,
+    assert_array_almost_equal(
-            1.07849396, -1.0995006 ,  1.08094452])
+        tp[ind1],
-    '''
+        np.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()
+    x = sea()
-    >>> y = rfcfilter(x[:,1], h=0, method=1)
+    y = rfcfilter(x[:, 1], h=0, method=1)
-    >>> y[0:5]
+    assert_array_almost_equal(
-    array([-1.2004945 ,  0.83950546, -0.09049454, -0.02049454, -0.09049454])
+        y[0:5],
+        np.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 = rfcfilter(x[:, 1], h=0.5)
-    >>> y1[0:5]
+    assert_array_almost_equal(
-    array([-1.2004945 ,  0.83950546, -0.43049454,  0.34950546, -0.51049454])
+        y1[0:5],
-
+        np.array([-1.2004945, 0.83950546, -0.43049454,
-    >>> t = linspace(0,7*pi,250)
+                  0.34950546, -0.51049454]))
-    >>> x = sin(t)+0.1*sin(50*t)
+
-    >>> ind = findextrema(x)
+    t = linspace(0, 7 * pi, 250)
-    >>> ind
+    x = sin(t) + 0.1 * sin(50 * t)
-    array([  1,   3,   4,   6,   7,   9,  11,  13,  14,  16,  18,  19,  21,
+    ind = findextrema(x)
+    assert_array_almost_equal(
+        ind,
+        np.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,
@@ -240,284 +234,237 @@ def test_rfcfilter():
              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])
+             241, 243, 245, 247, 248]))
-    >>> ti, tp = t[ind], x[ind]
+    _ti, tp = t[ind], x[ind]
-    >>> tp03 = rfcfilter(tp,0.3)
+    tp03 = rfcfilter(tp, 0.3)
-    >>> tp03
+    assert_array_almost_equal(
-    array([-0.00743352,  1.08753972, -1.07206545,  1.09550837, -1.07940458,
+        tp03,
-            1.07849396, -1.0995006 ,  1.08094452,  0.11983423])
+        np.array(
-    '''
+            [-0.00743352, 1.08753972, -1.07206545, 1.09550837, -1.07940458,
+             1.07849396, -1.0995006, 1.08094452, 0.11983423]))
 
 
 def test_findtp():
-    '''
+    x = sea()
-    >>> import numpy as np
+    x1 = x[0:200, :]
-    >>> x = sea()
+    itp = findtp(x1[:, 1], 0, 'Mw')
-    >>> x1 = x[0:200,:]
+    itph = findtp(x1[:, 1], 0.3, 'Mw')
-    >>> itp = findtp(x1[:,1],0,'Mw')
+    assert_array_almost_equal(
-    >>> itph = findtp(x1[:,1],0.3,'Mw')
+        itp,
-    >>> itp
+        np.array(
-    array([ 11,  21,  22,  24,  26,  28,  31,  39,  43,  45,  47,  51,  56,
+            [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])
+             149, 150, 159, 173, 184, 190, 199]))
-    >>> itph
+    assert_array_almost_equal(
-    array([ 11,  28,  31,  39,  47,  51,  56,  64,  70,  78,  89,  94, 101,
+        itph,
-           108, 119, 131, 141, 148, 159, 173, 184, 190, 199])
+        np.array(
-    '''
+            [11, 28, 31, 39, 47, 51, 56, 64, 70, 78, 89, 94, 101,
+             108, 119, 131, 141, 148, 159, 173, 184, 190, 199]))
 
 
 def test_findtc():
-    '''
+    x = sea()
-    >>> x = sea()
+    x1 = x[0:200, :]
-    >>> x1 = x[0:200,:]
+    itc, iv = findtc(x1[:, 1], 0, 'dw')
-    >>> itc, iv = findtc(x1[:,1],0,'dw')
+    assert_array_almost_equal(
-    >>> itc
+        itc,
-    array([ 28,  31,  39,  56,  64,  69,  78,  82,  83,  89,  94, 101, 108,
+        np.array(
-           119, 131, 140, 148, 159, 173, 184])
+            [28, 31, 39, 56, 64, 69, 78, 82, 83, 89, 94, 101, 108,
-    >>> iv
+             119, 131, 140, 148, 159, 173, 184]))
-    array([ 19,  29,  34,  53,  60,  67,  76,  81,  82,  84,  90,  99, 103,
+    assert_array_almost_equal(
-           112, 127, 137, 143, 154, 166, 180, 185])
+        iv,
-    '''
+        np.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()
-    >>> xx = sea()
+    dt = diff(xx[:2, 0])
-    >>> dt = diff(xx[:2,0])
+    dcrit = 5 * dt
-    >>> dcrit = 5*dt
+    ddcrit = 9.81 / 2 * dt * dt
-    >>> ddcrit = 9.81/2*dt*dt
+    zcrit = 0
-    >>> zcrit = 0
+    [inds, indg] = findoutliers(xx[:, 1], zcrit, dcrit, ddcrit, verbose=False)
-    >>> [inds, indg] = findoutliers(xx[:,1],zcrit,dcrit,ddcrit,verbose=True)
+    assert_array_almost_equal(inds[np.r_[0, 1, 2, -3, -2, -1]],
-    Found 0 spurious positive jumps of Dx
+                              np.array([6,   7,   8, 9509, 9510, 9511]))
-    Found 0 spurious negative jumps of Dx
+    assert_array_almost_equal(indg[np.r_[0, 1, 2, -3, -2, -1]],
-    Found 37 spurious positive jumps of D^2x
+                              np.array([0,   1,   2, 9521, 9522, 9523]))
-    Found 200 spurious negative jumps of D^2x
+
-    Found 244 consecutive equal values
+
-    Found the total of 1152 spurious points
+def test_hygfz():
-    >>> inds
+    #y = hyp2f1_taylor(-1, -4, 1, .9)
-    array([   6,    7,    8, ..., 9509, 9510, 9511])
+    assert_equal(4.6, hygfz(-1, -4, 1, .9))
-    >>> indg
+    assert_almost_equal(1.0464328112173522, hygfz(0.1, 0.2, 0.3, 0.5))
-    array([   0,    1,    2, ..., 9521, 9522, 9523])
+    assert_almost_equal(1.2027034401166194, hygfz(0.1, 0.2, 0.3, 0.95))
-    '''
+    #assert_equal(1.661006238211309e-07, hygfz(5, -300, 10, 0.5))
-
+    assert_equal(0.118311386286, hygfz(0.5, -99.0, 1.5, 0.5625))
+    assert_equal(0.0965606007742, hygfz(0.5, -149.0, 1.5, 0.5625))
+    assert_equal(0.49234384000963544+0.60513406166123973j, hygfz(1, 1, 4, 3+4j))
 
 def test_common_shape():
-    '''
+    A = np.ones((4, 1))
-    >>> import numpy as np
+    B = 2
-    >>> A = np.ones((4,1))
+    C = np.ones((1, 5)) * 5
-    >>> B = 2
+    assert_array_equal(common_shape(A, B, C), (4, 5))
-    >>> C = np.ones((1,5))*5
+    assert_array_equal(common_shape(A, B, C, shape=(3, 4, 1)), (3, 4, 5))
-    >>> common_shape(A,B,C)
+    A = np.ones((4, 1))
-    (4, 5)
+    B = 2
-    >>> common_shape(A,B,C,shape=(3,4,1))
+    C = np.ones((1, 5)) * 5
-    (3, 4, 5)
+    assert_array_equal(common_shape(A, B, C), (4, 5))
-    >>> A = np.ones((4,1))
+    assert_array_equal(common_shape(A, B, C, shape=(3, 4, 1)), (3, 4, 5))
-    >>> 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():
-    '''
+    A = linspace(0, 19, 20).reshape((4, 5))
-    >>> import numpy as np
+    B = 2
-    >>> rand = np.random.random_sample
+    C = range(5)
-    >>> A = linspace(0,19,20).reshape((4,5))
+    cond = np.ones(A.shape)
-    >>> B = 2
+    [_A1, B1, _C1] = argsreduce(cond, A, B, C)
-    >>> C = range(5)
+    assert_equal(B1.shape, (20,))
-    >>> cond = np.ones(A.shape)
+    cond[2, :] = 0
-    >>> [A1,B1,C1] = argsreduce(cond,A,B,C)
+    [A2, B2, C2] = argsreduce(cond, A, B, C)
-    >>> B1.shape
+    assert_equal(B2.shape, (15,))
-    (20,)
+    assert_array_equal(A2,
-    >>> cond[2,:] = 0
+                       np.array([0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,
-    >>> [A2,B2,C2] = argsreduce(cond,A,B,C)
+                                8.,  9., 15., 16., 17., 18., 19.]))
-    >>> B2.shape
+    assert_array_equal(
-    (15,)
+        B2, np.array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]))
-    >>> A2;B2;C2
+    assert_array_equal(
-    array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  15.,
+        C2, np.array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4]))
-            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():
-    '''
+    assert_array_almost_equal(stirlerr(range(5)),
-     >>> stirlerr(range(5))
+                       np.array([np.inf, 0.08106147, 0.0413407, 0.02767793,
-     array([        inf,  0.08106147,  0.0413407 ,  0.02767793,  0.02079067])
+                                 0.02079067]))
-    '''
 
 
 def test_getshipchar():
-    '''
+    sc = getshipchar(10, 'service_speed')
-    >>> sc = getshipchar(10,'service_speed')
+    true_sc = dict(beam=29,
-    >>> names = ['beam',    'beamSTD',     'draught',
+                   beamSTD=2.9,
-    ... 'draughtSTD',     'length',     'lengthSTD',
+                   draught=9.6,
-    ... 'max_deadweight',     'max_deadweightSTD',     'propeller_diameter',
+                   draughtSTD=2.112,
-    ... 'propeller_diameterSTD',     'service_speed',     'service_speedSTD']
+                   length=216,
-    >>> for name in names: print( '%s : %g' % (name, sc[name]))
+                   lengthSTD=2.011309883194276,
-    beam : 29
+                   max_deadweight=30969,
-    beamSTD : 2.9
+                   max_deadweightSTD=3096.9,
-    draught : 9.6
+                   propeller_diameter=6.761165385916601,
-    draughtSTD : 2.112
+                   propeller_diameterSTD=0.20267047566705432,
-    length : 216
+                   service_speed=10,
-    lengthSTD : 2.01131
+                   service_speedSTD=0)
-    max_deadweight : 30969
+
-    max_deadweightSTD : 3096.9
+    for name, val in true_sc.iteritems():
-    propeller_diameter : 6.76117
+        assert_almost_equal(val, sc[name])
-    propeller_diameterSTD : 0.20267
-    service_speed : 10
-    service_speedSTD : 0
-    '''
 
 
 def test_betaloge():
-    '''
+    assert_array_almost_equal(betaloge(3, arange(4)),
-    >>> betaloge(3, arange(4))
+    np.array([np.inf, -1.09861229, -2.48490665, -3.40119738]))
-    array([        inf, -1.09861229, -2.48490665, -3.40119738])
-    '''
 
 
 def test_gravity():
-    '''
+    phi = linspace(0, 45, 5)
-    >>> phi = linspace(0,45,5)
+    assert_array_almost_equal(gravity(phi),
-    >>> gravity(phi)
+                              np.array([9.78049, 9.78245014, 9.78803583,
-    array([ 9.78049   ,  9.78245014,  9.78803583,  9.79640552,  9.80629387])
+                                        9.79640552, 9.80629387]))
-    '''
 
 
 def test_nextpow2():
-    '''
+    assert_equal(nextpow2(10), 4)
-    >>> nextpow2(10)
+    assert_equal(nextpow2(np.arange(5)), 3)
-    4
-    >>> nextpow2(np.arange(5))
-    3
-    '''
 
 
 def test_discretize():
-    '''
+    x, y = discretize(np.cos, 0, np.pi)
-    >>> x, y = discretize(np.cos,0,np.pi)
+    assert_array_almost_equal(
-    >>> x; y
+        x,
-    array([ 0.        ,  0.19634954,  0.39269908,  0.58904862,  0.78539816,
+        np.array(
-            0.9817477 ,  1.17809725,  1.37444679,  1.57079633,  1.76714587,
+            [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])
+             2.94524311, 3.14159265]))
-    array([  1.00000000e+00,   9.80785280e-01,   9.23879533e-01,
+    assert_array_almost_equal(
+        y, np.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])
+                     -9.80785280e-01, -1.00000000e+00]))
-    '''
 
 
 def test_discretize_adaptive():
-    '''
+    x, y = discretize(np.cos, 0, np.pi, method='adaptive')
-    >>> x, y = discretize(np.cos,0,np.pi, method='adaptive')
+    assert_array_almost_equal(
-    >>> x; y
+        x,
-    array([ 0.        ,  0.19634954,  0.39269908,  0.58904862,  0.78539816,
+        np.array(
-            0.9817477 ,  1.17809725,  1.37444679,  1.57079633,  1.76714587,
+            [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])
+             2.94524311, 3.14159265]))
-    array([  1.00000000e+00,   9.80785280e-01,   9.23879533e-01,
+    assert_array_almost_equal(
+        y,
+        np.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])
+             -9.80785280e-01, -1.00000000e+00]))
-    '''
 
 
-def test_pol2cart_n_cart2pol():
+def test_polar2cart_n_cart2polar():
-    '''
+    r = 5
-    >>> r = 5
+    t = linspace(0, pi, 20)
-    >>> t = linspace(0,pi,20)
+    x, y = polar2cart(t, r)
-    >>> x, y = polar2cart(t,r)
+    assert_array_almost_equal(
-    >>> x; y
+        x,
-    array([ 5.        ,  4.93180652,  4.72908621,  4.39736876,  3.94570255,
+        np.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.        ])
+             -3.94570255, -4.39736876, -4.72908621, -4.93180652, -5.]))
-    array([  0.00000000e+00,   8.22972951e-01,   1.62349735e+00,
+    assert_array_almost_equal(
+        y,
+        np.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])
+             8.22972951e-01,  6.12323400e-16]))
-    >>> ti, ri = cart2polar(x,y)
+    ti, ri = cart2polar(x, y)
-    >>> ti;ri
+    assert_array_almost_equal(
-    array([ 0.        ,  0.16534698,  0.33069396,  0.49604095,  0.66138793,
+        ti,
+        np.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,
+             1.65346982, 1.8188168, 1.98416378, 2.14951076, 2.31485774,
-            2.48020473,  2.64555171,  2.81089869,  2.97624567,  3.14159265])
+             2.48020473, 2.64555171, 2.81089869, 2.97624567, 3.14159265]))
-    array([ 5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,
+    assert_array_almost_equal(
-            5.,  5.,  5.,  5.,  5.,  5.,  5.])
+        ri,
-    '''
+        np.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
-    >>> import wafo.transform.models as wtm
+    tr = wtm.TrHermite()
-    >>> tr = wtm.TrHermite()
+    x = linspace(-5, 5, 501)
-    >>> x = linspace(-5,5,501)
+    g = tr(x)
-    >>> g = tr(x)
+    y0, y1 = tranproc(x, g, range(5), ones(5))
-    >>> y0, y1 = tranproc(x, g, range(5), ones(5))
+    assert_array_almost_equal(
-    >>> y0;y1
+        y0,
-    array([ 0.02659612,  1.00115284,  1.92872532,  2.81453257,  3.66292878])
+        np.array([0.02659612, 1.00115284, 1.92872532,
-    array([ 1.00005295,  0.9501118 ,  0.90589954,  0.86643821,  0.83096482])
+                  2.81453257, 3.66292878]))
-    '''
+    assert_array_almost_equal(
+        y1,
+        np.array([1.00005295, 0.9501118, 0.90589954,
+                  0.86643821, 0.83096482]))
+
+
 if __name__ == '__main__':
-    import doctest
+    run_module_suite()
-    doctest.testmod()