From a2639888886e069fa040361544b8240501709930 Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Fri, 21 Dec 2018 15:54:26 +0100 Subject: [PATCH] Added lfind to misc.py Replaced call to valarray with np.full Refactored code from detrendma into a new function moving_average. Made doctests more robust. pep8 --- wafo/dctpack.py | 4 +- wafo/misc.py | 386 +++++++++++++++++++++++++++--------------- wafo/spectrum/core.py | 66 ++++---- wafo/stats/core.py | 95 +++++------ 4 files changed, 327 insertions(+), 224 deletions(-) diff --git a/wafo/dctpack.py b/wafo/dctpack.py index 9762e49..0679c8b 100644 --- a/wafo/dctpack.py +++ b/wafo/dctpack.py @@ -346,8 +346,8 @@ def shiftdim(x, n=None): if n is None: return x.reshape(no_leading_ones(x.shape)) elif n >= 0: - return x.transpose(np.roll(range(x.ndim), -n)) - return x.reshape((1,) * -n + x.shape) + return x.transpose(np.roll(range(x.ndim), -1 * n)) + return x.reshape((1,) * abs(n) + x.shape) def example_dct2(debug=True): diff --git a/wafo/misc.py b/wafo/misc.py index 9e09b0f..bd6abe6 100644 --- a/wafo/misc.py +++ b/wafo/misc.py @@ -7,14 +7,13 @@ import collections from wafo import numba_misc import fractions import numpy as np -from numpy import ( - meshgrid, - amax, logical_and, arange, linspace, atleast_1d, - asarray, ceil, floor, frexp, hypot, - sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff, - inf, pi, interp, isscalar, zeros, ones, - sign, unique, hstack, vstack, nonzero, where, extract) -from scipy.special import gammaln +from numpy import (amax, logical_and, arange, linspace, atleast_1d, + asarray, ceil, floor, frexp, hypot, + sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff, + inf, pi, interp, isscalar, zeros, ones, + sign, unique, hstack, vstack, nonzero, where, extract, + meshgrid) +from scipy.special import gammaln # pylint: disable=no-name-in-module from scipy.integrate import trapz, simps import warnings from time import strftime, gmtime @@ -25,7 +24,6 @@ try: except ImportError: warnings.warn('c_library not found. Check its compilation.') clib = None - FLOATINFO = np.finfo(float) _TINY = FLOATINFO.tiny _EPS = FLOATINFO.eps @@ -41,8 +39,10 @@ __all__ = ['now', 'spaceline', 'narg_smallest', 'args_flat', 'is_numlike', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere', 'lazyselect', + 'valarray', 'lfind', + 'moving_average', 'piecewise', - 'valarray', 'check_random_state'] + 'check_random_state'] def xor(a, b): @@ -52,6 +52,41 @@ def xor(a, b): return a ^ b +def lfind(haystack, needle, maxiter=None): + """Return indices to the maxiter first needles in the haystack as an iterable generator. + + Parameters + ---------- + haystack: list or tuple of items + needle: item to find + maxiter: scalar integer maximum number of occurences + + Returns + ------- + indices_gen: generator object + + Example + ------ + >>> haystack = (1, 3, 4, 3, 10, 3, 1, 2, 5, 7, 10) + >>> list(lfind(haystack, 3)) + [1, 3, 5] + >>> [i for i in lfind(haystack, 3, 2)] + [1, 3] + >>> [i for i in lfind(haystack, 0, 2)] + [] + """ + maxiter = inf if maxiter is None else maxiter + ix = -1 + i = 0 + while i < maxiter: + i += 1 + try: + ix = haystack.index(needle, ix + 1) + yield ix + except (ValueError, KeyError): + break + + def check_random_state(seed): """Turn seed into a np.random.RandomState instance @@ -82,16 +117,10 @@ def check_random_state(seed): raise ValueError(msg.format(seed)) -def valarray(shape, value=np.NaN, typecode=None): +def valarray(shape, value=np.nan, typecode=None): """Return an array of all value. """ - if typecode is None: - typecode = bool - out = ones(shape, dtype=typecode) * value - - if not isinstance(out, np.ndarray): - out = asarray(out) - return out + return np.full(shape, fill_value=value, dtype=typecode) def piecewise(condlist, funclist, xi=None, fillvalue=0.0, args=(), **kw): @@ -173,20 +202,22 @@ def piecewise(condlist, funclist, xi=None, fillvalue=0.0, args=(), **kw): Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for ``x >= 0``. - >>> piecewise([x < 0, x >= 0], [lambda x: -x, lambda x: x], xi=(x,)) - array([ 2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) + >>> np.allclose(piecewise([x < 0, x >= 0], [lambda x: -x, lambda x: x], xi=(x,)), + ... [ 2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) + True Define the absolute value, which is ``-x*y`` for ``x*y <0`` and ``x*y`` for ``x*y >= 0`` >>> X, Y = np.meshgrid(x, x) - >>> piecewise([X * Y < 0, ], [lambda x, y: -x * y, lambda x, y: x * y], - ... xi=(X, Y)) - array([[ 6.25, 3.75, 1.25, 1.25, 3.75, 6.25], - [ 3.75, 2.25, 0.75, 0.75, 2.25, 3.75], - [ 1.25, 0.75, 0.25, 0.25, 0.75, 1.25], - [ 1.25, 0.75, 0.25, 0.25, 0.75, 1.25], - [ 3.75, 2.25, 0.75, 0.75, 2.25, 3.75], - [ 6.25, 3.75, 1.25, 1.25, 3.75, 6.25]]) + >>> np.allclose(piecewise([X * Y < 0, ], [lambda x, y: -x * y, + ... lambda x, y: x * y], xi=(X, Y)), + ... [[ 6.25, 3.75, 1.25, 1.25, 3.75, 6.25], + ... [ 3.75, 2.25, 0.75, 0.75, 2.25, 3.75], + ... [ 1.25, 0.75, 0.25, 0.25, 0.75, 1.25], + ... [ 1.25, 0.75, 0.25, 0.25, 0.75, 1.25], + ... [ 3.75, 2.25, 0.75, 0.75, 2.25, 3.75], + ... [ 6.25, 3.75, 1.25, 1.25, 3.75, 6.25]]) + True """ def otherwise_condition(condlist): return ~np.logical_or.reduce(condlist, axis=0) @@ -201,6 +232,7 @@ def piecewise(condlist, funclist, xi=None, fillvalue=0.0, args=(), **kw): condlist = np.broadcast_arrays(*condlist) if len(condlist) == len(funclist) - 1: condlist.append(otherwise_condition(condlist)) + if xi is None: arrays = () dtype = np.result_type(*funclist) @@ -212,7 +244,7 @@ def piecewise(condlist, funclist, xi=None, fillvalue=0.0, args=(), **kw): dtype = np.result_type(*arrays) shape = arrays[0].shape - out = valarray(shape, fillvalue, dtype) + out = np.full(shape, fillvalue, dtype) for cond, func in zip(condlist, funclist): if cond.any(): if isinstance(func, collections.Callable): @@ -231,12 +263,15 @@ def lazywhere(cond, arrays, f, fillvalue=None, f2=None): >>> a, b = np.array([1, 2, 3, 4]), np.array([5, 6, 7, 8]) >>> def f(a, b): ... return a*b - >>> lazywhere(a > 2, (a, b), f, np.nan) - array([ nan, nan, 21., 32.]) + >>> r = lazywhere(a > 2, (a, b), f, np.nan) + >>> np.allclose(r[2:], [21., 32.]); np.all(np.isnan(r[:2])) + True + True >>> def f2(a, b): ... return (a*b)**2 - >>> lazywhere(a > 2, (a, b), f, f2=f2) - array([ 25., 144., 21., 32.]) + >>> np.allclose(lazywhere(a > 2, (a, b), f, f2=f2), + ... [ 25., 144., 21., 32.]) + True Notice it assumes that all `arrays` are of the same shape, or can be broadcasted together. @@ -250,7 +285,7 @@ def lazywhere(cond, arrays, f, fillvalue=None, f2=None): arrays = np.broadcast_arrays(*arrays) temp = tuple(np.extract(cond, arr) for arr in arrays) - out = valarray(np.shape(arrays[0]), value=fillvalue) + out = np.full(np.shape(arrays[0]), fill_value=fillvalue) np.place(out, cond, f(*temp)) if f2 is not None: temp = tuple(np.extract(~cond, arr) for arr in arrays) @@ -309,11 +344,11 @@ def rotation_matrix(heading, pitch, roll): Examples -------- >>> import numpy as np - >>> rotation_matrix(heading=0, pitch=0, roll=0) - array([[ 1., 0., 0.], - [ 0., 1., 0.], - [ 0., 0., 1.]]) - + >>> np.allclose(rotation_matrix(heading=0, pitch=0, roll=0), + ... [[ 1., 0., 0.], + ... [ 0., 1., 0.], + ... [ 0., 0., 1.]]) + True >>> np.allclose(rotation_matrix(heading=180, pitch=0, roll=0), ... [[ -1., 0., 0.], ... [ 0., -1., 0.], @@ -438,12 +473,20 @@ def spaceline(start_point, stop_point, num=10): Example ------- >>> import wafo.misc as pm - >>> pm.spaceline((2,0,0), (3,0,0), num=5) - array([[ 2. , 0. , 0. ], - [ 2.25, 0. , 0. ], - [ 2.5 , 0. , 0. ], - [ 2.75, 0. , 0. ], - [ 3. , 0. , 0. ]]) + >>> np.allclose(pm.spaceline((2,0,0), (3,0,0), num=5), + ... [[ 2. , 0. , 0. ], + ... [ 2.25, 0. , 0. ], + ... [ 2.5 , 0. , 0. ], + ... [ 2.75, 0. , 0. ], + ... [ 3. , 0. , 0. ]]) + True + >>> np.allclose(pm.spaceline((2,0,0), (0,0,3), num=5), + ... [[ 2. , 0. , 0. ], + ... [ 1.5 , 0. , 0.75], + ... [ 1. , 0. , 1.5 ], + ... [ 0.5 , 0. , 2.25], + ... [ 0. , 0. , 3. ]]) + True ''' num = int(num) start, stop = np.atleast_1d(start_point, stop_point) @@ -615,6 +658,8 @@ def is_numlike(obj): True >>> is_numlike('1') False + >>> is_numlike([1]) + False """ try: obj + 1 @@ -670,13 +715,13 @@ class Bunch(object): Example ------- >>> d = Bunch(test1=1,test2=3) - >>> d.test1 - 1 + >>> (d.test1, d.test2) + (1, 3) >>> sorted(d.keys()) == ['test1', 'test2'] True - >>> d.update(test1=2) - >>> d.test1 - 2 + >>> d.update(test1=5) + >>> (d.test1, d.test2) + (5, 3) ''' def __init__(self, **kwargs): @@ -703,7 +748,7 @@ def sub_dict_select(somedict, somekeys): # the keyword exists in options >>> opt = dict(arg1=2, arg2=3) >>> kwds = dict(arg2=100,arg3=1000) - >>> sub_dict = sub_dict_select(kwds,opt.keys()) + >>> sub_dict = sub_dict_select(kwds, opt) >>> opt.update(sub_dict) >>> opt == {'arg1': 2, 'arg2': 100} True @@ -713,7 +758,7 @@ def sub_dict_select(somedict, somekeys): dict_intersection ''' # slower: validKeys = set(somedict).intersection(somekeys) - return dict((k, somedict[k]) for k in somekeys if k in somedict) + return type(somedict)((k, somedict[k]) for k in somekeys if k in somedict) def parse_kwargs(options, **kwargs): @@ -731,13 +776,91 @@ def parse_kwargs(options, **kwargs): See also sub_dict_select ''' - newopts = sub_dict_select(kwargs, options.keys()) + newopts = sub_dict_select(kwargs, options) if len(newopts) > 0: options.update(newopts) return options -def detrendma(x, L): +def moving_average(x, L, axis=0): + """ + Return moving average from data using a window of size 2*L+1. + If 2*L+1 > len(x) then the mean is returned + + Parameters + ---------- + x : vector or matrix of column vectors + of data + L : scalar, integer + defines the size of the moving average window + axis: scalar integer + axis along which the moving average is computed. Default axis=0. + + Returns + ------- + y : ndarray + moving average + + Examples + -------- + >>> import matplotlib.pyplot as plt + >>> import numpy as np + >>> exp = np.exp; cos = np.cos; randn = np.random.randn + + >>> x = np.linspace(0,1,200) + >>> y = exp(x)+cos(5*2*pi*x)+1e-1*randn(x.size) + >>> y0 = moving_average(y, 20) + >>> np.allclose(y0[:4], [ 1.1189971, 1.1189971, 1.1189971, 1.1189971], atol=1e-1) + True + + >>> x2 = np.linspace(1, 5, 5) + >>> np.allclose(moving_average(x2, L=1), [2., 2., 3., 4., 4.]) + True + >>> np.allclose(moving_average(x2, L=10), [3., 3., 3., 3., 3.]) + True + >>> x3 = np.vstack((x2, x2+5)) + >>> np.allclose(moving_average(x3, L=1, axis=1), [[2., 2., 3., 4., 4.], + ... [7., 7., 8., 9., 9.]]) + True + >>> np.allclose(moving_average(x3, L=10, axis=1), [[3., 3., 3., 3., 3.], + ... [8., 8., 8., 8., 8.]]) + True + + h = plt.plot(x, y, x, y0, 'r', x, exp(x), 'k', x, tr, 'm') + plt.close('all') + + See also + -------- + Reconstruct + """ + _assert(0 < L, 'L must be positive') + _assert(L == np.round(L), 'L must be an integer') + + x1 = np.atleast_1d(x) + axes = np.arange(x.ndim) + + axes[0] = axis + axes[axis] = 0 + + x1 = np.transpose(x1, axes) + y = np.empty_like(x1) + + n = x1.shape[0] + if n < 2 * L + 1: # only able to remove the mean + y[:] = x1.mean(axis=0) + return np.transpose(y, axes) + + mn = x1[0:2 * L + 1].mean(axis=0) + y[0:L + 1] = mn + + ix = np.r_[L + 1:(n - L)] + y[ix] = ((x1[ix + L] - x1[ix - L]) / (2 * L)).cumsum(axis=0) + mn + y[n - L::] = y[n - L - 1] + + return np.transpose(y, axes) + + +def detrendma(x, L, axis=0): """ Removes a trend from data using a moving average of size 2*L+1. If 2*L+1 > len(x) then the mean is removed @@ -748,6 +871,8 @@ def detrendma(x, L): of data L : scalar, integer defines the size of the moving average window + axis: scalar integer + axis along which the moving average is computed. Default axis=0. Returns ------- @@ -760,11 +885,12 @@ def detrendma(x, L): >>> import numpy as np >>> import wafo.misc as wm >>> exp = np.exp; cos = np.cos; randn = np.random.randn + >>> x = np.linspace(0,1,200) >>> noise = 0.1*randn(x.size) >>> noise = 0.1*np.sin(100*x) >>> y = exp(x)+cos(5*2*pi*x) + noise - >>> y0 = wm.detrendma(y,20) + >>> y0 = wm.detrendma(y, 20) >>> tr = y-y0 >>> np.allclose(tr[:5], ... [ 1.14134814, 1.14134814, 1.14134814, 1.14134814, 1.14134814]) @@ -781,28 +907,11 @@ def detrendma(x, L): See also -------- - Reconstruct + Reconstruct, moving_average """ - _assert(0 < L, 'L must be positive') - _assert(L == np.round(L), 'L must be an integer') - - x_1 = np.atleast_1d(x) - if x_1.shape[0] == 1: - x_1 = x_1.ravel() - - n = x_1.shape[0] - if n < 2 * L + 1: # only able to remove the mean - return x_1 - x_1.mean(axis=0) - - mean = x_1[:2 * L + 1].mean(axis=0) - y = np.empty_like(x_1) - y[:L + 1] = x_1[:L + 1] - mean - - i = np.r_[L + 1:(n - L)] - trend = ((x_1[i + L] - x_1[i - L]) / (2 * L)).cumsum(axis=0) + mean - y[i] = x_1[i] - trend - y[n - L::] = x_1[n - L::] - trend[-1] - return y + x1 = np.atleast_1d(x) + trend = moving_average(x1, L, axis) + return x1 - trend def ecross(t, f, ind, v=0): @@ -860,7 +969,7 @@ def _findcross(x, method='numba'): '''Return indices to zero up and downcrossings of a vector ''' if clib is not None and method == 'clib': - ind, m = clib.findcross(x, 0.0) + ind, m = clib.findcross(x, 0.0) # pylint: disable=no-member return ind[:int(m)] return numba_misc.findcross(x) @@ -903,6 +1012,7 @@ def findcross(x, v=0.0, kind=None, method='clib'): >>> ind = wm.findcross(x,v) # all crossings >>> np.allclose(ind, [ 9, 25, 80, 97, 151, 168, 223, 239]) True + >>> ind2 = wm.findcross(x,v,'u') >>> np.allclose(ind2, [ 9, 80, 151, 223]) True @@ -940,7 +1050,6 @@ def findcross(x, v=0.0, kind=None, method='clib'): # if kind=='dw' or kind=='tw' # or that the first is a level v up-crossing # if kind=='uw' or kind=='cw' - first_is_down_crossing = int(xn[ind[0]] > xn[ind[0] + 1]) if xor(first_is_down_crossing, kind in ('dw', 'tw')): ind = ind[1::] @@ -1171,7 +1280,7 @@ def findrfc(tp, h=0.0, method='clib'): t_start = int(y[0] > y[1]) # first is a max, ignore it y = y[t_start::] n = len(y) - NC = np.floor(n / 2) - 1 + NC = np.floor(n // 2) - 1 if (NC < 1): return zeros(0, dtype=np.int) # No RFC cycles*/ @@ -1182,7 +1291,7 @@ def findrfc(tp, h=0.0, method='clib'): return zeros(0, dtype=np.int) if clib is not None and method == 'clib': - ind, ix = clib.findrfc(y, h) + ind, ix = clib.findrfc(y, h) # pylint: disable=no-member ix = int(ix) else: if isinstance(method, str): @@ -1218,25 +1327,23 @@ def rfcfilter(x, h, method=0): >>> import wafo.misc as wm >>> x = data.sea() >>> y = wm.rfcfilter(x[:,1], h=0, method=1) - >>> np.all(np.abs(y[0:5]-np.array([-1.2004945 , 0.83950546, -0.09049454, - ... -0.02049454, -0.09049454]))<1e-7) + >>> np.allclose(y[0:5], [-1.2004945 , 0.83950546, -0.09049454, -0.02049454, -0.09049454]) + True + >>> np.allclose(y.shape, (2172,)) True - >>> y.shape - (2172,) # 2. This removes all rainflow cycles with range less than 0.5. >>> y1 = wm.rfcfilter(x[:,1], h=0.5) - >>> y1.shape - (863,) - >>> np.all(np.abs(y1[0:5]-np.array([-1.2004945 , 0.83950546, -0.43049454, - ... 0.34950546, -0.51049454]))<1e-7) + >>> np.allclose(y1.shape, (863,)) + True + >>> np.allclose(y1[0:5], [-1.2004945 , 0.83950546, -0.43049454, 0.34950546, -0.51049454]) True >>> ind = wm.findtp(x[:,1], h=0.5) >>> y2 = x[ind,1] - >>> y2[0:5] - array([-1.2004945 , 0.83950546, -0.43049454, 0.34950546, -0.51049454]) - >>> y2[-5::] - array([ 0.83950546, -0.64049454, 0.65950546, -1.0004945 , 0.91950546]) + >>> np.allclose(y2[0:5], [-1.2004945 , 0.83950546, -0.43049454, 0.34950546, -0.51049454]) + True + >>> np.allclose(y2[-5::], [ 0.83950546, -0.64049454, 0.65950546, -1.0004945 , 0.91950546]) + True See also -------- @@ -1609,10 +1716,10 @@ def common_shape(*args, ** kwds): >>> A = np.ones((4,1)) >>> B = 2 >>> C = np.ones((1,5))*5 - >>> wm.common_shape(A,B,C) - (4, 5) - >>> wm.common_shape(A,B,C,shape=(3,4,1)) - (3, 4, 5) + >>> np.allclose(wm.common_shape(A,B,C), (4, 5)) + True + >>> np.allclose(wm.common_shape(A,B,C,shape=(3,4,1)), (3, 4, 5)) + True See also -------- @@ -1652,12 +1759,12 @@ def argsreduce(condition, * args): >>> C = rand((1,5)) >>> cond = np.ones(A.shape) >>> [A1,B1,C1] = wm.argsreduce(cond,A,B,C) - >>> B1.shape - (20,) + >>> np.allclose(B1.shape, (20,)) + True >>> cond[2,:] = 0 >>> [A2,B2,C2] = wm.argsreduce(cond,A,B,C) - >>> B2.shape - (15,) + >>> np.allclose(B2.shape, (15,)) + True See also -------- @@ -1794,7 +1901,7 @@ def getshipchar(**ship_property): ... 'lengthSTD': 2.0113098831942762, ... 'draught': 9.5999999999999996, ... 'propeller_diameterSTD': 0.20267047566705432, - ... 'max_deadweightSTD': 3096.9000000000001, + ... 'max_deadweight_std': 3096.9000000000001, ... 'beam': 29.0, 'length': 216.0, ... 'beamSTD': 2.9000000000000004, ... 'service_speed': 10.0, @@ -1818,7 +1925,7 @@ def getshipchar(**ship_property): ''' max_deadweight, prop = _get_max_deadweight(**ship_property) - propertySTD = prop + 'STD' + property_std = prop + 'STD' length = np.round(3.45 * max_deadweight ** 0.40) length_err = length ** 0.13 @@ -1837,18 +1944,16 @@ def getshipchar(**ship_property): p_diam_err = 0.12 * length_err ** (3.0 / 4.0) max_deadweight = np.round(max_deadweight) - max_deadweightSTD = 0.1 * max_deadweight + max_deadweight_std = 0.1 * max_deadweight shipchar = dict(beam=beam, beamSTD=beam_err, draught=draught, draughtSTD=draught_err, length=length, lengthSTD=length_err, - max_deadweight=max_deadweight, - max_deadweightSTD=max_deadweightSTD, - propeller_diameter=p_diam, - propeller_diameterSTD=p_diam_err, + max_deadweight=max_deadweight, max_deadweightSTD=max_deadweight_std, + propeller_diameter=p_diam, propeller_diameterSTD=p_diam_err, service_speed=speed, service_speedSTD=speed_err) - shipchar[propertySTD] = 0 + shipchar[property_std] = 0 return shipchar @@ -1874,8 +1979,8 @@ def binomln(z, w): Example ------- - >>> np.abs(binomln(3,2)- 1.09861229)<1e-7 - array([ True], dtype=bool) + >>> np.allclose(binomln(3,2), 1.09861229) + True See also -------- @@ -1918,12 +2023,11 @@ def betaloge(z, w): -------- betaln, beta ''' - # y = gammaln(z)+gammaln(w)-gammaln(z+w) zpw = z + w - return (stirlerr(z) + stirlerr(w) + 0.5 * log(2 * pi) + - (w - 0.5) * log(w) + (z - 0.5) * log(z) - stirlerr(zpw) - - (zpw - 0.5) * log(zpw)) + return (stirlerr(z) + stirlerr(w) + 0.5 * log(2 * pi) + (w - 0.5) * log(w) + + (z - 0.5) * log(z) - stirlerr(zpw) - (zpw - 0.5) * log(zpw)) + # y = gammaln(z)+gammaln(w)-gammaln(z+w) # stirlings approximation: # (-(zpw-0.5).*log(zpw) +(w-0.5).*log(w)+(z-0.5).*log(z) +0.5*log(2*pi)) # return y @@ -1973,8 +2077,8 @@ def gravity(phi=45): ''' phir = phi * pi / 180. # change from degrees to radians - return 9.78049 * (1. + 0.0052884 * sin(phir) ** 2. - - 0.0000059 * sin(2 * phir) ** 2.) + return 9.78049 * (1. + 0.0052884 * sin(phir) ** 2. + - 0.0000059 * sin(2 * phir) ** 2.) def nextpow2(x): @@ -2035,7 +2139,6 @@ def discretize(fun, a, b, tol=0.005, n=5, method='linear'): >>> np.allclose(xa[:5], [0., 0.19634954, 0.39269908, 0.58904862, 0.78539816]) True - t = plt.plot(x, y, xa, ya, 'r.') plt.show() plt.close('all') @@ -2099,7 +2202,6 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5): abserr = np.abs(fy0 - fy) erri = 0.5 * (abserr / (np.abs(fy0) + np.abs(fy) + _TINY + tol)) - # abserr = np.abs(fy0 - fy) # converged = abserr <= np.maximum(abseps, releps * abs(fy)) # converged = abserr <= np.maximum(tol, tol * abs(fy)) err = erri.max() @@ -2295,7 +2397,7 @@ def tranproc(x, f, x0, *xi): x,f : array-like [x,f(x)], transform function, y = f(x). x0, x1,...,xn : vectors - where xi is the i'th time derivative of x0. 0<=N<=4. + where xi is the i'th time derivative of x0. 0<=n<=4. Returns ------- @@ -2337,21 +2439,21 @@ def tranproc(x, f, x0, *xi): -------- trangood. """ - def _default_step(xo, N): + def _default_step(xo, num_derivatives): hn = xo[1] - xo[0] - if hn ** N < sqrt(_EPS): + if hn ** num_derivatives < sqrt(_EPS): msg = ('Numerical problems may occur for the derivatives in ' + 'tranproc.\n' + 'The sampling of the transformation may be too small.') warnings.warn(msg) return hn - def _diff(xo, fo, x0, N): - hn = _default_step(xo, N) + def _diff(xo, fo, x0, num_derivatives): + hn = _default_step(xo, num_derivatives) # Transform X with the derivatives of f. fder = vstack((xo, fo)) - fxder = zeros((N, x0.size)) - for k in range(N): # Derivation of f(x) using a difference method. + fxder = zeros((num_derivatives, x0.size)) + for k in range(num_derivatives): # Derivation of f(x) using a difference method. n = fder.shape[-1] fder = vstack([(fder[0, 0:n - 1] + fder[0, 1:n]) / 2, diff(fder[1, :]) / hn]) @@ -2386,8 +2488,8 @@ def tranproc(x, f, x0, *xi): xi = atleast_1d(*xi) if not isinstance(xi, list): xi = [xi, ] - N = len(xi) # N = number of derivatives - nmax = ceil((xo.ptp()) * 10 ** (7. / max(N, 1))) + num_derivatives = len(xi) # num_derivatives = number of derivatives + nmax = ceil((xo.ptp()) * 10 ** (7. / max(num_derivatives, 1))) xo, fo = trangood(xo, fo, min_x=min(x0), max_x=max(x0), max_n=nmax) n = f.shape[0] @@ -2400,16 +2502,16 @@ def tranproc(x, f, x0, *xi): y0 = fo[fi] + (fo[fi + 1] - fo[fi]) * xu y = y0 - if N > 4: + if num_derivatives > 4: warnings.warn('Transformation of derivatives of order>4 is ' + 'not supported.') - N = 4 - if N > 0: + num_derivatives = 4 + if num_derivatives > 0: y = [y0] - fxder = _diff(xo, fo, x0, N) + fxder = _diff(xo, fo, x0, num_derivatives) # Calculate the transforms of the derivatives of X. dfuns = [_der_1, _der_2, _der_3, _der_4] - for dfun in dfuns[:N]: + for dfun in dfuns[:num_derivatives]: y.append(dfun(fxder, xi)) return y @@ -2437,14 +2539,18 @@ def good_bins(data=None, range=None, num_bins=None, odd=False, loose=True): # @ Example ------- >>> import wafo.misc as wm - >>> wm.good_bins(range=(0,5), num_bins=6) - array([-1., 0., 1., 2., 3., 4., 5., 6.]) - >>> wm.good_bins(range=(0,5), num_bins=6, loose=False) - array([ 0., 1., 2., 3., 4., 5.]) - >>> wm.good_bins(range=(0,5), num_bins=6, odd=True) - array([-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]) - >>> wm.good_bins(range=(0,5), num_bins=6, odd=True, loose=False) - array([-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5]) + >>> np.allclose(wm.good_bins(range=(0,5), num_bins=6), + ... [-1., 0., 1., 2., 3., 4., 5., 6.]) + True + >>> np.allclose(wm.good_bins(range=(0,5), num_bins=6, loose=False), + ... [ 0., 1., 2., 3., 4., 5.]) + True + >>> np.allclose(wm.good_bins(range=(0,5), num_bins=6, odd=True), + ... [-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5]) + True + >>> np.allclose(wm.good_bins(range=(0,5), num_bins=6, odd=True, loose=False), + ... [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5]) + True ''' def _default_range(range_, x): return range_ if range_ else (x.min(), x.max()) diff --git a/wafo/spectrum/core.py b/wafo/spectrum/core.py index 1d95db4..0b2ddda 100644 --- a/wafo/spectrum/core.py +++ b/wafo/spectrum/core.py @@ -135,7 +135,7 @@ def qtf(w, h=inf, g=9.81, method='winterstein', rtol=1e-7, atol=0): >>> np.allclose(hd2, [[-2.50000000e-101, 2.54841998e-022, 2.54841998e-012, 2.54841998e-004], ... [2.54841998e-022, -2.50000000e-101, 2.54836901e-012, 2.54841997e-004], ... [2.54841998e-012, 2.54836901e-012, -2.50000000e-101, 2.54791032e-004], - ... [2.54841998e-004, 2.54841997e-004, 2.54791032e-004, -2.50000000e-101]], + ... [2.54841998e-004, 2.54841997e-004, 2.54791032e-004, -2.500000e-101]], ... atol=0) True @@ -175,9 +175,8 @@ def qtf(w, h=inf, g=9.81, method='winterstein', rtol=1e-7, atol=0): h_dii = zeros(num_w) return h_s, h_d, h_dii - w1 = w + _TINY **(1./10) * (np.sign(w) * np.int_(np.abs(w) < _EPS) + np.int_(w==0)) - # dw = w1-w - # print(dw) + w1 = w + _TINY ** (1. / 10) * (np.sign(w) * np.int_(np.abs(w) < _EPS) + np.int_(w == 0)) + w = w1 # k_w += _TINY ** (1./3) * (np.sign(k_w) * np.int_(np.abs(k_w) < _EPS) + np.int_(k_w==0)) k_w = w2k(w, theta=0, h=h, g=g, rtol=rtol, atol=atol)[0] @@ -194,10 +193,10 @@ def qtf(w, h=inf, g=9.81, method='winterstein', rtol=1e-7, atol=0): if method.startswith('langley'): p_1 = (-2 * w1p2 * (k12 * g ** 2. - w12 ** 2.) + w_1 * (w_2 ** 4. - g ** 2 * k_2 ** 2) + - w_2 * (w_1 ** 4 - g * 2. * k_1 ** 2)) / (4. * w12 +_TINY) + w_2 * (w_1 ** 4 - g * 2. * k_1 ** 2)) / (4. * w12 + _TINY) p_2 = w1p2 ** 2. * cosh((k1p2) * h) - g * (k1p2) * sinh((k1p2) * h) - p_2 += _TINY * np.int_(p_2==0) + p_2 += _TINY * np.int_(p_2 == 0) h_s = (-p_1 / p_2 * w1p2 * cosh((k1p2) * h) / g - (k12 * g ** 2 - w12 ** 2.) / (4 * g * w12 + _TINY) + @@ -207,7 +206,7 @@ def qtf(w, h=inf, g=9.81, method='winterstein', rtol=1e-7, atol=0): w_1 * (w_2 ** 4 - g ** 2 * k_2 ** 2) + w_2 * (w_1 ** 4 - g ** 2 * k_1 ** 2)) / (4. * w12 + _TINY) p_4 = w1m2 ** 2. * cosh(k1m2 * h) - g * (k1m2) * sinh((k1m2) * h) - p_4 += _TINY * np.int_(p_4==0) + p_4 += _TINY * np.int_(p_4 == 0) h_d = (-p_3 / p_4 * (w1m2) * cosh((k1m2) * h) / g - (k12 * g ** 2 + w12 ** 2) / (4 * g * w12 + _TINY) + @@ -218,13 +217,14 @@ def qtf(w, h=inf, g=9.81, method='winterstein', rtol=1e-7, atol=0): tmp2 = (w_1 ** 2. + w_2 ** 2. + w12) / g h_s = 0.25 * ((tmp1 - tmp2 - + g * (w_1 * k_2 ** 2. + w_2 * k_1 ** 2) / (w12 * w1p2 + 0)) - / (1. - g * h * k1p2 / (w1p2 ** 2. + 0) * tanh(k1p2)) - + tmp2 - 0.5 * tmp1) # OK + + g * (w_1 * k_2 ** 2. + w_2 * k_1 ** 2) / (w12 * w1p2 + 0)) + / (1. - g * h * k1p2 / (w1p2 ** 2. + 0) * tanh(k1p2)) + + tmp2 - 0.5 * tmp1) # OK tiny_diag = _TINY * np.diag(np.ones(num_w)) # Avoid division by zero on diagonal tmp3 = (w_1 ** 2 + w_2 ** 2 - w12) / g # OK - numerator = (tmp1 - tmp3 - g * (w_1 * k_2 ** 2 - w_2 * k_1 ** 2) / (w12 * w1m2 + tiny_diag)) + numerator = (tmp1 - tmp3 - g * (w_1 * k_2 ** 2 + - w_2 * k_1 ** 2) / (w12 * w1m2 + tiny_diag)) h_d = 0.25 * (numerator / (1. - g * h * k1m2 / (w1m2 ** 2. + tiny_diag) * tanh(k1m2)) + tmp3 - 0.5 * tmp1) # OK @@ -825,7 +825,7 @@ class SpecData1D(PlotData): nt = rate * (n_f - 1) else: # check if Nt is ok nt = minimum(nt, rate * (n_f - 1)) - # nr, nt = int(nr), int(nt) + spec = self.copy() spec.resample(dt) @@ -990,10 +990,10 @@ class SpecData1D(PlotData): iseed : scalar integer starting seed number for the random number generator (default none is set) - fnLimit : real scalar + fn_limit : real scalar normalized upper frequency limit of spectrum for 2'nd order components. The frequency is normalized with - sqrt(gravity*tanh(kbar*water_depth)/Amax)/(2*pi) + sqrt(gravity*tanh(kbar*water_depth)/a_max)/(2*pi) (default sqrt(2), i.e., Convergence criterion). Generally this should be the same as used in the final non-linear simulation (see example below). @@ -1066,7 +1066,7 @@ class SpecData1D(PlotData): max_sim = 30 tolerance = 5e-4 - L = 200 # maximum lag size of the window function used in estimate + max_lag = 200 # maximum lag size of the window function used in estimate # ftype = self.freqtype #options are 'f' and 'w' and 'k' # switch ftype # case 'f', @@ -1094,9 +1094,9 @@ class SpecData1D(PlotData): # Expected maximum amplitude for 10000 waves seastate num_waves = 10000 # Typical number of waves in 30 hour seastate - Amax = sqrt(2 * log(num_waves)) * Hm0 / 4 + a_max = sqrt(2 * log(num_waves)) * Hm0 / 4 - fLimitLo = sqrt(gravity * tanh(kbar * water_depth) * Amax / water_depth ** 3) + f_limit_lo = sqrt(gravity * tanh(kbar * water_depth) * a_max / water_depth ** 3) freq = S.args eps = finfo(float).eps @@ -1105,11 +1105,11 @@ class SpecData1D(PlotData): SL = S - indZero = nonzero(freq < fLimitLo)[0] - if len(indZero): - SL.data[indZero] = 0 + ind_zero = nonzero(freq < f_limit_lo)[0] + if len(ind_zero): + SL.data[ind_zero] = 0 - maxS = max(S.data) + max_spec = max(S.data) # Fs = 2*freq(end)+eps # sampling frequency for ix in range(max_sim): @@ -1117,12 +1117,12 @@ class SpecData1D(PlotData): method=method, fnlimit=fn_limit, output='timeseries') x2.data -= x1.data # x2(:,2:end) = x2(:,2:end) -x1(:,2:end) - S2 = x2.tospecdata(L) - S1 = x1.tospecdata(L) + S2 = x2.tospecdata(max_lag) + S1 = x1.tospecdata(max_lag) # TODO: Finish spec.to_linspec -# S2 = dat2spec(x2, L) -# S1 = dat2spec(x1, L) +# S2 = dat2spec(x2, max_lag) +# S1 = dat2spec(x1, max_lag) # %[tf21,fi] = tfe(x2(:,2),x1(:,2),1024,Fs,[],512) # %Hw11 = interp1q(fi,tf21.*conj(tf21),freq) if True: @@ -1144,8 +1144,8 @@ class SpecData1D(PlotData): SL.data = (Hw1 * S.data) - if len(indZero): - SL.data[indZero] = 0 + if len(ind_zero): + SL.data[ind_zero] = 0 # end k = nonzero(SL.data < 0)[0] if len(k): # Make sure that the current guess is larger than zero @@ -1176,17 +1176,17 @@ class SpecData1D(PlotData): # figtile # end - print('Iteration : %d, Hw12 : %g Hw12/maxS : %g' % - (ix, maxHw12, (maxHw12 / maxS))) - if (maxHw12 < maxS * tolerance) and (Hw1[-1] < Hw2[-1]): + print('Iteration : %d, Hw12 : %g Hw12/max_spec : %g' % + (ix, maxHw12, (maxHw12 / max_spec))) + if (maxHw12 < max_spec * tolerance) and (Hw1[-1] < Hw2[-1]): break # end Hw2 = Hw1 # end # Hw1(end) - # maxS*1e-3 - # if Hw1[-1]*S.data>maxS*1e-3, + # max_spec*1e-3 + # if Hw1[-1]*S.data>max_spec*1e-3, # warning('The Nyquist frequency of the spectrum may be too low') # end @@ -1837,7 +1837,7 @@ class SpecData1D(PlotData): # end try: - f.cl, f.pl = qlevels(f.f, [10, 30, 50, 70, 90, 95, 99, 99.9], + f.cl, f.pl = qlevels(f.data, [10, 30, 50, 70, 90, 95, 99, 99.9], f.args[0], f.args[1]) except Exception: warnings.warn('Singularity likely in pdf') diff --git a/wafo/stats/core.py b/wafo/stats/core.py index c1d31ad..b3e92cc 100644 --- a/wafo/stats/core.py +++ b/wafo/stats/core.py @@ -28,7 +28,6 @@ def now(): def valarray(shape, value=nan, typecode=None): """Return an array of all value. """ - #out = reshape(repeat([value], product(shape, axis=0), axis=0), shape) out = ones(shape, dtype=bool) * value if typecode is not None: out = out.astype(typecode) @@ -37,7 +36,7 @@ def valarray(shape, value=nan, typecode=None): return out -def _cdff(self, x, dfn, dfd): +def _cdff(x, dfn, dfd): return special.fdtr(dfn, dfd, x) @@ -91,17 +90,15 @@ def edf(x, method=2): z = atleast_1d(x) z.sort() - N = len(z) + n = len(z) if method == 1: - Fz1 = arange(0.5, N) / N + e_cdf = arange(0.5, n) / n elif method == 3: - Fz1 = arange(1, N + 1) / N + e_cdf = arange(1, n + 1) / n else: - Fz1 = arange(1, N + 1) / (N + 1) + e_cdf = arange(1, n + 1) / (n + 1) - F = PlotData(Fz1, z, xlab='x', ylab='F(x)') - F.setplotter('step') - return F + return PlotData(e_cdf, z, xlab='x', ylab='F(x)', plotmethod='step') def edfcnd(x, c=None, method=2): @@ -113,9 +110,9 @@ def edfcnd(x, c=None, method=2): x : array-like data vector method : integer scalar - 1. Interpolation so that F(X_(k)) == (k-0.5)/n. - 2. Interpolation so that F(X_(k)) == k/(n+1). (default) - 3. The empirical distribution. F(X_(k)) = k/n + 1. Interpolation so that cdf(X_(k)) == (k-0.5)/n. + 2. Interpolation so that cdf(X_(k)) == k/(n+1). (default) + 3. The empirical distribution. cdf(X_(k)) = k/n Example ------- @@ -124,8 +121,8 @@ def edfcnd(x, c=None, method=2): >>> R = ws.rayleigh.rvs(scale=2,size=100) >>> Fc = ws.edfcnd(R, 1) >>> hc = Fc.plot() - >>> F = ws.edf(R) - >>> h = F.plot() + >>> cdf = ws.edf(R) + >>> h = cdf.plot() See also edf, pdfplot, cumtrapz """ @@ -133,15 +130,13 @@ def edfcnd(x, c=None, method=2): if c is None: c = floor(min(z.min(), 0)) - try: - F = edf(z[c <= z], method=method) - except: - ValueError('No data points above c=%d' % int(c)) + cdf = edf(z[c <= z], method=method) if - inf < c: - F.labels.ylab = 'F(x| X>=%g)' % c + cdf.labels.ylab = 'F(x| X>=%g)' % c + + return cdf - return F def _check_nmin(nmin, n): nmin = max(nmin, 1) @@ -149,6 +144,7 @@ def _check_nmin(nmin, n): warnings.warn('nmin possibly too large!') return nmin + def _check_umin_umax(data, umin, umax, nmin): sd = np.sort(data) sdmax, sdmin = sd[-nmin], sd[0] @@ -156,11 +152,13 @@ def _check_umin_umax(data, umin, umax, nmin): umin = sdmin if umin is None else max(umin, sdmin) return umin, umax + def _check_nu(nu, nmin, n): if nu is None: nu = min(n - nmin, 100) return nu + def reslife(data, u=None, umin=None, umax=None, nu=None, nmin=3, alpha=0.05, plotflag=False): """ @@ -227,7 +225,9 @@ def reslife(data, u=None, umin=None, umax=None, nu=None, nmin=3, alpha=0.05, # srl = valarray(nu) # num = valarray(nu) - mean_and_std = lambda data1: (data1.mean(), data1.std(), data1.size) + def mean_and_std(data1): + return data1.mean(), data1.std(), data1.size + dat = arr(data) tmp = arr([mean_and_std(dat[dat > tresh] - tresh) for tresh in u.tolist()]) @@ -236,20 +236,18 @@ def reslife(data, u=None, umin=None, umax=None, nu=None, nmin=3, alpha=0.05, alpha2 = alpha / 2 # Approximate P% confidence interval - # Za = -invnorm(alpha2); % known mean - Za = -_invt(alpha2, num - 1) # unknown mean - mrlu = mrl + Za * srl / sqrt(num) - mrll = mrl - Za * srl / sqrt(num) - - # options.CI = [mrll,mrlu]; - # options.numdata = num; - titleTxt = 'Mean residual life with %d%s CI' % (100 * p, '%') + # z_a = -invnorm(alpha2); % known mean + z_a = -_invt(alpha2, num - 1) # unknown mean + mrlu = mrl + z_a * srl / sqrt(num) + mrll = mrl - z_a * srl / sqrt(num) + + title_txt = 'Mean residual life with %d%s CI' % (100 * p, '%') res = PlotData(mrl, u, xlab='Threshold', - ylab='Mean Excess', title=titleTxt) + ylab='Mean Excess', title=title_txt) res.workspace = dict( numdata=num, umin=umin, umax=umax, nu=nu, nmin=nmin, alpha=alpha) res.children = [ - PlotData(vstack([mrll, mrlu]).T, u, xlab='Threshold', title=titleTxt)] + PlotData(vstack([mrll, mrlu]).T, u, xlab='Threshold', title=title_txt)] res.plot_args_children = [':r'] if plotflag: res.plot() @@ -390,11 +388,11 @@ def dispersion_idx( b_u, ok_u = _find_appropriate_threshold(u, di, di_low, di_up) ci_txt = '{0:d}{1} CI'.format(100 * p, '%') - titleTxt = 'Dispersion Index plot' + title_txt = 'Dispersion Index plot' - res = PlotData(di, u, title=titleTxt, + res = PlotData(di, u, title=title_txt, labx='Threshold', laby='Dispersion Index') - #'caption',ci_txt); + res.workspace = dict(umin=umin, umax=umax, nu=nu, nmin=nmin, alpha=alpha) res.children = [ PlotData(vstack([di_low * ones(nu), di_up * ones(nu)]).T, u, @@ -448,7 +446,6 @@ def decluster(data, t=None, thresh=None, tmin=1): return data[i], t[i] - def _remove_index_to_data_too_close_to_each_other(ix_e, is_too_small, di_e, ti_e, tmin): is_too_close = np.hstack((is_too_small[0], is_too_small[:-1] | is_too_small[1:], is_too_small[-1])) @@ -759,18 +756,18 @@ class RegLogit(object): See also regglm, reglm, regnonlm """ - #% Original for MATLAB written by Gordon K Smyth , - #% U of Queensland, Australia, on Nov 19, 1990. Last revision Aug 3, - #% 1992. + # Original for MATLAB written by Gordon K Smyth , + # U of Queensland, Australia, on Nov 19, 1990. Last revision Aug 3, + # 1992. # - #% Author: Gordon K Smyth , - #% Revised by: pab - #% -renamed from oridinal to reglogit - #% -added predict, summary and compare - #% Description: Ordinal logistic regression + # Author: Gordon K Smyth , + # Revised by: pab + # -renamed from oridinal to reglogit + # -added predict, summary and compare + # Description: Ordinal logistic regression # - #% Uses the auxiliary functions logistic_regression_derivatives and - #% logistic_regression_likelihood. + # Uses the auxiliary functions logistic_regression_derivatives and + # logistic_regression_likelihood. def __init__(self, maxiter=500, accuracy=1e-6, alpha=0.05, deletecolinear=True, verbose=False): @@ -939,7 +936,7 @@ class RegLogit(object): break # end %while - #% tidy up output + # tidy up output theta = tb[:nz, ] beta = tb[nz:(nz + nx)] @@ -1204,7 +1201,7 @@ class RegLogit(object): R = np.dot(U, np.dot(np.diag(sqrt(S)), V)) ib = np.r_[0, nz:np1] - #% Var(eta_i) = var(theta_i+Xnew*b) + # Var(eta_i) = var(theta_i+Xnew*b) vareta = zeros((n, nz)) u = np.hstack((one, Xnew)) for i in range(nz): @@ -1317,7 +1314,7 @@ def test_reglogit(): x = np.arange(1, 11).reshape(-1, 1) b = RegLogit() b.fit(y, x) - # b.display() #% members and methods + # b.display() # members and methods b.summary() [mu, plo, pup] = b.predict(fulloutput=True) # @UnusedVariable @@ -1331,7 +1328,7 @@ def test_reglogit2(): y = (np.cos(x) > 2 * np.random.rand(n, 1) - 1) b = RegLogit() b.fit(y, x) - # b.display() #% members and methods + # b.display() # members and methods b.summary() [mu, plo, pup] = b.predict(fulloutput=True) import matplotlib.pyplot as pl