From 745237a4c039eb7ea84b6efb179949f3c2fa3f5a Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Sat, 1 Dec 2012 14:07:26 +0000 Subject: [PATCH] Ongoing work --- pywafo/src/wafo/bitwise.py | 7 + pywafo/src/wafo/covariance/core.py | 6 +- pywafo/src/wafo/interpolate.py | 6 +- pywafo/src/wafo/kdetools.py | 63 +-- pywafo/src/wafo/misc.py | 10 +- pywafo/src/wafo/objects.py | 338 +++++++-------- pywafo/src/wafo/polynomial.py | 2 +- pywafo/src/wafo/polynomial_old.py | 4 +- pywafo/src/wafo/spectrum/__init__.py | 2 +- pywafo/src/wafo/spectrum/core.py | 406 +++++++++--------- .../src/wafo/spectrum/dispersion_relation.py | 206 --------- .../spectrum/test/test_dispersion_relation.py | 4 +- pywafo/src/wafo/stats/__init__.py | 2 +- pywafo/src/wafo/stats/core.py | 44 +- pywafo/src/wafo/stats/distributions.py | 47 +- pywafo/src/wafo/stats/estimation.py | 27 +- pywafo/src/wafo/stats/misc.py | 2 +- pywafo/src/wafo/test/test_gaussian.py | 6 +- pywafo/src/wafo/test/test_kdetools.py | 10 +- pywafo/src/wafo/test/test_misc.py | 28 +- pywafo/src/wafo/transform/core.py | 9 +- pywafo/src/wafo/transform/models.py | 11 +- pywafo/src/wafo/wafodata.py | 24 +- pywafo/src/wafo/wave_theory/core.py | 170 ++++---- .../wafo/wave_theory/dispersion_relation.py | 21 +- 25 files changed, 638 insertions(+), 817 deletions(-) delete mode 100644 pywafo/src/wafo/spectrum/dispersion_relation.py diff --git a/pywafo/src/wafo/bitwise.py b/pywafo/src/wafo/bitwise.py index 225bcec..fd53dd9 100644 --- a/pywafo/src/wafo/bitwise.py +++ b/pywafo/src/wafo/bitwise.py @@ -66,6 +66,13 @@ def setbit(i, bit, value=1): def setbits(bitlist): """ Set bits of val to values in bitlist + + Example + ------- + >>> setbits([1,1]) + 3 + >>> setbits([1,0]) + 1 """ # return bitlist[7]<<7 | bitlist[6]<<6 | bitlist[5]<<5 | bitlist[4]<<4 | \ # bitlist[3]<<3 | bitlist[2]<<2 | bitlist[1]<<1 | bitlist[0] diff --git a/pywafo/src/wafo/covariance/core.py b/pywafo/src/wafo/covariance/core.py index 37bc56a..cd245ba 100644 --- a/pywafo/src/wafo/covariance/core.py +++ b/pywafo/src/wafo/covariance/core.py @@ -27,7 +27,7 @@ from scipy.linalg import toeplitz, sqrtm, svd, cholesky, diagsvd, pinv from scipy import sparse from pylab import stineman_interp -from wafo.wafodata import WafoData +from wafo.wafodata import PlotData from wafo.misc import sub_dict_select, nextpow2 #, JITImport import wafo.spectrum as _wafospec #_wafospec = JITImport('wafo.spectrum') @@ -103,7 +103,7 @@ def _set_seed(iseed): # return dot(std,random.randn(n, cases)) + mu[:,newaxis] -class CovData1D(WafoData): +class CovData1D(PlotData): """ Container class for 1D covariance data objects in WAFO Member variables @@ -127,7 +127,7 @@ class CovData1D(WafoData): See also -------- - WafoData + PlotData CovData """ diff --git a/pywafo/src/wafo/interpolate.py b/pywafo/src/wafo/interpolate.py index 2d4c977..e8e855e 100644 --- a/pywafo/src/wafo/interpolate.py +++ b/pywafo/src/wafo/interpolate.py @@ -52,7 +52,7 @@ def savitzky_golay(y, window_size, order, deriv=0): Notes ----- The Savitzky-Golay is a type of low-pass filter, particularly - suited for smoothing noisy data. The main idea behind this + suited for smoothing noisy data. The test_doctstrings idea behind this approach is to make for each point a least-square fit with a polynomial of high order over a odd-sized window centered at the point. @@ -1088,7 +1088,7 @@ def demo_monoticity(): plt.ioff() plt.show() -def main(): +def test_doctstrings(): from scipy import interpolate import matplotlib.pyplot as plt import matplotlib @@ -1150,7 +1150,7 @@ def test_docstrings(): if __name__ == '__main__': # test_docstrings() - #main() + #test_doctstrings() #test_smoothing_spline() #compare_methods() demo_monoticity() \ No newline at end of file diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py index 166a38c..bdd8868 100644 --- a/pywafo/src/wafo/kdetools.py +++ b/pywafo/src/wafo/kdetools.py @@ -21,11 +21,13 @@ from scipy.ndimage.morphology import distance_transform_edt from numpy import pi, sqrt, atleast_2d, exp, newaxis #@UnresolvedImport from wafo.misc import meshgrid, nextpow2, tranproc #, trangood -from wafo.wafodata import WafoData +from wafo.wafodata import PlotData from wafo.dctpack import dct, dctn, idctn from wafo.plotbackend import plotbackend as plt from wafo import fig +_TINY = np.finfo(float).machar.tiny + def _invnorm(q): return special.ndtri(q) @@ -65,6 +67,7 @@ def sphere_volume(d, r=1.0): Chapman and Hall, pp 105 """ return (r ** d) * 2.0 * pi ** (d / 2.0) / (d * gamma(d / 2.0)) + class KDEgauss(object): """ Kernel-Density Estimator base class. @@ -136,7 +139,7 @@ class KDEgauss(object): for i in ind.tolist(): # h[i] = get_smoothing(self.dataset[i]) deth = h.prod() - self.inv_hs = linalg.diag(1.0 / h) + self.inv_hs = np.diag(1.0 / h) else: #fully general smoothing matrix deth = linalg.det(h) if deth <= 0: @@ -836,7 +839,7 @@ class KDE(_KDE): for i in ind.tolist(): # h[i] = get_smoothing(self.dataset[i]) deth = h.prod() - self.inv_hs = linalg.diag(1.0 / h) + self.inv_hs = np.diag(1.0 / h) else: #fully general smoothing matrix deth = linalg.det(h) if deth <= 0: @@ -1025,7 +1028,7 @@ class KRegression(_KDE): >>> y = 2*np.exp(-x**2/(2*0.3**2))+3*np.exp(-(x-1)**2/(2*0.7**2)) + ei >>> kreg = wk.KRegression(x, y) >>> f = kreg(output='plotobj', title='Kernel regression', plotflag=1) - >>> f.plot(label='p=0') + >>> h = f.plot(label='p=0') """ def __init__(self, data, y, p=0, hs=None, kernel=None, alpha=0.0, xmin=None, xmax=None, inc=128, L2=None): @@ -3045,22 +3048,24 @@ def evar(y): Examples: -------- 1D signal - >>> n = 1e6; + >>> n = 1e6 >>> x = np.linspace(0,100,n); >>> y = np.cos(x/10)+(x/50) >>> var0 = 0.02 # noise variance >>> yn = y + sqrt(var0)*np.random.randn(*y.shape) - >>> evar(yn) #estimated variance - 0.020018619214933957 - + >>> s = evar(yn) #estimated variance + >>> np.abs(s-var0)/var0 < 3.5/np.sqrt(n) + True + 2D function >>> xp = np.linspace(0,1,50) >>> x, y = np.meshgrid(xp,xp) >>> f = np.exp(x+y) + np.sin((x-2*y)*3) >>> var0 = 0.04 # noise variance >>> fn = f + sqrt(var0)*np.random.randn(*f.shape) - >>> evar(fn) # estimated variance - 0.022806636928561375 + >>> s = evar(fn) # estimated variance + >>> np.abs(s-var0)/var0 < 3.5/np.sqrt(50) + True 3D function >>> yp = np.linspace(-2,2,50) @@ -3068,8 +3073,10 @@ def evar(y): >>> f = x*exp(-x**2-y**2-z**2) >>> var0 = 0.5 # noise variance >>> fn = f + sqrt(var0)*np.random.randn(*f.shape) - >>> evar(fn) # estimated variance - 0.47375136534336421 + >>> s = evar(fn) # estimated variance + >>> np.abs(s-var0)/var0 < 3.5/np.sqrt(50) + True + Other example ------------- @@ -3161,7 +3168,7 @@ def smoothn(data, s=None, weight=None, robust=False, z0=None, tolz=1e-3, maxiter 1-D example >>> import matplotlib.pyplot as plt >>> x = np.linspace(0,100,2**8) - >>> y = cos(x/10)+(x/50)**2 + np.random.randn(x.shape)/10 + >>> y = np.cos(x/10)+(x/50)**2 + np.random.randn(*x.shape)/10 >>> y[np.r_[70, 75, 80]] = np.array([5.5, 5, 6]) >>> z = smoothn(y) # Regular smoothing >>> zr = smoothn(y,robust=True) # Robust smoothing @@ -3839,8 +3846,8 @@ def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): # ref Casella and Berger (1990) "Statistical inference" pp444 - a = 2*pi + z0**2/(ciii+1e-16) - b = 2*(1+z0**2/(ciii+1e-16)) +# a = 2*pi + z0**2/(ciii+1e-16) +# b = 2*(1+z0**2/(ciii+1e-16)) # plo2 = ((a-sqrt(a**2-2*pi**2*b))/b).clip(min=0,max=1) # pup2 = ((a+sqrt(a**2-2*pi**2*b))/b).clip(min=0,max=1) @@ -3986,15 +3993,15 @@ def check_kreg_demo4(): #kde_gauss_demo() #kreg_demo2(n=120,symmetric=True,fun='hste', plotlog=True) k = 0 - for i, n in enumerate([100,300,600,4000]): + for i, n in enumerate([100,300,600,4000]): #@UnusedVariable x,y, fun1 = _get_data(n, symmetric=True,loc1=0.1, scale1=0.6, scale2=0.75) - k0 = k - hopt1, h1,h2 = _get_regression_smooting(x,y,fun='hos') - hopt2, h1,h2 = _get_regression_smooting(x,y,fun='hste') + #k0 = k + hopt1, _h1, _h2 = _get_regression_smooting(x,y,fun='hos') + hopt2, _h1, _h2 = _get_regression_smooting(x,y,fun='hste') hopt = sqrt(hopt1*hopt2) #hopt = _get_regression_smooting(x,y,fun='hos')[0] - for j, fun in enumerate(['hste']): # , 'hisj', 'hns', 'hstt' - hsmax, hs1, hs2 =_get_regression_smooting(x,y,fun=fun) + for j, fun in enumerate(['hste']): # , 'hisj', 'hns', 'hstt' @UnusedVariable + hsmax, _hs1, _hs2 =_get_regression_smooting(x,y,fun=fun) fmax = kreg_demo4(x, y, hsmax+0.1, hopt) for hi in np.linspace(hsmax*0.1,hsmax,55): @@ -4019,7 +4026,7 @@ def check_regression_bin(): #kde_gauss_demo() #kreg_demo2(n=120,symmetric=True,fun='hste', plotlog=True) k = 0 - for i, n in enumerate([100,300,600,4000]): + for i, n in enumerate([100,300,600,4000]): #@UnusedVariable x,y, fun1 = _get_data(n, symmetric=True,loc1=0.1, scale1=0.6, scale2=0.75) fbest = regressionbin(x, y, alpha=0.05, color='g', label='Transit_D') @@ -4040,7 +4047,7 @@ def check_regression_bin(): def check_bkregression(): plt.ion() k = 0 - for i, n in enumerate([50, 100,300,600]): + for i, n in enumerate([50, 100,300,600]): #@UnusedVariable x,y, fun1 = _get_data(n, symmetric=True,loc1=0.1, scale1=0.6, scale2=0.75) bkreg = BKRegression(x,y) fbest = bkreg.prb_search_best(hsfun='hste', alpha=0.05, color='g', label='Transit_D') @@ -4204,7 +4211,7 @@ def regressionbin(x,y, alpha=0.05, color='r', label=''): f = smoothed_bin_prb(x, y, hi, hopt, alpha, color, label, bin_prb) if f.aicc<=fbest.aicc: fbest = f - hbest = hi + #hbest = hi return fbest def kde_gauss_demo(n=50): ''' @@ -4227,7 +4234,7 @@ def kde_gauss_demo(n=50): #dist = st.norm dist = st.expon data = dist.rvs(loc=0, scale=1.0, size=n) - d, N = np.atleast_2d(data).shape + d, N = np.atleast_2d(data).shape #@UnusedVariable if d==1: plot_options = [dict(color='red'), dict(color='green'), dict(color='black')] @@ -4265,7 +4272,9 @@ def test_docstrings(): doctest.testmod() if __name__ == '__main__': - check_bkregression() + test_docstrings() + + #check_bkregression() #check_regression_bin() #check_kreg_demo3() #check_kreg_demo4() @@ -4274,7 +4283,7 @@ if __name__ == '__main__': #test_smoothn_2d() #test_smoothn_cardioid() - #test_docstrings() + #kde_demo2() #kreg_demo1(fast=True) #kde_gauss_demo() diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index 2a4a91a..e833c98 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -30,7 +30,7 @@ __all__ = ['is_numlike', 'JITImport', 'DotDict', 'Bunch', 'printf', 'sub_dict_se 'findextrema', 'findpeaks', 'findrfc', 'rfcfilter', 'findtp', 'findtc', 'findoutliers', 'common_shape', 'argsreduce', 'stirlerr', 'getshipchar', 'betaloge', 'gravity', 'nextpow2', - 'discretize', 'pol2cart', 'cart2pol', 'meshgrid', 'ndgrid', + 'discretize', 'polar2cart', 'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'plot_histgrm', 'num2pistr', 'test_docstrings'] @@ -1770,7 +1770,7 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5): return x, fx -def pol2cart(theta, rho, z=None): +def polar2cart(theta, rho, z=None): ''' Transform polar coordinates into 2D cartesian coordinates. @@ -1781,7 +1781,7 @@ def pol2cart(theta, rho, z=None): See also -------- - cart2pol + cart2polar ''' x, y = rho * cos(theta), rho * sin(theta) if z is None: @@ -1790,7 +1790,7 @@ def pol2cart(theta, rho, z=None): return x, y, z -def cart2pol(x, y, z=None): +def cart2polar(x, y, z=None): ''' Transform 2D cartesian coordinates into polar coordinates. Returns @@ -1802,7 +1802,7 @@ def cart2pol(x, y, z=None): See also -------- - pol2cart + polar2cart ''' t, r = arctan2(y, x), hypot(x, y) if z is None: diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index 8d82c3a..0f37d69 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -18,7 +18,7 @@ from wafo.transform.models import TrHermite, TrOchi, TrLinear from wafo.stats import edf, distributions from wafo.misc import (nextpow2, findtp, findrfc, findtc, findcross, ecross, JITImport, DotDict, gravity, findrfc_astm) -from wafodata import WafoData +from wafodata import PlotData from wafo.interpolate import SmoothSpline from scipy.interpolate.interpolate import interp1d from scipy.integrate.quadrature import cumtrapz #@UnresolvedImport @@ -56,7 +56,7 @@ __all__ = ['TimeSeries', 'LevelCrossings', 'CyclePairs', 'TurningPoints', def _invchi2(q, df): return special.chdtri(df, q) -class LevelCrossings(WafoData): +class LevelCrossings(PlotData): ''' Container class for Level crossing data objects in WAFO @@ -109,7 +109,7 @@ class LevelCrossings(WafoData): cmax = self.data[icmax] x = (self.args - self.mean) / self.sigma y = cmax * exp(-x ** 2 / 2.0) - self.children = [WafoData(y, self.args)] + self.children = [PlotData(y, self.args)] def extrapolate(self, u_min=None, u_max=None, method='ml', dist='genpar', plotflag=0): ''' @@ -663,7 +663,7 @@ def test_levelcrossings_extrapolate(): s = x[:,1].std() lc_gpd = lc.extrapolate(-2*s, 2*s, dist='rayleigh') #@UnusedVariable -class CyclePairs(WafoData): +class CyclePairs(PlotData): ''' Container class for Cycle Pairs data objects in WAFO @@ -839,7 +839,7 @@ class CyclePairs(WafoData): ylab = 'Intensity [count/sec]' return LevelCrossings(dcount, levels, mean=self.mean, sigma=self.sigma, ylab=ylab, intensity=intensity) -class TurningPoints(WafoData): +class TurningPoints(PlotData): ''' Container class for Turning Points data objects in WAFO @@ -1029,7 +1029,7 @@ def mat2timeseries(x): """ return TimeSeries(x[:, 1::], x[:, 0].ravel()) -class TimeSeries(WafoData): +class TimeSeries(PlotData): ''' Container class for 1D TimeSeries data objects in WAFO @@ -1177,7 +1177,7 @@ class TimeSeries(WafoData): #cumsum = np.cumsum acf = _wafocov.CovData1D(R[lags], t) acf.sigma = sqrt(r_[ 0, r0 ** 2 , r0 ** 2 + 2 * cumsum(R[1:] ** 2)] / Ncens) - acf.children = [WafoData(-2. * acf.sigma[lags], t), WafoData(2. * acf.sigma[lags], t)] + acf.children = [PlotData(-2. * acf.sigma[lags], t), PlotData(2. * acf.sigma[lags], t)] acf.plot_args_children = ['r:'] acf.norm = norm return acf @@ -2244,167 +2244,167 @@ class TimeSeries(WafoData): #wafostamp return figs -def hyperbolic_ratio(a, b, sa, sb): - ''' - Return ratio of hyperbolic functions - to allow extreme variations of arguments. - - Parameters - ---------- - a, b : array-like - arguments vectors of the same size - sa, sb : scalar integers - defining the hyperbolic function used, i.e., f(x,1)=cosh(x), f(x,-1)=sinh(x) - - Returns - ------- - r : ndarray - f(a,sa)/f(b,sb), ratio of hyperbolic functions of same - size as a and b - Examples - -------- - >>> x = [-2,0,2] - >>> hyperbolic_ratio(x,1,1,1) # gives r=cosh(x)/cosh(1) - array([ 2.438107 , 0.64805427, 2.438107 ]) - >>> hyperbolic_ratio(x,1,1,-1) # gives r=cosh(x)/sinh(1) - array([ 3.20132052, 0.85091813, 3.20132052]) - >>> hyperbolic_ratio(x,1,-1,1) # gives r=sinh(x)/cosh(1) - array([-2.35040239, 0. , 2.35040239]) - >>> hyperbolic_ratio(x,1,-1,-1) # gives r=sinh(x)/sinh(1) - array([-3.08616127, 0. , 3.08616127]) - >>> hyperbolic_ratio(1,x,1,1) # gives r=cosh(1)/cosh(x) - array([ 0.41015427, 1.54308063, 0.41015427]) - >>> hyperbolic_ratio(1,x,1,-1) # gives r=cosh(1)/sinh(x) - array([-0.42545906, inf, 0.42545906]) - >>> hyperbolic_ratio(1,x,-1,1) # gives r=sinh(1)/cosh(x) - array([ 0.3123711 , 1.17520119, 0.3123711 ]) - >>> hyperbolic_ratio(1,x,-1,-1) # gives r=sinh(1)/sinh(x) - array([-0.32402714, inf, 0.32402714]) - - See also - -------- - tran - ''' - ak, bk, sak, sbk = np.atleast_1d(a, b, sign(sa), sign(sb)) - # old call - #return exp(ak-bk)*(1+sak*exp(-2*ak))/(1+sbk*exp(-2*bk)) - # TODO: Does not always handle division by zero correctly - - signRatio = np.where(sak * ak < 0, sak, 1) - signRatio = np.where(sbk * bk < 0, sbk * signRatio, signRatio) - - bk = np.abs(bk) - ak = np.abs(ak) - - num = np.where(sak < 0, expm1(-2 * ak), 1 + exp(-2 * ak)) - den = np.where(sbk < 0, expm1(-2 * bk), 1 + exp(-2 * bk)) - iden = np.ones(den.shape) * inf - ind = np.flatnonzero(den != 0) - iden.flat[ind] = 1.0 / den[ind] - val = np.where(num == den, 1, num * iden) - return signRatio * exp(ak - bk) * val #((sak+exp(-2*ak))/(sbk+exp(-2*bk))) - -def sensor_typeid(*sensortypes): - ''' Return ID for sensortype name - - Parameter - --------- - sensortypes : list of strings defining the sensortype - - Returns - ------- - sensorids : list of integers defining the sensortype - - Valid senor-ids and -types for time series are as follows: - 0, 'n' : Surface elevation (n=Eta) - 1, 'n_t' : Vertical surface velocity - 2, 'n_tt' : Vertical surface acceleration - 3, 'n_x' : Surface slope in x-direction - 4, 'n_y' : Surface slope in y-direction - 5, 'n_xx' : Surface curvature in x-direction - 6, 'n_yy' : Surface curvature in y-direction - 7, 'n_xy' : Surface curvature in xy-direction - 8, 'P' : Pressure fluctuation about static MWL pressure - 9, 'U' : Water particle velocity in x-direction - 10, 'V' : Water particle velocity in y-direction - 11, 'W' : Water particle velocity in z-direction - 12, 'U_t' : Water particle acceleration in x-direction - 13, 'V_t' : Water particle acceleration in y-direction - 14, 'W_t' : Water particle acceleration in z-direction - 15, 'X_p' : Water particle displacement in x-direction from its mean position - 16, 'Y_p' : Water particle displacement in y-direction from its mean position - 17, 'Z_p' : Water particle displacement in z-direction from its mean position - - Example: - >>> sensor_typeid('W','v') - [11, 10] - >>> sensor_typeid('rubbish') - [nan] - - See also - -------- - sensor_type - ''' - - sensorid_table = dict(n=0, n_t=1, n_tt=2, n_x=3, n_y=4, n_xx=5, - n_yy=6, n_xy=7, p=8, u=9, v=10, w=11, u_t=12, - v_t=13, w_t=14, x_p=15, y_p=16, z_p=17) - try: - return [sensorid_table.get(name.lower(), nan) for name in sensortypes] - except: - raise ValueError('Input must be a string!') - - - -def sensor_type(*sensorids): - ''' - Return sensortype name - - Parameter - --------- - sensorids : vector or list of integers defining the sensortype - - Returns - ------- - sensornames : tuple of strings defining the sensortype - Valid senor-ids and -types for time series are as follows: - 0, 'n' : Surface elevation (n=Eta) - 1, 'n_t' : Vertical surface velocity - 2, 'n_tt' : Vertical surface acceleration - 3, 'n_x' : Surface slope in x-direction - 4, 'n_y' : Surface slope in y-direction - 5, 'n_xx' : Surface curvature in x-direction - 6, 'n_yy' : Surface curvature in y-direction - 7, 'n_xy' : Surface curvature in xy-direction - 8, 'P' : Pressure fluctuation about static MWL pressure - 9, 'U' : Water particle velocity in x-direction - 10, 'V' : Water particle velocity in y-direction - 11, 'W' : Water particle velocity in z-direction - 12, 'U_t' : Water particle acceleration in x-direction - 13, 'V_t' : Water particle acceleration in y-direction - 14, 'W_t' : Water particle acceleration in z-direction - 15, 'X_p' : Water particle displacement in x-direction from its mean position - 16, 'Y_p' : Water particle displacement in y-direction from its mean position - 17, 'Z_p' : Water particle displacement in z-direction from its mean position - - Example: - >>> sensor_type(range(3)) - ('n', 'n_t', 'n_tt') - - See also - -------- - sensor_typeid, tran - ''' - valid_names = ('n', 'n_t', 'n_tt', 'n_x', 'n_y', 'n_xx', 'n_yy', 'n_xy', - 'p', 'u', 'v', 'w', 'u_t', 'v_t', 'w_t', 'x_p', 'y_p', 'z_p', - nan) - ids = atleast_1d(*sensorids) - if isinstance(ids, list): - ids = hstack(ids) - n = len(valid_names) - 1 - ids = where(((ids < 0) | (n < ids)), n , ids) - return tuple(valid_names[i] for i in ids) - +#def hyperbolic_ratio(a, b, sa, sb): +# ''' +# Return ratio of hyperbolic functions +# to allow extreme variations of arguments. +# +# Parameters +# ---------- +# a, b : array-like +# arguments vectors of the same size +# sa, sb : scalar integers +# defining the hyperbolic function used, i.e., f(x,1)=cosh(x), f(x,-1)=sinh(x) +# +# Returns +# ------- +# r : ndarray +# f(a,sa)/f(b,sb), ratio of hyperbolic functions of same +# size as a and b +# Examples +# -------- +# >>> x = [-2,0,2] +# >>> hyperbolic_ratio(x,1,1,1) # gives r=cosh(x)/cosh(1) +# array([ 2.438107 , 0.64805427, 2.438107 ]) +# >>> hyperbolic_ratio(x,1,1,-1) # gives r=cosh(x)/sinh(1) +# array([ 3.20132052, 0.85091813, 3.20132052]) +# >>> hyperbolic_ratio(x,1,-1,1) # gives r=sinh(x)/cosh(1) +# array([-2.35040239, 0. , 2.35040239]) +# >>> hyperbolic_ratio(x,1,-1,-1) # gives r=sinh(x)/sinh(1) +# array([-3.08616127, 0. , 3.08616127]) +# >>> hyperbolic_ratio(1,x,1,1) # gives r=cosh(1)/cosh(x) +# array([ 0.41015427, 1.54308063, 0.41015427]) +# >>> hyperbolic_ratio(1,x,1,-1) # gives r=cosh(1)/sinh(x) +# array([-0.42545906, inf, 0.42545906]) +# >>> hyperbolic_ratio(1,x,-1,1) # gives r=sinh(1)/cosh(x) +# array([ 0.3123711 , 1.17520119, 0.3123711 ]) +# >>> hyperbolic_ratio(1,x,-1,-1) # gives r=sinh(1)/sinh(x) +# array([-0.32402714, inf, 0.32402714]) +# +# See also +# -------- +# tran +# ''' +# ak, bk, sak, sbk = np.atleast_1d(a, b, sign(sa), sign(sb)) +# # old call +# #return exp(ak-bk)*(1+sak*exp(-2*ak))/(1+sbk*exp(-2*bk)) +# # TODO: Does not always handle division by zero correctly +# +# signRatio = np.where(sak * ak < 0, sak, 1) +# signRatio = np.where(sbk * bk < 0, sbk * signRatio, signRatio) +# +# bk = np.abs(bk) +# ak = np.abs(ak) +# +# num = np.where(sak < 0, expm1(-2 * ak), 1 + exp(-2 * ak)) +# den = np.where(sbk < 0, expm1(-2 * bk), 1 + exp(-2 * bk)) +# iden = np.ones(den.shape) * inf +# ind = np.flatnonzero(den != 0) +# iden.flat[ind] = 1.0 / den[ind] +# val = np.where(num == den, 1, num * iden) +# return signRatio * exp(ak - bk) * val #((sak+exp(-2*ak))/(sbk+exp(-2*bk))) +# +#def sensor_typeid(*sensortypes): +# ''' Return ID for sensortype name +# +# Parameter +# --------- +# sensortypes : list of strings defining the sensortype +# +# Returns +# ------- +# sensorids : list of integers defining the sensortype +# +# Valid senor-ids and -types for time series are as follows: +# 0, 'n' : Surface elevation (n=Eta) +# 1, 'n_t' : Vertical surface velocity +# 2, 'n_tt' : Vertical surface acceleration +# 3, 'n_x' : Surface slope in x-direction +# 4, 'n_y' : Surface slope in y-direction +# 5, 'n_xx' : Surface curvature in x-direction +# 6, 'n_yy' : Surface curvature in y-direction +# 7, 'n_xy' : Surface curvature in xy-direction +# 8, 'P' : Pressure fluctuation about static MWL pressure +# 9, 'U' : Water particle velocity in x-direction +# 10, 'V' : Water particle velocity in y-direction +# 11, 'W' : Water particle velocity in z-direction +# 12, 'U_t' : Water particle acceleration in x-direction +# 13, 'V_t' : Water particle acceleration in y-direction +# 14, 'W_t' : Water particle acceleration in z-direction +# 15, 'X_p' : Water particle displacement in x-direction from its mean position +# 16, 'Y_p' : Water particle displacement in y-direction from its mean position +# 17, 'Z_p' : Water particle displacement in z-direction from its mean position +# +# Example: +# >>> sensor_typeid('W','v') +# [11, 10] +# >>> sensor_typeid('rubbish') +# [nan] +# +# See also +# -------- +# sensor_type +# ''' +# +# sensorid_table = dict(n=0, n_t=1, n_tt=2, n_x=3, n_y=4, n_xx=5, +# n_yy=6, n_xy=7, p=8, u=9, v=10, w=11, u_t=12, +# v_t=13, w_t=14, x_p=15, y_p=16, z_p=17) +# try: +# return [sensorid_table.get(name.lower(), nan) for name in sensortypes] +# except: +# raise ValueError('Input must be a string!') +# +# +# +#def sensor_type(*sensorids): +# ''' +# Return sensortype name +# +# Parameter +# --------- +# sensorids : vector or list of integers defining the sensortype +# +# Returns +# ------- +# sensornames : tuple of strings defining the sensortype +# Valid senor-ids and -types for time series are as follows: +# 0, 'n' : Surface elevation (n=Eta) +# 1, 'n_t' : Vertical surface velocity +# 2, 'n_tt' : Vertical surface acceleration +# 3, 'n_x' : Surface slope in x-direction +# 4, 'n_y' : Surface slope in y-direction +# 5, 'n_xx' : Surface curvature in x-direction +# 6, 'n_yy' : Surface curvature in y-direction +# 7, 'n_xy' : Surface curvature in xy-direction +# 8, 'P' : Pressure fluctuation about static MWL pressure +# 9, 'U' : Water particle velocity in x-direction +# 10, 'V' : Water particle velocity in y-direction +# 11, 'W' : Water particle velocity in z-direction +# 12, 'U_t' : Water particle acceleration in x-direction +# 13, 'V_t' : Water particle acceleration in y-direction +# 14, 'W_t' : Water particle acceleration in z-direction +# 15, 'X_p' : Water particle displacement in x-direction from its mean position +# 16, 'Y_p' : Water particle displacement in y-direction from its mean position +# 17, 'Z_p' : Water particle displacement in z-direction from its mean position +# +# Example: +# >>> sensor_type(range(3)) +# ('n', 'n_t', 'n_tt') +# +# See also +# -------- +# sensor_typeid, tran +# ''' +# valid_names = ('n', 'n_t', 'n_tt', 'n_x', 'n_y', 'n_xx', 'n_yy', 'n_xy', +# 'p', 'u', 'v', 'w', 'u_t', 'v_t', 'w_t', 'x_p', 'y_p', 'z_p', +# nan) +# ids = atleast_1d(*sensorids) +# if isinstance(ids, list): +# ids = hstack(ids) +# n = len(valid_names) - 1 +# ids = where(((ids < 0) | (n < ids)), n , ids) +# return tuple(valid_names[i] for i in ids) +# #class TransferFunction(object): # ''' # Class for computing transfer functions based on linear wave theory @@ -2822,7 +2822,7 @@ def main(): x = np.arange(-2, 2, 0.2) # Plot 2 objects in one call - d2 = WafoData(np.sin(x), x, xlab='x', ylab='sin', title='sinus') + d2 = PlotData(np.sin(x), x, xlab='x', ylab='sin', title='sinus') d0 = d2.copy() diff --git a/pywafo/src/wafo/polynomial.py b/pywafo/src/wafo/polynomial.py index c836ce9..b483d5f 100644 --- a/pywafo/src/wafo/polynomial.py +++ b/pywafo/src/wafo/polynomial.py @@ -21,7 +21,7 @@ from plotbackend import plotbackend as plt import numpy as np from numpy.fft import fft, ifft from numpy import (zeros, ones, zeros_like, array, asarray, newaxis, arange, #@UnresolvedImport - logical_or, any, pi, cos, round, diff, all, r_, exp, #atleast_1d, hstack,#@UnresolvedImport + logical_or, any, pi, cos, round, diff, all, r_, exp, atleast_1d, # hstack,#@UnresolvedImport where, extract, dot, linalg, sign, concatenate, floor, isreal, conj, remainder, #@UnresolvedImport linspace) #@UnresolvedImport from numpy.lib.polynomial import * #@UnusedWildImport diff --git a/pywafo/src/wafo/polynomial_old.py b/pywafo/src/wafo/polynomial_old.py index fbe9c70..4f78206 100644 --- a/pywafo/src/wafo/polynomial_old.py +++ b/pywafo/src/wafo/polynomial_old.py @@ -1171,7 +1171,7 @@ class Cheb1d(object): """ return Cheb1d(chebder(self.coeffs,self.a,self.b)) -def main(): +def test_doctstrings(): if False: #True: # x = np.arange(4) dx = dct(x) @@ -1204,4 +1204,4 @@ def main(): import doctest doctest.testmod() if __name__== '__main__': - main() + test_doctstrings() diff --git a/pywafo/src/wafo/spectrum/__init__.py b/pywafo/src/wafo/spectrum/__init__.py index 13b56b1..307e615 100644 --- a/pywafo/src/wafo/spectrum/__init__.py +++ b/pywafo/src/wafo/spectrum/__init__.py @@ -8,4 +8,4 @@ Spectrum package in WAFO Toolbox. from core import * #SpecData1D, SpecData2D, cltext import models -import dispersion_relation +from wafo.wave_theory import dispersion_relation diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py index 3c39661..f32bfd6 100644 --- a/pywafo/src/wafo/spectrum/core.py +++ b/pywafo/src/wafo/spectrum/core.py @@ -1,15 +1,15 @@ from __future__ import division -from wafo.misc import meshgrid, gravity, cart2pol, pol2cart -from wafo.objects import mat2timeseries, TimeSeries +from wafo.misc import meshgrid, gravity, cart2polar, polar2cart +from wafo.objects import TimeSeries #mat2timeseries, import warnings import numpy as np from numpy import (pi, inf, zeros, ones, where, nonzero, #@UnresolvedImport flatnonzero, ceil, sqrt, exp, log, arctan2, log10, #@UnresolvedImport - tanh, cosh, sinh, random, atleast_1d, maximum, #@UnresolvedImport + tanh, cosh, sinh, random, atleast_1d, #@UnresolvedImport minimum, diff, isnan, any, r_, conj, mod, #@UnresolvedImport hstack, vstack, interp, ravel, finfo, linspace, #@UnresolvedImport - arange, array, nan, newaxis, fliplr, sign) #@UnresolvedImport + arange, array, nan, newaxis, sign) #, fliplr, maximum) #@UnresolvedImport from numpy.fft import fft from scipy.integrate import simps, trapz from scipy.special import erf @@ -17,8 +17,8 @@ from scipy.linalg import toeplitz import scipy.interpolate as interpolate from wafo.interpolate import stineman_interp -from dispersion_relation import w2k #, k2w -from wafo.wafodata import WafoData, now +from wafo.wave_theory.dispersion_relation import w2k #, k2w +from wafo.wafodata import PlotData, now from wafo.misc import sub_dict_select, nextpow2, discretize, JITImport, findpeaks #, tranproc from wafo.graphutil import cltext from wafo.kdetools import qlevels @@ -45,7 +45,7 @@ from wafo.plotbackend import plotbackend # Trick to avoid error due to circular import -#_WAFOCOV = ppimport('wafo.covariance') + _WAFOCOV = JITImport('wafo.covariance') @@ -156,185 +156,185 @@ def qtf(w, h=inf, g=9.81): return h_s, h_d , h_dii -def plotspec(specdata, linetype='b-', flag=1): - ''' - PLOTSPEC Plot a spectral density - - Parameters - ---------- - S : SpecData1D or SpecData2D object - defining spectral density. - linetype : string - defining color and linetype, see plot for possibilities - flag : scalar integer - defining the type of plot - 1D: - 1 plots the density, S, (default) - 2 plot 10log10(S) - 3 plots both the above plots - 2D: - Directional spectra: S(w,theta), S(f,theta) - 1 polar plot S (default) - 2 plots spectral density and the directional - spreading, int S(w,theta) dw or int S(f,theta) df - 3 plots spectral density and the directional - spreading, int S(w,theta)/S(w) dw or int S(f,theta)/S(f) df - 4 mesh of S - 5 mesh of S in polar coordinates - 6 contour plot of S - 7 filled contour plot of S - Wavenumber spectra: S(k1,k2) - 1 contour plot of S (default) - 2 filled contour plot of S - - Example - ------- - >>> import numpy as np - >>> import wafo.spectrum as ws - >>> Sj = ws.models.Jonswap(Hm0=3, Tp=7) - >>> S = Sj.tospecdata() - >>> ws.plotspec(S,flag=1) - - S = demospec('dir'); S2 = mkdspec(jonswap,spreading); - plotspec(S,2), hold on - plotspec(S,3,'g') % Same as previous fig. due to frequency independent spreading - plotspec(S2,2,'r') % Not the same as previous figs. due to frequency dependent spreading - plotspec(S2,3,'m') - % transform from angular frequency and radians to frequency and degrees - Sf = ttspec(S,'f','d'); clf - plotspec(Sf,2), - - See also dat2spec, createspec, simpson - ''' - - # label the contour levels - txtFlag = 0; - LegendOn = 1; - - - ftype = specdata.freqtype #options are 'f' and 'w' and 'k' - data = specdata.data - if data.ndim == 2: - freq = specdata.args[1] - theta = specdata.args[0] - else: - freq = specdata.args - #if isinstance(specdata.args, (list, tuple)): - - if ftype == 'w': - xlbl_txt = 'Frequency [rad/s]'; - ylbl1_txt = 'S(w) [m^2 s / rad]'; - ylbl3_txt = 'Directional Spectrum'; - zlbl_txt = 'S(w,\theta) [m^2 s / rad^2]'; - funit = ' [rad/s]'; - #Sunit = ' [m^2 s / rad]'; - elif ftype == 'f': - xlbl_txt = 'Frequency [Hz]'; - ylbl1_txt = 'S(f) [m^2 s]'; - ylbl3_txt = 'Directional Spectrum'; - zlbl_txt = 'S(f,\theta) [m^2 s / rad]'; - funit = ' [Hz]'; - #Sunit = ' [m^2 s ]'; - elif ftype == 'k': - xlbl_txt = 'Wave number [rad/m]'; - ylbl1_txt = 'S(k) [m^3/ rad]'; - funit = ' [rad/m]'; - #Sunit = ' [m^3 / rad]'; - ylbl4_txt = 'Wave Number Spectrum'; - - else: - raise ValueError('Frequency type unknown') - - - if hasattr(specdata, 'norm') and specdata.norm : - #Sunit=[]; - funit = []; - ylbl1_txt = 'Normalized Spectral density'; - ylbl3_txt = 'Normalized Directional Spectrum'; - ylbl4_txt = 'Normalized Wave Number Spectrum'; - if ftype == 'k': - xlbl_txt = 'Normalized Wave number'; - else: - xlbl_txt = 'Normalized Frequency'; - - ylbl2_txt = 'Power spectrum (dB)'; - - phi = specdata.phi - - spectype = specdata.type.lower() - stype = spectype[-3::] - if stype in ('enc', 'req', 'k1d') : #1D plot - Fn = freq[-1] # Nyquist frequency - indm = findpeaks(data, n=4) - maxS = data.max() -# if isfield(S,'CI') && ~isempty(S.CI), -# maxS = maxS*S.CI(2); -# txtCI = [num2str(100*S.p), '% CI']; -# end - - Fp = freq[indm]# %peak frequency/wave number - - if len(indm) == 1: - txt = [('fp = %0.2g' % Fp) + funit] - else: - txt = [] - for i, fp in enumerate(Fp.tolist()): - txt.append(('fp%d = %0.2g' % (i, fp)) + funit) - - txt = ''.join(txt) - if (flag == 3): - plotbackend.subplot(2, 1, 1) - if (flag == 1) or (flag == 3):#% Plot in normal scale - plotbackend.plot(np.vstack([Fp, Fp]), np.vstack([zeros(len(indm)), data.take(indm)]), ':', label=txt) - plotbackend.plot(freq, data, linetype) - specdata.labels.labelfig() -# if isfield(S,'CI'), -# plot(freq,S.S*S.CI(1), 'r:' ) -# plot(freq,S.S*S.CI(2), 'r:' ) - - a = plotbackend.axis() - - a1 = Fn - if (Fp > 0): - a1 = max(min(Fn, 10 * max(Fp)), a[1]); - - plotbackend.axis([0, a1 , 0, max(1.01 * maxS, a[3])]) - #plotbackend.title('Spectral density') - #plotbackend.xlabel(xlbl_txt) - #plotbackend.ylabel(ylbl1_txt) - - - if (flag == 3): - plotbackend.subplot(2, 1, 2) - - if (flag == 2) or (flag == 3) : # Plot in logaritmic scale - ind = np.flatnonzero(data > 0) - - plotbackend.plot(np.vstack([Fp, Fp]), - np.vstack((min(10 * log10(data.take(ind) / maxS)).repeat(len(Fp)), - 10 * log10(data.take(indm) / maxS))), ':',label=txt) -# hold on -# if isfield(S,'CI'), -# plot(freq(ind),10*log10(S.S(ind)*S.CI(1)/maxS), 'r:' ) -# plot(freq(ind),10*log10(S.S(ind)*S.CI(2)/maxS), 'r:' ) -# end - plotbackend.plot(freq[ind], 10 * log10(data[ind] / maxS), linetype) - - a = plotbackend.axis() - - a1 = Fn - if (Fp > 0): - a1 = max(min(Fn, 10 * max(Fp)), a[1]); - - plotbackend.axis([0, a1 , -20, max(1.01 * 10 * log10(1), a[3])]) - - specdata.labels.labelfig() - #plotbackend.title('Spectral density') - #plotbackend.xlabel(xlbl_txt) - plotbackend.ylabel(ylbl2_txt) - - if LegendOn: - plotbackend.legend() +#def plotspec(specdata, linetype='b-', flag=1): +# ''' +# PLOTSPEC Plot a spectral density +# +# Parameters +# ---------- +# S : SpecData1D or SpecData2D object +# defining spectral density. +# linetype : string +# defining color and linetype, see plot for possibilities +# flag : scalar integer +# defining the type of plot +# 1D: +# 1 plots the density, S, (default) +# 2 plot 10log10(S) +# 3 plots both the above plots +# 2D: +# Directional spectra: S(w,theta), S(f,theta) +# 1 polar plot S (default) +# 2 plots spectral density and the directional +# spreading, int S(w,theta) dw or int S(f,theta) df +# 3 plots spectral density and the directional +# spreading, int S(w,theta)/S(w) dw or int S(f,theta)/S(f) df +# 4 mesh of S +# 5 mesh of S in polar coordinates +# 6 contour plot of S +# 7 filled contour plot of S +# Wavenumber spectra: S(k1,k2) +# 1 contour plot of S (default) +# 2 filled contour plot of S +# +# Example +# ------- +# >>> import numpy as np +# >>> import wafo.spectrum as ws +# >>> Sj = ws.models.Jonswap(Hm0=3, Tp=7) +# >>> S = Sj.tospecdata() +# >>> ws.plotspec(S,flag=1) +# +# S = demospec('dir'); S2 = mkdspec(jonswap,spreading); +# plotspec(S,2), hold on +# plotspec(S,3,'g') % Same as previous fig. due to frequency independent spreading +# plotspec(S2,2,'r') % Not the same as previous figs. due to frequency dependent spreading +# plotspec(S2,3,'m') +# % transform from angular frequency and radians to frequency and degrees +# Sf = ttspec(S,'f','d'); clf +# plotspec(Sf,2), +# +# See also dat2spec, createspec, simpson +# ''' +# +# # label the contour levels +# txtFlag = 0; +# LegendOn = 1; +# +# +# ftype = specdata.freqtype #options are 'f' and 'w' and 'k' +# data = specdata.data +# if data.ndim == 2: +# freq = specdata.args[1] +# theta = specdata.args[0] +# else: +# freq = specdata.args +# #if isinstance(specdata.args, (list, tuple)): +# +# if ftype == 'w': +# xlbl_txt = 'Frequency [rad/s]'; +# ylbl1_txt = 'S(w) [m^2 s / rad]'; +# ylbl3_txt = 'Directional Spectrum'; +# zlbl_txt = 'S(w,\theta) [m^2 s / rad^2]'; +# funit = ' [rad/s]'; +# #Sunit = ' [m^2 s / rad]'; +# elif ftype == 'f': +# xlbl_txt = 'Frequency [Hz]'; +# ylbl1_txt = 'S(f) [m^2 s]'; +# ylbl3_txt = 'Directional Spectrum'; +# zlbl_txt = 'S(f,\theta) [m^2 s / rad]'; +# funit = ' [Hz]'; +# #Sunit = ' [m^2 s ]'; +# elif ftype == 'k': +# xlbl_txt = 'Wave number [rad/m]'; +# ylbl1_txt = 'S(k) [m^3/ rad]'; +# funit = ' [rad/m]'; +# #Sunit = ' [m^3 / rad]'; +# ylbl4_txt = 'Wave Number Spectrum'; +# +# else: +# raise ValueError('Frequency type unknown') +# +# +# if hasattr(specdata, 'norm') and specdata.norm : +# #Sunit=[]; +# funit = []; +# ylbl1_txt = 'Normalized Spectral density'; +# ylbl3_txt = 'Normalized Directional Spectrum'; +# ylbl4_txt = 'Normalized Wave Number Spectrum'; +# if ftype == 'k': +# xlbl_txt = 'Normalized Wave number'; +# else: +# xlbl_txt = 'Normalized Frequency'; +# +# ylbl2_txt = 'Power spectrum (dB)'; +# +# phi = specdata.phi +# +# spectype = specdata.type.lower() +# stype = spectype[-3::] +# if stype in ('enc', 'req', 'k1d') : #1D plot +# Fn = freq[-1] # Nyquist frequency +# indm = findpeaks(data, n=4) +# maxS = data.max() +## if isfield(S,'CI') && ~isempty(S.CI), +## maxS = maxS*S.CI(2); +## txtCI = [num2str(100*S.p), '% CI']; +## end +# +# Fp = freq[indm]# %peak frequency/wave number +# +# if len(indm) == 1: +# txt = [('fp = %0.2g' % Fp) + funit] +# else: +# txt = [] +# for i, fp in enumerate(Fp.tolist()): +# txt.append(('fp%d = %0.2g' % (i, fp)) + funit) +# +# txt = ''.join(txt) +# if (flag == 3): +# plotbackend.subplot(2, 1, 1) +# if (flag == 1) or (flag == 3):#% Plot in normal scale +# plotbackend.plot(np.vstack([Fp, Fp]), np.vstack([zeros(len(indm)), data.take(indm)]), ':', label=txt) +# plotbackend.plot(freq, data, linetype) +# specdata.labels.labelfig() +## if isfield(S,'CI'), +## plot(freq,S.S*S.CI(1), 'r:' ) +## plot(freq,S.S*S.CI(2), 'r:' ) +# +# a = plotbackend.axis() +# +# a1 = Fn +# if (Fp > 0): +# a1 = max(min(Fn, 10 * max(Fp)), a[1]); +# +# plotbackend.axis([0, a1 , 0, max(1.01 * maxS, a[3])]) +# #plotbackend.title('Spectral density') +# #plotbackend.xlabel(xlbl_txt) +# #plotbackend.ylabel(ylbl1_txt) +# +# +# if (flag == 3): +# plotbackend.subplot(2, 1, 2) +# +# if (flag == 2) or (flag == 3) : # Plot in logaritmic scale +# ind = np.flatnonzero(data > 0) +# +# plotbackend.plot(np.vstack([Fp, Fp]), +# np.vstack((min(10 * log10(data.take(ind) / maxS)).repeat(len(Fp)), +# 10 * log10(data.take(indm) / maxS))), ':',label=txt) +## hold on +## if isfield(S,'CI'), +## plot(freq(ind),10*log10(S.S(ind)*S.CI(1)/maxS), 'r:' ) +## plot(freq(ind),10*log10(S.S(ind)*S.CI(2)/maxS), 'r:' ) +## end +# plotbackend.plot(freq[ind], 10 * log10(data[ind] / maxS), linetype) +# +# a = plotbackend.axis() +# +# a1 = Fn +# if (Fp > 0): +# a1 = max(min(Fn, 10 * max(Fp)), a[1]); +# +# plotbackend.axis([0, a1 , -20, max(1.01 * 10 * log10(1), a[3])]) +# +# specdata.labels.labelfig() +# #plotbackend.title('Spectral density') +# #plotbackend.xlabel(xlbl_txt) +# plotbackend.ylabel(ylbl2_txt) +# +# if LegendOn: +# plotbackend.legend() # if isfield(S,'CI'), # legend(txt{:},txtCI,1) # else @@ -384,7 +384,7 @@ def plotspec(specdata, linetype='b-', flag=1): # h = polar([0 2*pi],[0 freq(end)]); # delete(h);hold on # [X,Y]=meshgrid(S.theta,freq); -# [X,Y]=pol2cart(X,Y); +# [X,Y]=polar2cart(X,Y); # contour(X,Y,S.S',lintype) # else # if (abs(thmax-thmin)<3*pi), % angle given in radians @@ -458,7 +458,7 @@ def plotspec(specdata, linetype='b-', flag=1): # %h=polar([0 2*pi],[0 freq(end)]); # %delete(h);hold on # [X,Y]=meshgrid(S.theta-phi,freq); -# [X,Y]=pol2cart(X,Y); +# [X,Y]=polar2cart(X,Y); # mesh(X,Y,S.S') # % display the unit circle beneath the surface # hold on, mesh(X,Y,zeros(size(S.S'))),hold off @@ -626,7 +626,7 @@ def plotspec(specdata, linetype='b-', flag=1): -class SpecData1D(WafoData): +class SpecData1D(PlotData): """ Container class for 1D spectrum data objects in WAFO @@ -661,7 +661,7 @@ class SpecData1D(WafoData): See also -------- - WafoData + PlotData CovData """ @@ -1219,7 +1219,7 @@ class SpecData1D(WafoData): dh = h[1] - h[0] uvdens *= dh * dh - mmpdf = WafoData(uvdens, args=(h, h), xlab='max [m]', ylab='min [m]', + mmpdf = PlotData(uvdens, args=(h, h), xlab='max [m]', ylab='min [m]', title='Joint density of maximum and minimum') try: pl = [10, 30, 50, 70, 90, 95, 99, 99.9] @@ -1388,7 +1388,7 @@ class SpecData1D(WafoData): xtxt = 'period [s]' Htxt = '%s_{v =%2.5g}' % (Htxt, u) - pdf = WafoData(f / A, t * A, title=Htxt, xlab=xtxt) + pdf = PlotData(f / A, t * A, title=Htxt, xlab=xtxt) pdf.err = err / A pdf.u = u pdf.options = opts @@ -3196,7 +3196,7 @@ class SpecData1D(WafoData): self.labels.ylab = labels[1] self.labels.zlab = labels[2] -class SpecData2D(WafoData): +class SpecData2D(PlotData): """ Container class for 2D spectrum data objects in WAFO Member variables @@ -3221,7 +3221,7 @@ class SpecData2D(WafoData): See also -------- - WafoData + PlotData CovData """ @@ -3342,8 +3342,8 @@ class SpecData2D(WafoData): # Do a physical rotation of spectrum [k, k2] = meshgrid(*self.args) - [th, r] = cart2pol(k, k2) - [k, k2] = pol2cart(th + self.phi, r) + [th, r] = cart2polar(k, k2) + [k, k2] = polar2cart(th + self.phi, r) ki1, ki2 = self.args Sn = interp2(ki1, ki2, self.data, k, k2, method) self.data = np.where(np.isnan(Sn), 0, Sn) @@ -3632,7 +3632,7 @@ def main(): print('done') def test_mm_pdf(): - import numpy as np + import wafo.spectrum.models as sm Sj = sm.Jonswap(Hm0=7, Tp=11) w = np.linspace(0, 4, 256) @@ -3640,10 +3640,12 @@ def test_mm_pdf(): S = sm.SpecData1D(Sj(w), w) # Alternatively do it manually mm = S.to_mm_pdf() +def test_docstrings(): + import doctest + doctest.testmod() + + if __name__ == '__main__': - if False : #True: # - import doctest - doctest.testmod() - else: - test_mm_pdf() + test_docstrings() + #test_mm_pdf() #main() diff --git a/pywafo/src/wafo/spectrum/dispersion_relation.py b/pywafo/src/wafo/spectrum/dispersion_relation.py deleted file mode 100644 index d921a4e..0000000 --- a/pywafo/src/wafo/spectrum/dispersion_relation.py +++ /dev/null @@ -1,206 +0,0 @@ -""" -Dispersion relation module --------------------------- -k2w - Translates from wave number to frequency -w2k - Translates from frequency to wave number -""" -import warnings -#import numpy as np -from numpy import (atleast_1d, sqrt, ones_like, zeros_like, arctan2, where, tanh, any, #@UnresolvedImport - sin, cos, sign, inf, flatnonzero, finfo, cosh, abs) #@UnresolvedImport - -__all__ = ['k2w', 'w2k'] - -def k2w(k1, k2=0e0, h=inf, g=9.81, u1=0e0, u2=0e0): - ''' Translates from wave number to frequency - using the dispersion relation - - Parameters - ---------- - k1 : array-like - wave numbers [rad/m]. - k2 : array-like, optional - second dimension wave number - h : real scalar, optional - water depth [m]. - g : real scalar, optional - acceleration of gravity, see gravity - u1, u2 : real scalars, optional - current velocity [m/s] along dimension 1 and 2. - note: when u1!=0 | u2!=0 then theta is not calculated correctly - - Returns - ------- - w : ndarray - angular frequency [rad/s]. - theta : ndarray - direction [rad]. - - Dispersion relation - ------------------- - w = sqrt(g*K*tanh(K*h)) ( 0 < w < inf) - theta = arctan2(k2,k1) (-pi < theta < pi) - where - K = sqrt(k1**2+k2**2) - - The shape of w and theta is the common shape of k1 and k2 according to the - numpy broadcasting rules. - - See also - -------- - w2k - - Example - ------- - >>> from numpy import arange - >>> import wafo.spectrum.dispersion_relation as wsd - >>> wsd.k2w(arange(0.01,.5,0.2))[0] - array([ 0.3132092 , 1.43530485, 2.00551739]) - >>> wsd.k2w(arange(0.01,.5,0.2),h=20)[0] - array([ 0.13914927, 1.43498213, 2.00551724]) - ''' - - k1i, k2i, hi, gi, u1i, u2i = atleast_1d(k1, k2, h, g, u1, u2) - - if k1i.size == 0: - return zeros_like(k1i) - ku1 = k1i*u1i - ku2 = k2i*u2i - - theta = arctan2(k2, k1) - - k = sqrt(k1i**2+k2i**2) - w = where(k>0, ku1+ku2+sqrt(gi*k*tanh(k*hi)), 0.0) - - cond = (w<0) - if any(cond): - txt0 = ''' - Waves and current are in opposite directions - making some of the frequencies negative. - Here we are forcing the negative frequencies to zero. - ''' - warnings.warn(txt0) - w = where(cond, 0.0, w) # force w to zero - - return w, theta - -def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): - ''' - Translates from frequency to wave number - using the dispersion relation - - Parameters - ---------- - w : array-like - angular frequency [rad/s]. - theta : array-like, optional - direction [rad]. - h : real scalar, optional - water depth [m]. - g : real scalar or array-like of size 2. - constant of gravity [m/s**2] or 3D normalizing constant - - Returns - ------- - k1, k2 : ndarray - wave numbers [rad/m] along dimension 1 and 2. - - Description - ----------- - Uses Newton Raphson method to find the wave number k in the dispersion relation - w**2= g*k*tanh(k*h). - The solution k(w) => k1 = k(w)*cos(theta) - k2 = k(w)*sin(theta) - The size of k1,k2 is the common shape of w and theta according to numpy - broadcasting rules. If w or theta is scalar it functions as a constant - matrix of the same shape as the other. - - Example - ------- - >>> import pylab as plb - >>> import wafo.spectrum.dispersion_relation as wsd - >>> w = plb.linspace(0,3); - >>> h = plb.plot(w,w2k(w)[0]) - >>> wsd.w2k(range(4))[0] - array([ 0. , 0.1019368 , 0.4077472 , 0.91743119]) - >>> wsd.w2k(range(4),h=20)[0] - array([ 0. , 0.10503601, 0.40774726, 0.91743119]) - - >>> plb.close('all') - - See also - -------- - k2w - ''' - wi, th, hi, gi = atleast_1d(w, theta, h, g) - - if wi.size == 0: - return zeros_like(wi) - - k = 1.0*sign(wi)*wi**2.0 / gi[0] # deep water - if (hi > 10. ** 25).all(): - k2 = k*sin(th)*gi[0]/gi[-1] #size np x nf - k1 = k*cos(th) - return k1, k2 - - - if gi.size > 1: - txt0 = ''' - Finite depth in combination with 3D normalization (len(g)=2) is not implemented yet. - ''' - raise ValueError(txt0) - - - find = flatnonzero - eps = finfo(float).eps - - oshape = k.shape - wi, k, hi = wi.ravel(), k.ravel(), hi.ravel() - - # Newton's Method - # Permit no more than count_limit iterations. - hi = hi * ones_like(k) - hn = zeros_like(k) - ix = find((wi<0) | (00 and count < count_limit): - ki = k[ix] - kh = ki * hi[ix] - hn[ix] = (ki*tanh(kh)-wi[ix]**2.0/gi)/(tanh(kh)+kh/(cosh(kh)**2.0)) - knew = ki - hn[ix] - # Make sure that the current guess is not zero. - # When Newton's Method suggests steps that lead to zero guesses - # take a step 9/10ths of the way to zero: - ksmall = find(abs(knew)==0) - if ksmall.size>0: - knew[ksmall] = ki[ksmall] / 10.0 - hn[ix[ksmall]] = ki[ksmall]-knew[ksmall] - - k[ix] = knew - # disp(['Iteration ',num2str(count),' Number of points left: ' num2str(length(ix)) ]), - - ix = find((abs(hn) > sqrt(eps)*abs(k)) * abs(hn) > sqrt(eps)) - count += 1 - - if count == count_limit: - txt1 = ''' W2K did not converge. - The maximum error in the last step was: %13.8f''' % max(hn[ix]) - warnings.warn(txt1) - - k.shape = oshape - - k2 = k*sin(th) - k1 = k*cos(th) - return k1, k2 - -def main(): - import doctest - doctest.testmod() - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/pywafo/src/wafo/spectrum/test/test_dispersion_relation.py b/pywafo/src/wafo/spectrum/test/test_dispersion_relation.py index 8dac652..7a380bf 100644 --- a/pywafo/src/wafo/spectrum/test/test_dispersion_relation.py +++ b/pywafo/src/wafo/spectrum/test/test_dispersion_relation.py @@ -24,9 +24,9 @@ def test_w2k(): array([ 0. , 0.10503601, 0.40774726, 0.91743119]) ''' -def main(): +def test_doctstrings(): import doctest doctest.testmod() if __name__ == '__main__': - main() \ No newline at end of file + test_doctstrings() \ No newline at end of file diff --git a/pywafo/src/wafo/stats/__init__.py b/pywafo/src/wafo/stats/__init__.py index b263fe3..57fccda 100644 --- a/pywafo/src/wafo/stats/__init__.py +++ b/pywafo/src/wafo/stats/__init__.py @@ -6,7 +6,7 @@ Statistics package in WAFO Toolbox. """ from scipy.stats import * from core import * -import distributions +import distributions #@Reimport from wafo.stats.distributions import * import estimation \ No newline at end of file diff --git a/pywafo/src/wafo/stats/core.py b/pywafo/src/wafo/stats/core.py index 53b9a9e..1a0fedc 100644 --- a/pywafo/src/wafo/stats/core.py +++ b/pywafo/src/wafo/stats/core.py @@ -1,6 +1,6 @@ from __future__ import division import warnings -from wafo.wafodata import WafoData +from wafo.wafodata import PlotData from wafo.misc import findextrema from scipy import special import numpy as np @@ -80,7 +80,7 @@ def edf(x, method=2): else: Fz1 = arange(1, N + 1) / (N + 1) - F = WafoData(Fz1, z, xlab='x', ylab='F(x)') + F = PlotData(Fz1, z, xlab='x', ylab='F(x)') F.setplotter('step') return F @@ -147,7 +147,7 @@ def reslife(data, u=None, umin=None, umax=None, nu=None, nmin=3, alpha=0.05, plo Returns ------- - mrl : WafoData object + mrl : PlotData object Mean residual life values, i.e., mean excesses over thresholds, u. Notes @@ -217,9 +217,9 @@ def reslife(data, u=None, umin=None, umax=None, nu=None, nmin=3, alpha=0.05, plo #options.CI = [mrll,mrlu]; #options.numdata = num; titleTxt = 'Mean residual life with %d%s CI' % (100 * p, '%') - res = WafoData(mrl, u, xlab='Threshold', ylab='Mean Excess', title=titleTxt) + res = PlotData(mrl, u, xlab='Threshold', ylab='Mean Excess', title=titleTxt) res.workspace = dict(numdata=num, umin=umin, umax=umax, nu=nu, nmin=nmin, alpha=alpha) - res.children = [WafoData(vstack([mrll, mrlu]).T, u, xlab='Threshold', title=titleTxt)] + res.children = [PlotData(vstack([mrll, mrlu]).T, u, xlab='Threshold', title=titleTxt)] res.plot_args_children = [':r'] if plotflag: res.plot() @@ -249,7 +249,7 @@ def dispersion_idx(data, t=None, u=None, umin=None, umax=None, nu=None, nmin=10, Returns ------- - DI : WafoData object + DI : PlotData object Dispersion index b_u : real scalar threshold where the number of exceedances in a fixed period (Tb) is @@ -381,10 +381,10 @@ def dispersion_idx(data, t=None, u=None, umin=None, umax=None, nu=None, nmin=10, CItxt = '%d%s CI' % (100 * p, '%') titleTxt = 'Dispersion Index plot'; - res = WafoData(di, u, title=titleTxt, labx='Threshold', laby='Dispersion Index') + res = PlotData(di, u, title=titleTxt, labx='Threshold', laby='Dispersion Index') #'caption',CItxt); res.workspace = dict(umin=umin, umax=umax, nu=nu, nmin=nmin, alpha=alpha) - res.children = [WafoData(vstack([diLo * ones(nu), diUp * ones(nu)]).T, u, xlab='Threshold', title=CItxt)] + res.children = [PlotData(vstack([diLo * ones(nu), diUp * ones(nu)]).T, u, xlab='Threshold', title=CItxt)] res.plot_args_children = ['--r'] if plotflag: res.plot(di) @@ -848,9 +848,9 @@ class RegLogit(object): z1 = (y * ones((1, yrange))) == ((y * 0 + 1) * np.arange(ymin + 1, ymax+1)) z = z[:, np.flatnonzero(z.any(axis=0))]; z1 = z1[:, np.flatnonzero(z1.any(axis=0))] - [mz, nz] = z.shape - [mx, nx] = X.shape - [my, ny] = y.shape + [_mz, nz] = z.shape + [_mx, nx] = X.shape + [my, _ny] = y.shape g = (z.sum(axis=0).cumsum() / my).reshape(-1,1) theta00 = np.log(g / (1 - g)).ravel() @@ -1120,7 +1120,7 @@ class RegLogit(object): .size : size if binomial family (default 1). ''' - [mx, nx] = self.X.shape + [_mx, nx] = self.X.shape if Xnew is None: Xnew = self.X; else: @@ -1229,7 +1229,7 @@ def _test_dispersion_idx(): xn = wafo.data.sea() t, data = xn.T Ie = findpot(data,t,0,5); - di, u, ok_u = dispersion_idx(data[Ie],t[Ie],tb=100) + di, _u, _ok_u = dispersion_idx(data[Ie],t[Ie],tb=100) di.plot() # a threshold around 1 seems appropriate. di.show() pass @@ -1240,7 +1240,7 @@ def _test_findpot(): from wafo.misc import findtc x = wafo.data.sea() t, data = x[:, :].T - itc, iv = findtc(data, 0, 'dw') + itc, _iv = findtc(data, 0, 'dw') ytc, ttc = data[itc], t[itc] ymin = 2 * data.std() tmin = 10 # sec @@ -1248,7 +1248,7 @@ def _test_findpot(): yp, tp = data[I], t[I] Ie = findpot(yp, tp, ymin, tmin) ye, te = yp[Ie], tp[Ie] - h = pylab.plot(t, data, ttc,ytc,'ro', t, zeros(len(t)), ':', te, ye, 'kx', tp, yp, '+') + pylab.plot(t, data, ttc,ytc,'ro', t, zeros(len(t)), ':', te, ye, 'kx', tp, yp, '+') pylab.show() # pass @@ -1266,7 +1266,7 @@ def test_reglogit(): #b.display() #% members and methods b.summary() - [mu,plo,pup] = b.predict(fulloutput=True); + [mu,plo,pup] = b.predict(fulloutput=True) #@UnusedVariable pass #plot(x,mu,'g',x,plo,'r:',x,pup,'r:') def test_reglogit2(): @@ -1284,12 +1284,12 @@ def test_reglogit2(): def test_sklearn0(): from sklearn.linear_model import LogisticRegression - from sklearn import datasets + from sklearn import datasets #@UnusedImport # FIXME: the iris dataset has only 4 features! - iris = datasets.load_iris() - X = iris.data - y = iris.target +# iris = datasets.load_iris() +# X = iris.data +# y = iris.target X = np.sort(5*np.random.rand(40, 1)-2.5, axis=0) y = (2*(np.cos(X)>2*np.random.rand(40, 1)-1)-1).ravel() @@ -1302,7 +1302,7 @@ def test_sklearn0(): clf_LR.fit(X, y) score.append(clf_LR.score(X,y)) - plot(cvals, score) + #plot(cvals, score) def test_sklearn(): X = np.sort(5*np.random.rand(40, 1)-2.5, axis=0) @@ -1335,7 +1335,7 @@ def test_sklearn1(): y = (2*(np.cos(X)>2*np.random.rand(40, 1)-1)-1).ravel() from sklearn.svm import SVR - cvals= np.logspace(-1,4,10) +# cvals= np.logspace(-1,4,10) svr_rbf = SVR(kernel='rbf', C=1e4, gamma=0.1, probability=True) svr_lin = SVR(kernel='linear', C=1e4, probability=True) svr_poly = SVR(kernel='poly', C=1e4, degree=2, probability=True) diff --git a/pywafo/src/wafo/stats/distributions.py b/pywafo/src/wafo/stats/distributions.py index c3032d8..2bc06f7 100644 --- a/pywafo/src/wafo/stats/distributions.py +++ b/pywafo/src/wafo/stats/distributions.py @@ -971,7 +971,7 @@ class rv_generic(object): return args, loc, scale def _fix_loc(self, args, loc): - args, loc, scale = self._fix_loc_scale(args, loc) + args, loc, _scale = self._fix_loc_scale(args, loc) return args, loc # These are actually called, and should not be overwritten if you @@ -1523,7 +1523,7 @@ class rv_continuous(rv_generic): for item in ['callparams', 'default', 'before_notes']: tempdict[item] = tempdict[item].replace(\ "\n%(shapes)s : array_like\n shape parameters", "") - for i in range(2): + for i in range(2): #@UnusedVariable if self.shapes is None: # necessary because we use %(shapes)s in two forms (w w/o ", ") self.__doc__ = self.__doc__.replace("%(shapes)s, ", "") @@ -2752,7 +2752,7 @@ def _norm_logpdf(x): def _norm_cdf(x): return special.ndtr(x) def _norm_logcdf(x): - return special.log_ndtr(x) + return log(special.ndtr(x)) def _norm_ppf(q): return special.ndtri(q) class norm_gen(rv_continuous): @@ -3166,8 +3166,8 @@ class cauchy_gen(rv_continuous): return inf, inf, nan, nan def _entropy(self): return log(4*pi) - def _fitstart(data, args=None): - return (0, 1) + def _fitstart(self, data, args=None): + return (0, 1) cauchy = cauchy_gen(name='cauchy') @@ -4129,8 +4129,8 @@ class genextreme_gen(rv_continuous): """ def _argcheck(self, c): - min = np.minimum - max = np.maximum + min = np.minimum #@ReservedAssignment + max = np.maximum #@ReservedAssignment sml = floatinfo.machar.xmin #self.b = where(c > 0, 1.0 / c,inf) #self.a = where(c < 0, 1.0 / c, -inf) @@ -6440,7 +6440,7 @@ def _drv2_ppfsingle(self, q, *args): # Use basic bisection algorithm else: return c -def reverse_dict(dict): +def reverse_dict(dict): #@ReservedAssignment newdict = {} sorted_keys = copy(dict.keys()) sorted_keys.sort() @@ -6738,7 +6738,7 @@ class rv_discrete(rv_generic): for item in ['callparams', 'default', 'before_notes']: tempdict[item] = tempdict[item].replace(\ "\n%(shapes)s : array_like\n shape parameters", "") - for i in range(2): + for i in range(2): #@UnusedVariable if self.shapes is None: # necessary because we use %(shapes)s in two forms (w w/o ", ") self.__doc__ = self.__doc__.replace("%(shapes)s, ", "") @@ -7286,9 +7286,9 @@ class rv_discrete(rv_generic): if (n > 0) and (n < 5): signature = inspect.getargspec(self._stats.im_func) if (signature[2] is not None) or ('moments' in signature[0]): - dict = {'moments':{1:'m',2:'v',3:'vs',4:'vk'}[n]} + dict = {'moments':{1:'m',2:'v',3:'vs',4:'vk'}[n]} #@ReservedAssignment else: - dict = {} + dict = {} #@ReservedAssignment mu, mu2, g1, g2 = self._stats(*args,**dict) val = _moment_from_stats(n, mu, mu2, g1, g2, self._munp, args) @@ -7422,7 +7422,7 @@ class rv_discrete(rv_generic): else: invfac = 1.0 - tot = 0.0 +# tot = 0.0 low, upp = self._ppf(0.001, *args), self._ppf(0.999, *args) low = max(min(-suppnmin, low), lb) upp = min(max(suppnmin, upp), ub) @@ -7564,8 +7564,8 @@ class bernoulli_gen(binom_gen): return binom._stats(1, pr) def _entropy(self, pr): return -pr*log(pr)-(1-pr)*log1p(-pr) -) +bernoulli = bernoulli_gen(name='bernoulli',shapes="pr") # Negative binomial class nbinom_gen(rv_discrete): """A negative binomial discrete random variable. @@ -7984,22 +7984,22 @@ class randint_gen(rv_discrete): %(example)s """ - def _argcheck(self, min, max): + def _argcheck(self, min, max): #@ReservedAssignment self.a = min self.b = max-1 return (max > min) - def _pmf(self, k, min, max): - fact = 1.0 / (max - min) + def _pmf(self, k, min, max): #@ReservedAssignment + fact = 1.0 / (max - min) return fact - def _cdf(self, x, min, max): + def _cdf(self, x, min, max): #@ReservedAssignment k = floor(x) return (k-min+1)*1.0/(max-min) - def _ppf(self, q, min, max): + def _ppf(self, q, min, max): #@ReservedAssignment vals = ceil(q*(max-min)+min)-1 vals1 = (vals-1).clip(min, max) temp = self._cdf(vals1, min, max) return where(temp >= q, vals1, vals) - def _stats(self, min, max): + def _stats(self, min, max): #@ReservedAssignment m2, m1 = arr(max), arr(min) mu = (m2 + m1 - 1.0) / 2 d = m2 - m1 @@ -8007,14 +8007,14 @@ class randint_gen(rv_discrete): g1 = 0.0 g2 = -6.0/5.0*(d*d+1.0)/(d-1.0)*(d+1.0) return mu, var, g1, g2 - def _rvs(self, min, max=None): + def _rvs(self, min, max=None): #@ReservedAssignment """An array of *size* random integers >= min and < max. If max is None, then range is >=0 and < min """ return mtrand.randint(min, max, self._size) - def _entropy(self, min, max): + def _entropy(self, min, max): #@ReservedAssignment return log(max-min) randint = randint_gen(name='randint',longname='A discrete uniform '\ '(random integer)', shapes="min, max") @@ -8212,11 +8212,12 @@ def test_truncrayleigh(): phat = truncrayleigh.fit2(R, method='ml'); -def main(): +def test_doctstrings(): import matplotlib matplotlib.interactive(True) R = norm.rvs(size=100) phat = norm.fit(R) + phat = genpareto.fit(R[R > 0.7], f0=0.1, floc=0.7) #nbinom(10, 0.75).rvs(3) @@ -8273,7 +8274,7 @@ def test_genpareto(): print(phat.par) if __name__ == '__main__': - #main() + #test_doctstrings() test_genpareto() #test_truncrayleigh() #test_lognorm() diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py index 24fb821..c3adc4b 100644 --- a/pywafo/src/wafo/stats/estimation.py +++ b/pywafo/src/wafo/stats/estimation.py @@ -20,8 +20,8 @@ from scipy import optimize import numpy import numpy as np from numpy import alltrue, arange, ravel, sum, zeros, log, sqrt, exp -from numpy import (atleast_1d, any, asarray, nan, pi, reshape, repeat, - product, ndarray, isfinite) +from numpy import (atleast_1d, any, asarray, nan, pi, #reshape, #repeat, product, ndarray, + isfinite) from numpy import flatnonzero as nonzero @@ -32,7 +32,7 @@ __all__ = [ floatinfo = np.finfo(float) #arr = atleast_1d arr = asarray -all = alltrue +all = alltrue #@ReservedAssignment def chi2isf(p, df): return special.chdtri(df, p) @@ -845,11 +845,9 @@ class FitDistribution(rv_frozen): if not self.par_fix == None: numfix = len(self.i_fixed) if numfix > 0: - format = '%d,'*numfix - format = format[:-1] - format1 = '%g,'*numfix - format1 = format1[:-1] - phatistr = format % tuple(self.i_fixed) + format0 = ', '.join(['%d']*numfix) + format1 = ', '.join(['%g']*numfix) + phatistr = format0 % tuple(self.i_fixed) phatvstr = format1 % tuple(self.par[self.i_fixed]) fixstr = 'Fixed: phat[%s] = %s ' % (phatistr, phatvstr) @@ -917,7 +915,7 @@ class FitDistribution(rv_frozen): mn = (np.floor(mn / d) - 1) * d - odd * d / 2 mx = (np.ceil(mx / d) + 1) * d + odd * d / 2 limits = np.arange(mn, mx, d) - bin, limits = numpy.histogram(self.data, bins=limits, normed=True) #, new=True) + bin, limits = np.histogram(self.data, bins=limits, normed=True) #, new=True) @ReservedAssignment limits.shape = (-1, 1) xx = limits.repeat(3, axis=1) xx.shape = (-1,) @@ -986,7 +984,7 @@ class FitDistribution(rv_frozen): dx = numpy.diff(x, axis=0) tie = (dx == 0) if any(tie): - print('P-value is on the conservative side (i.e. too large) due to ties in the data!') + warnings.warn('P-value is on the conservative side (i.e. too large) due to ties in the data!') T = self.dist.nlogps(theta, x) @@ -1009,16 +1007,17 @@ class FitDistribution(rv_frozen): -def main(): +def test_doctstrings(): import doctest doctest.testmod() + def test1(): import wafo.stats as ws R = ws.weibull_min.rvs(1,size=100); phat = FitDistribution(ws.weibull_min, R, method='mps') # # Better CI for phat.par[i=0] - Lp1 = Profile(phat, i=1) + Lp1 = Profile(phat, i=1) #@UnusedVariable # Lp2 = Profile(phat, i=2) # SF = 1./990 # x = phat.isf(SF) @@ -1085,5 +1084,5 @@ def test1(): # lp = pht.profile() if __name__ == '__main__': - test1() - main() + #test1() + test_doctstrings() diff --git a/pywafo/src/wafo/stats/misc.py b/pywafo/src/wafo/stats/misc.py index 01f6d0b..0cb9108 100644 --- a/pywafo/src/wafo/stats/misc.py +++ b/pywafo/src/wafo/stats/misc.py @@ -1,6 +1,6 @@ from numpy import asarray, ndarray, ones, nan #, reshape, repeat, product -def valarray(shape,value=nan,typecode=None): +def valarray(shape, value=nan, typecode=None): """Return an array of all value. """ #out = reshape(repeat([value],product(shape,axis=0),axis=0),shape) diff --git a/pywafo/src/wafo/test/test_gaussian.py b/pywafo/src/wafo/test/test_gaussian.py index 5bed473..d83b084 100644 --- a/pywafo/src/wafo/test/test_gaussian.py +++ b/pywafo/src/wafo/test/test_gaussian.py @@ -3,9 +3,9 @@ Created on 17. juli 2010 @author: pab ''' -import numpy as np -from numpy import pi, inf -from wafo.gaussian import Rind, prbnormtndpc, prbnormndpc, prbnormnd, cdfnorm2d, prbnorm2d +import numpy as np #@UnusedImport +from numpy import pi, inf #@UnusedImport +from wafo.gaussian import Rind, prbnormtndpc, prbnormndpc, prbnormnd, cdfnorm2d, prbnorm2d #@UnusedImport def test_rind(): ''' diff --git a/pywafo/src/wafo/test/test_kdetools.py b/pywafo/src/wafo/test/test_kdetools.py index 0c5e480..11b829e 100644 --- a/pywafo/src/wafo/test/test_kdetools.py +++ b/pywafo/src/wafo/test/test_kdetools.py @@ -4,9 +4,9 @@ Created on 20. nov. 2010 @author: pab ''' -import numpy as np -from numpy import array -import wafo.kdetools as wk +import numpy as np #@UnusedImport +from numpy import array #@UnusedImport +import wafo.kdetools as wk #@UnusedImport #import pylab as plb def test0_KDE1D(): @@ -215,7 +215,7 @@ def test_smooth_params(): array([ 0.1732289 , 0.33159097, 0.3107633 ]) - >>> gauss.hisj(data) + >>> gauss.hisj(data) array([ 0.24222479, 0.74277133, 0.15492661]) >>> data = np.array([0.753557920000000, 0.727791940000000, 0.941491690000000, @@ -223,6 +223,8 @@ def test_smooth_params(): ... 0.602882730000000, 1.368836350000000, 1.747543260000000, 1.095475610000000, ... 1.016711330000000, 0.732111430000000, 0.618917190000000, 0.759034870000000, ... 1.891946900000000, 0.724338080000000, 1.929730940000000, 0.447498380000000, 1.365084520000000]) + + ''' def test_gridcount_1D(): diff --git a/pywafo/src/wafo/test/test_misc.py b/pywafo/src/wafo/test/test_misc.py index ab7a055..85db211 100644 --- a/pywafo/src/wafo/test/test_misc.py +++ b/pywafo/src/wafo/test/test_misc.py @@ -1,12 +1,12 @@ -import numpy as np -from numpy import cos, exp, linspace, pi, sin, diff, arange, ones -from numpy.random import randn -from wafo.data import sea -from wafo.misc import (JITImport, Bunch, detrendma, DotDict, findcross, ecross, findextrema, - findrfc, rfcfilter, findtp, findtc, findoutliers, - common_shape, argsreduce, stirlerr, getshipchar, betaloge, - gravity, nextpow2, discretize, discretize2, pol2cart, - cart2pol, meshgrid, tranproc) +import numpy as np #@UnusedImport +from numpy import cos, exp, linspace, pi, sin, diff, arange, ones #@UnusedImport +from numpy.random import randn #@UnusedImport +from wafo.data import sea #@UnusedImport +from wafo.misc import (JITImport, Bunch, detrendma, DotDict, findcross, ecross, findextrema, #@UnusedImport + findrfc, rfcfilter, findtp, findtc, findoutliers, #@UnusedImport + common_shape, argsreduce, stirlerr, getshipchar, betaloge, #@UnusedImport + gravity, nextpow2, discretize, polar2cart, #@UnusedImport + cart2polar, meshgrid, tranproc)#@UnusedImport def test_JITImport(): ''' @@ -231,7 +231,7 @@ def test_rfcfilter(): >>> tp03 = rfcfilter(tp,0.3) >>> tp03 array([-0.00743352, 1.08753972, -1.07206545, 1.09550837, -1.07940458, - 1.07849396, -1.0995006 , 1.08094452]) + 1.07849396, -1.0995006 , 1.08094452, 0.11983423]) ''' def test_findtp(): ''' @@ -379,9 +379,9 @@ def test_discretize(): -7.07106781e-01, -8.31469612e-01, -9.23879533e-01, -9.80785280e-01, -1.00000000e+00]) ''' -def test_discretize2(): +def test_discretize_adaptive(): ''' - >>> x, y = discretize2(np.cos,0,np.pi) + >>> x, y = discretize(np.cos,0,np.pi, method='adaptive') >>> x; y array([ 0. , 0.19634954, 0.39269908, 0.58904862, 0.78539816, 0.9817477 , 1.17809725, 1.37444679, 1.57079633, 1.76714587, @@ -398,7 +398,7 @@ def test_pol2cart_n_cart2pol(): ''' >>> r = 5 >>> t = linspace(0,pi,20) - >>> x, y = pol2cart(t,r) + >>> x, y = polar2cart(t,r) >>> x; y array([ 5. , 4.93180652, 4.72908621, 4.39736876, 3.94570255, 3.38640786, 2.73474079, 2.00847712, 1.22742744, 0.41289673, @@ -411,7 +411,7 @@ def test_pol2cart_n_cart2pol(): 4.57886663e+00, 4.18583239e+00, 3.67861955e+00, 3.07106356e+00, 2.37973697e+00, 1.62349735e+00, 8.22972951e-01, 6.12323400e-16]) - >>> ti, ri = cart2pol(x,y) + >>> ti, ri = cart2polar(x,y) >>> ti;ri array([ 0. , 0.16534698, 0.33069396, 0.49604095, 0.66138793, 0.82673491, 0.99208189, 1.15742887, 1.32277585, 1.48812284, diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py index 47d5cfc..14e9bba 100644 --- a/pywafo/src/wafo/transform/core.py +++ b/pywafo/src/wafo/transform/core.py @@ -4,7 +4,7 @@ from __future__ import division #import numpy as np from numpy import trapz, sqrt, linspace #@UnresolvedImport -from wafo.wafodata import WafoData +from wafo.wafodata import PlotData from wafo.misc import tranproc #, trangood __all__ = ['TrData', 'TrCommon'] @@ -127,8 +127,9 @@ class TrCommon(object): def _dat2gauss(self, x, *xi): pass -class TrData(WafoData, TrCommon): - __doc__ = TrCommon.__doc__.split('mean')[0].replace('','Data') + """ +class TrData(PlotData, TrCommon): + __doc__ = TrCommon.__doc__.split('mean')[0].replace('','Data' #@ReservedAssignment + ) + """ data : array-like Gaussian values, Y args : array-like @@ -186,7 +187,7 @@ class TrData(WafoData, TrCommon): ym = self.ymean-self.ysigma self.sigma = (self.gauss2dat(yp)-self.gauss2dat(ym))/2. - self.children = [WafoData((self.args-self.mean)/self.sigma, self.args)] + self.children = [PlotData((self.args-self.mean)/self.sigma, self.args)] def trdata(self): return self diff --git a/pywafo/src/wafo/transform/models.py b/pywafo/src/wafo/transform/models.py index 4134645..1f79c25 100644 --- a/pywafo/src/wafo/transform/models.py +++ b/pywafo/src/wafo/transform/models.py @@ -43,7 +43,7 @@ _example = ''' >>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world ''' class TrCommon2(TrCommon): - __doc__ = TrCommon.__doc__ + __doc__ = TrCommon.__doc__ #@ReservedAssignment def trdata(self, x=None, xnmin= -5, xnmax=5, n=513): """ Return a discretized transformation model. @@ -75,7 +75,8 @@ class TrCommon2(TrCommon): return TrData(yn, x, mean=self.mean, sigma=self.sigma) class TrHermite(TrCommon2): - __doc__ = TrCommon2.__doc__.replace('', 'Hermite') + """ + __doc__ = TrCommon2.__doc__.replace('', 'Hermite' #@ReservedAssignment + ) + """ pardef : scalar, integer 1 Winterstein et. al. (1994) parametrization [1]_ (default) 2 Winterstein (1988) parametrization [2]_ @@ -328,7 +329,8 @@ class TrHermite(TrCommon2): class TrLinear(TrCommon2): - __doc__ = TrCommon2.__doc__.replace('', 'Linear') + """ + __doc__ = TrCommon2.__doc__.replace('', 'Linear' #@ReservedAssignment + ) + """ Description ----------- The linear transformation model is monotonic linear polynomial, calibrated @@ -369,7 +371,8 @@ class TrLinear(TrCommon2): class TrOchi(TrCommon2): - __doc__ = TrCommon2.__doc__.replace('', 'Ochi') + """ + __doc__ = TrCommon2.__doc__.replace('', 'Ochi' #@ReservedAssignment + ) + """ Description ----------- diff --git a/pywafo/src/wafo/wafodata.py b/pywafo/src/wafo/wafodata.py index ff20c62..bdccbcc 100644 --- a/pywafo/src/wafo/wafodata.py +++ b/pywafo/src/wafo/wafodata.py @@ -7,7 +7,7 @@ from scipy.integrate.quadrature import cumtrapz #@UnresolvedImport from scipy.interpolate import griddata from scipy import integrate -__all__ = ['WafoData', 'AxisLabels'] +__all__ = ['PlotData', 'AxisLabels'] def empty_copy(obj): class Empty(obj.__class__): @@ -29,7 +29,7 @@ def now(): ''' return strftime("%a, %d %b %Y %H:%M:%S", gmtime()) -class WafoData(object): +class PlotData(object): ''' Container class for data objects in WAFO @@ -38,7 +38,7 @@ class WafoData(object): data : array_like args : vector for 1D, list of vectors for 2D, 3D, ... labels : AxisLabels - children : list of WafoData objects + children : list of PlotData objects Member methods -------------- @@ -52,12 +52,12 @@ class WafoData(object): >>> x = np.arange(-2, 2, 0.2) # Plot 2 objects in one call - >>> d2 = WafoData(np.sin(x), x, xlab='x', ylab='sin', title='sinus') + >>> d2 = PlotData(np.sin(x), x, xlab='x', ylab='sin', title='sinus') >>> h = d2.plot() Plot with confidence interval - >>> d3 = WafoData(np.sin(x),x) - >>> d3.children = [WafoData(np.vstack([np.sin(x)*0.9, np.sin(x)*1.2]).T,x)] + >>> d3 = PlotData(np.sin(x),x) + >>> d3.children = [PlotData(np.vstack([np.sin(x)*0.9, np.sin(x)*1.2]).T,x)] >>> d3.plot_args_children=[':r'] >>> h = d3.plot() @@ -110,9 +110,9 @@ class WafoData(object): def eval_points(self, *args, **kwds): ''' >>> x = np.linspace(0,5,20) - >>> d = WafoData(np.sin(x),x) + >>> d = PlotData(np.sin(x),x) >>> xi = np.linspace(0,5,60) - >>> di = WafoData(d.eval_points(xi, method='cubic'),xi) + >>> di = PlotData(d.eval_points(xi, method='cubic'),xi) >>> h = d.plot('.') >>> hi = di.plot() @@ -134,7 +134,7 @@ class WafoData(object): def integrate(self, a, b, **kwds): ''' >>> x = np.linspace(0,5,60) - >>> d = WafoData(np.sin(x), x) + >>> d = PlotData(np.sin(x), x) >>> d.integrate(0,np.pi/2) 0.99940054759302188 @@ -490,15 +490,15 @@ def plot2d(axis, wdata, plotflag, *args, **kwds): def test_eval_points(): plotbackend.ioff() x = np.linspace(0,5,21) - d = WafoData(np.sin(x),x) + d = PlotData(np.sin(x),x) xi = np.linspace(0,5,61) - di = WafoData(d.eval_points(xi,method='cubic'),xi) + di = PlotData(d.eval_points(xi,method='cubic'),xi) d.plot('.') di.plot() di.show() def test_integrate(): x = np.linspace(0,5,60) - d = WafoData(np.sin(x), x) + d = PlotData(np.sin(x), x) print(d.integrate(0,np.pi/2,method='simps')) def test_docstrings(): import doctest diff --git a/pywafo/src/wafo/wave_theory/core.py b/pywafo/src/wafo/wave_theory/core.py index b2e2171..d26a5e5 100644 --- a/pywafo/src/wafo/wave_theory/core.py +++ b/pywafo/src/wafo/wave_theory/core.py @@ -5,7 +5,9 @@ Created on 3. juni 2011 ''' import numpy as np from numpy import exp, expm1, inf, nan, pi, hstack, where, atleast_1d, cos, sin -from dispersion_relation import w2k, k2w +from dispersion_relation import w2k, k2w #@UnusedImport + +__all__ =['w2k', 'k2w', 'sensor_typeid', 'sensor_type', 'TransferFunction'] def hyperbolic_ratio(a, b, sa, sb): ''' @@ -474,90 +476,90 @@ class TransferFunction(object): zk = self._get_zk(kw) return hyperbolic_ratio(zk, hk, -1, -1), ee # sinh(zk)./sinh(hk), ee -def wave_pressure(z, Hm0, h=10000, g=9.81, rho=1028): - ''' - Calculate pressure amplitude due to water waves. - - Parameters - ---------- - z : array-like - depth where pressure is calculated [m] - Hm0 : array-like - significant wave height (same as the average of the 1/3'rd highest - waves in a seastate. [m] - h : real scalar - waterdepth (default 10000 [m]) - g : real scalar - acceleration of gravity (default 9.81 m/s**2) - rho : real scalar - water density (default 1028 kg/m**3) - - - Returns - ------- - p : ndarray - pressure amplitude due to water waves at water depth z. [Pa] - - PRESSURE calculate pressure amplitude due to water waves according to - linear theory. - - Example - ----- - >>> import pylab as plt - >>> z = -np.linspace(10,20) - >>> fh = plt.plot(z, wave_pressure(z, Hm0=1, h=20)) - >>> plt.show() - - See also - -------- - w2k - - - u = psweep.Fn*sqrt(mgf.length*9.81) - z = -10; h = inf; - Hm0 = 1.5;Tp = 4*sqrt(Hm0); - S = jonswap([],[Hm0,Tp]); - Hw = tran(S.w,0,[0 0 -z],'P',h) - Sm = S; - Sm.S = Hw.'.*S.S; - x1 = spec2sdat(Sm,1000); - pwave = pressure(z,Hm0,h) - - plot(psweep.x{1}/u, psweep.f) - hold on - plot(x1(1:100,1)-30,x1(1:100,2),'r') - ''' - - - # Assume seastate with jonswap spectrum: - - Tp = 4 * np.sqrt(Hm0) - gam = jonswap_peakfact(Hm0, Tp) - Tm02 = Tp / (1.30301 - 0.01698 * gam + 0.12102 / gam) - w = 2 * np.pi / Tm02 - kw, unused_kw2 = w2k(w, 0, h) - - hk = kw * h - zk1 = kw * z - zk = hk + zk1 # z measured positive upward from mean water level (default) - #zk = hk-zk1; % z measured positive downward from mean water level - #zk1 = -zk1; - #zk = zk1; % z measured positive upward from sea floor - - # cosh(zk)/cosh(hk) approx exp(zk) for large h - # hyperbolic_ratio(zk,hk,1,1) = cosh(zk)/cosh(hk) - # pr = np.where(np.pi < hk, np.exp(zk1), hyperbolic_ratio(zk, hk, 1, 1)) - pr = hyperbolic_ratio(zk, hk, 1, 1) - pressure = (rho * g * Hm0 / 2) * pr - -# pos = [np.zeros_like(z),np.zeros_like(z),z] -# tf = TransferFunction(pos=pos, sensortype='p', h=h, rho=rho, g=g) -# Hw, Gwt = tf.tran(w,0) -# pressure2 = np.abs(Hw) * Hm0 / 2 - - return pressure +#def wave_pressure(z, Hm0, h=10000, g=9.81, rho=1028): +# ''' +# Calculate pressure amplitude due to water waves. +# +# Parameters +# ---------- +# z : array-like +# depth where pressure is calculated [m] +# Hm0 : array-like +# significant wave height (same as the average of the 1/3'rd highest +# waves in a seastate. [m] +# h : real scalar +# waterdepth (default 10000 [m]) +# g : real scalar +# acceleration of gravity (default 9.81 m/s**2) +# rho : real scalar +# water density (default 1028 kg/m**3) +# +# +# Returns +# ------- +# p : ndarray +# pressure amplitude due to water waves at water depth z. [Pa] +# +# PRESSURE calculate pressure amplitude due to water waves according to +# linear theory. +# +# Example +# ----- +# >>> import pylab as plt +# >>> z = -np.linspace(10,20) +# >>> fh = plt.plot(z, wave_pressure(z, Hm0=1, h=20)) +# >>> plt.show() +# +# See also +# -------- +# w2k +# +# +# u = psweep.Fn*sqrt(mgf.length*9.81) +# z = -10; h = inf; +# Hm0 = 1.5;Tp = 4*sqrt(Hm0); +# S = jonswap([],[Hm0,Tp]); +# Hw = tran(S.w,0,[0 0 -z],'P',h) +# Sm = S; +# Sm.S = Hw.'.*S.S; +# x1 = spec2sdat(Sm,1000); +# pwave = pressure(z,Hm0,h) +# +# plot(psweep.x{1}/u, psweep.f) +# hold on +# plot(x1(1:100,1)-30,x1(1:100,2),'r') +# ''' +# +# +# # Assume seastate with jonswap spectrum: +# +# Tp = 4 * np.sqrt(Hm0) +# gam = jonswap_peakfact(Hm0, Tp) +# Tm02 = Tp / (1.30301 - 0.01698 * gam + 0.12102 / gam) +# w = 2 * np.pi / Tm02 +# kw, unused_kw2 = w2k(w, 0, h) +# +# hk = kw * h +# zk1 = kw * z +# zk = hk + zk1 # z measured positive upward from mean water level (default) +# #zk = hk-zk1; % z measured positive downward from mean water level +# #zk1 = -zk1; +# #zk = zk1; % z measured positive upward from sea floor +# +# # cosh(zk)/cosh(hk) approx exp(zk) for large h +# # hyperbolic_ratio(zk,hk,1,1) = cosh(zk)/cosh(hk) +# # pr = np.where(np.pi < hk, np.exp(zk1), hyperbolic_ratio(zk, hk, 1, 1)) +# pr = hyperbolic_ratio(zk, hk, 1, 1) +# pressure = (rho * g * Hm0 / 2) * pr +# +## pos = [np.zeros_like(z),np.zeros_like(z),z] +## tf = TransferFunction(pos=pos, sensortype='p', h=h, rho=rho, g=g) +## Hw, Gwt = tf.tran(w,0) +## pressure2 = np.abs(Hw) * Hm0 / 2 +# +# return pressure def main(): - sensortype(range(21)) + sensor_type(range(21)) if __name__ == '__main__': pass \ No newline at end of file diff --git a/pywafo/src/wafo/wave_theory/dispersion_relation.py b/pywafo/src/wafo/wave_theory/dispersion_relation.py index 45fed47..d921a4e 100644 --- a/pywafo/src/wafo/wave_theory/dispersion_relation.py +++ b/pywafo/src/wafo/wave_theory/dispersion_relation.py @@ -6,7 +6,7 @@ w2k - Translates from frequency to wave number """ import warnings #import numpy as np -from numpy import (atleast_1d, sqrt, zeros_like, arctan2, where, tanh, any, #@UnresolvedImport +from numpy import (atleast_1d, sqrt, ones_like, zeros_like, arctan2, where, tanh, any, #@UnresolvedImport sin, cos, sign, inf, flatnonzero, finfo, cosh, abs) #@UnresolvedImport __all__ = ['k2w', 'w2k'] @@ -132,15 +132,15 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): -------- k2w ''' - wi, th, gi = atleast_1d(w, theta, g) + wi, th, hi, gi = atleast_1d(w, theta, h, g) if wi.size == 0: return zeros_like(wi) - k = 1.0*sign(wi)*wi**2.0 #% deep water - if h > 10. ** 25: - k2 = k*sin(th)/gi[-1] #%size np x nf - k1 = k*cos(th)/gi[0] + k = 1.0*sign(wi)*wi**2.0 / gi[0] # deep water + if (hi > 10. ** 25).all(): + k2 = k*sin(th)*gi[0]/gi[-1] #size np x nf + k1 = k*cos(th) return k1, k2 @@ -155,11 +155,11 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): eps = finfo(float).eps oshape = k.shape - wi, k = wi.ravel(), k.ravel() + wi, k, hi = wi.ravel(), k.ravel(), hi.ravel() # Newton's Method # Permit no more than count_limit iterations. - + hi = hi * ones_like(k) hn = zeros_like(k) ix = find((wi<0) | (00 and count < count_limit): ki = k[ix] - hn[ix] = (ki*tanh(ki*h)-wi[ix]**2.0/gi)/(tanh(ki*h)+ki*h/(cosh(ki*h)**2.0)) + kh = ki * hi[ix] + hn[ix] = (ki*tanh(kh)-wi[ix]**2.0/gi)/(tanh(kh)+kh/(cosh(kh)**2.0)) knew = ki - hn[ix] # Make sure that the current guess is not zero. # When Newton's Method suggests steps that lead to zero guesses @@ -181,7 +182,7 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): hn[ix[ksmall]] = ki[ksmall]-knew[ksmall] k[ix] = knew - # disp(['Iteration ',num2str(count),' Number of points left: ' num2str(length(ix)) ]), + # disp(['Iteration ',num2str(count),' Number of points left: ' num2str(length(ix)) ]), ix = find((abs(hn) > sqrt(eps)*abs(k)) * abs(hn) > sqrt(eps)) count += 1