updated doctests + fitting

wafo/misc.py

@@ -302,20 +302,20 @@ def rotation_matrix(heading, pitch, roll):
            [ 0.,  1.,  0.],
            [ 0.,  0.,  1.]])
 
-    >>> np.all(np.abs(rotation_matrix(heading=180, pitch=0, roll=0)-
-    ... np.array([[ -1.000000e+00,  -1.224647e-16,   0.000000e+00],
-    ...       [  1.224647e-16,  -1.000000e+00,   0.000000e+00],
-    ...       [ -0.000000e+00,   0.000000e+00,   1.000000e+00]]))<1e-7)
+    >>> np.allclose(rotation_matrix(heading=180, pitch=0, roll=0),
+    ...      [[ -1.,   0.,   0.],
+    ...       [  0.,  -1.,   0.],
+    ...       [  0.,   0.,   1.]])
     True
-    >>> np.all(np.abs(rotation_matrix(heading=0, pitch=180, roll=0)-
-    ... np.array([[ -1.000000e+00,   0.000000e+00,   1.224647e-16],
-    ...       [ -0.000000e+00,   1.000000e+00,   0.000000e+00],
-    ...       [ -1.224647e-16,  -0.000000e+00,  -1.000000e+00]]))<1e-7)
+    >>> np.allclose(rotation_matrix(heading=0, pitch=180, roll=0),
+    ...      [[ -1.,  0.,   0.],
+    ...       [  0.,  1.,   0.],
+    ...       [  0.,  0.,  -1.]])
     True
-    >>> np.all(np.abs(rotation_matrix(heading=0, pitch=0, roll=180)-
-    ... np.array([[  1.000000e+00,   0.000000e+00,   0.000000e+00],
-    ...       [  0.000000e+00,  -1.000000e+00,  -1.224647e-16],
-    ...       [ -0.000000e+00,   1.224647e-16,  -1.000000e+00]]))<1e-7)
+    >>> np.allclose(rotation_matrix(heading=0, pitch=0, roll=180),
+    ...      [[  1.,  0.,   0.],
+    ...       [ 0.,  -1.,   0.],
+    ...       [ 0.,   0.,  -1.]])
     True
     '''
     data = np.diag(np.ones(3))  # No transform if H=P=R=0
@@ -369,14 +369,14 @@ def rotate_2d(x, y, angle_deg):
 
     Examples
     --------
-    >>> rotate_2d(x=1, y=0, angle_deg=0)
-    (1.0, 0.0)
-    >>> rotate_2d(x=1, y=0, angle_deg=90)
-    (6.123233995736766e-17, 1.0)
-    >>> rotate_2d(x=1, y=0, angle_deg=180)
-    (-1.0, 1.2246467991473532e-16)
-    >>> rotate_2d(x=1, y=0, angle_deg=360)
-    (1.0, -2.4492935982947064e-16)
+    >>> np.allclose(rotate_2d(x=1, y=0, angle_deg=0), (1.0, 0.0))
+    True
+    >>> np.allclose(rotate_2d(x=1, y=0, angle_deg=90), (0, 1.0))
+    True
+    >>> np.allclose(rotate_2d(x=1, y=0, angle_deg=180), (-1.0, 0))
+    True
+    >>> np.allclose(rotate_2d(x=1, y=0, angle_deg=360), (1.0, 0))
+    True
     '''
     angle_rad = angle_deg * pi / 180
     ch = cos(angle_rad)
@@ -2455,7 +2455,8 @@ def _discretize_linear(fun, a, b, tol=0.005, n=5):
     err0 = inf
     err = 10000
     nmax = 2 ** 20
-    while (err != err0 and err > tol and n < nmax):
+    num_tries = 0
+    while (num_tries < 5 and err > tol and n < nmax):
         err0 = err
         x0 = x
         y0 = y
@@ -2464,6 +2465,7 @@ def _discretize_linear(fun, a, b, tol=0.005, n=5):
         y = fun(x)
         y00 = interp(x, x0, y0)
         err = 0.5 * amax(np.abs((y00 - y) / (np.abs(y00 + y) + _TINY)))
