From 9d37eb4f679172222ab0d07090729e97fb40380f Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Fri, 18 Nov 2016 20:31:50 +0100 Subject: [PATCH] updated doctests + fitting --- wafo/misc.py | 53 ++++++++++++++++++------------------ wafo/numba_misc.py | 4 +-- wafo/stats/estimation.py | 58 +++++++++++++++++++++++++--------------- 3 files changed, 66 insertions(+), 49 deletions(-) diff --git a/wafo/misc.py b/wafo/misc.py index df7aa1f..a323121 100644 --- a/wafo/misc.py +++ b/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: diff --git a/wafo/numba_misc.py b/wafo/numba_misc.py index d16690c..496e504 100644 --- a/wafo/numba_misc.py +++ b/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] diff --git a/wafo/stats/estimation.py b/wafo/stats/estimation.py index def01f9..dde85ea 100644 --- a/wafo/stats/estimation.py +++ b/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):