+        num_tries += int(abs (err - err0) <= tol/2)
     return x, y
 
 
@@ -2479,9 +2481,9 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5):
     erri = hstack((zeros((n2, 1)), ones((n2, 1)))).ravel()
     err = erri.max()
     err0 = inf
-    # while (err != err0 and err > tol and n < nmax):
+    num_tries = 0
     for j in range(50):
-        if err != err0 and np.any(erri > tol):
+        if num_tries < 5 and np.any(erri > tol):
             err0 = err
             # find top errors
 
@@ -2490,10 +2492,9 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5):
             y = (vstack(((x[I] + x[I - 1]) / 2,
                          (x[I + 1] + x[I]) / 2)).T).ravel()
             fy = fun(y)
-
             fy0 = interp(y, x, fx)
-            erri = 0.5 * (np.abs((fy0 - fy) / (np.abs(fy0 + fy) + _TINY)))
 
+            erri = 0.5 * (np.abs((fy0 - fy) / (np.abs(fy0 + fy) + _TINY)))
             err = erri.max()
 
             x = hstack((x, y))
@@ -2502,7 +2503,7 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5):
             x = x[I]
             erri = hstack((zeros(len(fx)), erri))[I]
             fx = hstack((fx, fy))[I]
-
+            num_tries += int(abs (err - err0) <= tol/2)
         else:
             break
     else:

wafo/numba_misc.py

@@ -8,7 +8,7 @@ from numba import jit, float64, int64, int32, int8, void
 import numpy as np
 
 
-@jit(int32(int32[:], int8[:]))
+@jit(int64(int64[:], int8[:]))
 def _findcross(ind, y):
     ix, dcross, start, v = 0, 0, 0, 0
     n = len(y)
@@ -44,7 +44,7 @@ def _findcross(ind, y):
 def findcross(xn):
     '''Return indices to zero up and downcrossings of a vector
     '''
-    ind = np.empty(len(xn), dtype=np.int32)
+    ind = np.empty(len(xn), dtype=np.int64)
     m = _findcross(ind, xn)
     return ind[:m]
 

wafo/stats/estimation.py

@@ -329,18 +329,21 @@ class Profile(object):
     def _correct_Lmax(self, Lmax, free_par, fix_par):
 
         if Lmax > self.Lmax:  # foundNewphat = True
-
-            dL = self.Lmax - Lmax
-            self.alpha_cross_level -= dL
+            par_old = str(self._par)
+            dL = Lmax - self.Lmax
+            self.alpha_cross_level += dL
             self.Lmax = Lmax
             par = self._par.copy()
             par[self.i_free] = free_par
             par[self.i_fixed] = fix_par
             self.best_par = par
             self._par = par
+
             warnings.warn(
                 'The fitted parameters does not provide the optimum fit. ' +
-                'Something wrong with fit (par = {})'.format(str(par)))
+                'Something wrong with fit ' +
+                '(par = {}, par_old= {}, dl = {})'.format(str(par), par_old,
+                                                          dL))
 
     def _profile_optimum(self, phatfree0, p_opt):
         phatfree = optimize.fmin(self._profile_fun, phatfree0, args=(p_opt,),
@@ -438,17 +441,16 @@ class Profile(object):
     def _search_pmin(self, phatfree0, p_min0, p_opt):
         phatfree = phatfree0.copy()
 
-        dp = (p_opt - p_min0)/40
-        if dp < 1e-2:
-            dp = 0.1
-        p_min_opt = p_opt - dp
+        dp = np.maximum((p_opt - p_min0)/40, 0.01)*10
         Lmax, phatfree = self._profile_optimum(phatfree, p_opt)
-        for _j in range(50):
+        p_min_opt = p_min0
+        for j in range(51):
             p_min = p_opt - dp
             Lmax, phatfree = self._profile_optimum(phatfree, p_min)
+            # print((dp, p_min, p_min_opt, Lmax))
             if np.isnan(Lmax):
                 dp *= 0.33
-            elif Lmax < self.alpha_cross_level - self.alpha_Lrange * 5:
+            elif Lmax < self.alpha_cross_level - self.alpha_Lrange*5*(j+1):
                 p_min_opt = p_min
                 dp *= 0.33
             elif Lmax < self.alpha_cross_level:
@@ -456,6 +458,10 @@ class Profile(object):
                 break
             else:
                 dp *= 1.67
+        else:
+            msg = 'Exceeded max iterations. (p_min0={}, p_min={}, p={})'
+            warnings.warn(msg.format(p_min0, p_min_opt, p_opt))
+        # print('search_pmin iterations={}'.format(j))
         return p_min_opt
 
     def _search_pmax(self, phatfree0, p_max0, p_opt):
@@ -464,14 +470,14 @@ class Profile(object):
         dp = (p_max0 - p_opt)/40
         if dp < 1e-2:
             dp = 0.1
-        p_max_opt = p_opt + dp
         Lmax, phatfree = self._profile_optimum(phatfree, p_opt)
-        for _j in range(50):
+        p_max_opt = p_opt
+        for j in range(51):
             p_max = p_opt + dp
             Lmax, phatfree = self._profile_optimum(phatfree, p_max)
             if np.isnan(Lmax):
                 dp *= 0.33
-            elif Lmax < self.alpha_cross_level - self.alpha_Lrange * 2:
+            elif Lmax < self.alpha_cross_level - self.alpha_Lrange*5*(j+1):
                 p_max_opt = p_max
                 dp *= 0.33
             elif Lmax < self.alpha_cross_level:
@@ -479,6 +485,10 @@ class Profile(object):
                 break
             else:
                 dp *= 1.67
+        else:
+            msg = 'Exceeded max iterations. (p={}, p_max={}, p_max0 = {})'
+            warnings.warn(msg.format(p_opt, p_max_opt, p_max0))
+        # print('search_pmax iterations={}'.format(j))
         return p_max_opt
 
     def _profile_fun(self, free_par, fix_par):
@@ -498,9 +508,8 @@ class Profile(object):
         '''Return confidence interval for profiled parameter
         '''
         if alpha < self.alpha:
-            warnings.warn(
-                'Might not be able to return bounds with alpha less than %g' %
-                self.alpha)
+            msg = 'Might not be able to return bounds with alpha less than {}'
+            warnings.warn(msg.format(self.alpha))
         cross_level = self.Lmax - 0.5 * chi2isf(alpha, 1)
         ind = findcross(self.data, cross_level)
         n = len(ind)
@@ -511,8 +520,8 @@ class Profile(object):
 
         elif n == 1:
             x0 = ecross(self.args, self.data, ind, cross_level)
-            isUpcrossing = self.data[ind] > self.data[ind + 1]
-            if isUpcrossing:
+            is_upcrossing = self.data[ind] < self.data[ind + 1]
+            if is_upcrossing:
                 bounds = (x0, self.pmax)
                 warnings.warn('Upper bound is larger')
             else:
@@ -1399,12 +1408,19 @@ class FitDistribution(rv_frozen):
                 fixstr = 'Fixed: phat[{0:s}] = {1:s} '.format(phatistr,
                                                               phatvstr)
 
-        subtxt = 'Fit method: {0:s}, Fit p-value: {1:2.2f} {2:s}, phat=[{3:s}]'
+        subtxt = ('Fit method: {0:s}, Fit p-value: {1:2.2f} {2:s}, ' +
+                  'phat=[{3:s}], {4:s}')
         par_txt = ('{:1.2g}, '*len(self.par))[:-2].format(*self.par)
         try:
-            fig.text(0.05, 0.01, subtxt.format(self.method.upper(),
-                                               self.pvalue, fixstr, par_txt))
+            LL_txt = 'Lps_max={:2.2g}, Ll_max={:2.2g}'.format(self.LPSmax,
+                                                              self.LLmax)
         except Exception:
+            LL_txt = 'Lps_max={}, Ll_max={}'.format(self.LPSmax, self.LLmax)
+        try:
+            fig.text(0.05, 0.01, subtxt.format(self.method.upper(),
+                                               self.pvalue, fixstr, par_txt,
+                                               LL_txt))
+        except AttributeError:
             pass
 
     def plotesf(self, symb1='r-', symb2='b.', axis=None, plot_ci=False):