diff --git a/wafo/containers.py b/wafo/containers.py index 70a7133..1771053 100644 --- a/wafo/containers.py +++ b/wafo/containers.py @@ -160,6 +160,12 @@ class PlotData(object): return interpolate.griddata( self.args, self.data, points, **options) + def to_cdf(self): + if isinstance(self.args, (list, tuple)): # Multidimensional data + raise NotImplementedError('integration for ndim>1 not implemented') + cdf = np.hstack((0, integrate.cumtrapz(self.data, self.args))) + return PlotData(cdf, np.copy(self.args), xlab='x', ylab='F(x)') + def integrate(self, a, b, **kwds): ''' >>> x = np.linspace(0,5,60) @@ -551,7 +557,7 @@ def plot2d(axis, wdata, plotflag, *args, **kwds): clvec = np.sort(CL) if plotflag in [1, 8, 9]: - h = axis.contour(*args1, levels=CL, **kwds) + h = axis.contour(*args1, levels=clvec, **kwds) # else: # [cs hcs] = contour3(f.x{:},f.f,CL,sym); @@ -563,9 +569,7 @@ def plot2d(axis, wdata, plotflag, *args, **kwds): 'Only the first 12 levels will be listed in table.') clvals = PL[:ncl] if isPL else clvec[:ncl] - unused_axcl = cltext( - clvals, - percent=isPL) # print contour level text + unused_axcl = cltext(clvals, percent=isPL) elif any(plotflag == [7, 9]): axis.clabel(h) else: diff --git a/wafo/doc/tutorial_scripts/chapter1.py b/wafo/doc/tutorial_scripts/chapter1.py index 96a84ae..a67dd73 100644 --- a/wafo/doc/tutorial_scripts/chapter1.py +++ b/wafo/doc/tutorial_scripts/chapter1.py @@ -1,6 +1,5 @@ -from scipy import * -from pylab import * - +import wafo.plotbackend.plotbackend as plt +import numpy as np # pyreport -o chapter1.html chapter1.py #! CHAPTER1 demonstrates some applications of WAFO @@ -15,16 +14,16 @@ from pylab import * #! Section 1.4 Some applications of WAFO #!--------------------------------------- -#! Section 1.4.1 Simulation from spectrum, estimation of spectrum +#! Section 1.4.1 Simulation from spectrum, estimation of spectrum #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #! Simulation of the sea surface from spectrum #! The following code generates 200 seconds of data sampled with 10Hz from #! the Torsethaugen spectrum import wafo.spectrum.models as wsm -S = wsm.Torsethaugen(Hm0=6, Tp=8); +S = wsm.Torsethaugen(Hm0=6, Tp=8) S1 = S.tospecdata() S1.plot() -show() +plt.show() ## @@ -32,23 +31,23 @@ import wafo.objects as wo xs = S1.sim(ns=2000, dt=0.1) ts = wo.mat2timeseries(xs) ts.plot_wave('-') -show() +plt.show() -#! Estimation of spectrum +#! Estimation of spectrum #!~~~~~~~~~~~~~~~~~~~~~~~ #! A common situation is that one wants to estimate the spectrum for wave #! measurements. The following code simulate 20 minutes signal sampled at 4Hz #! and compare the spectral estimate with the original Torsethaugen spectum. -clf() -Fs = 4; -xs = S1.sim(ns=fix(20 * 60 * Fs), dt=1. / Fs) -ts = wo.mat2timeseries(xs) +plt.clf() +Fs = 4 +xs = S1.sim(ns=np.fix(20 * 60 * Fs), dt=1. / Fs) +ts = wo.mat2timeseries(xs) Sest = ts.tospecdata(L=400) S1.plot() Sest.plot('--') -axis([0, 3, 0, 5]) # This may depend on the simulation -show() +plt.axis([0, 3, 0, 5]) +plt.show() #! Section 1.4.2 Probability distributions of wave characteristics. #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -58,7 +57,7 @@ show() #! In the following example we study the trough period extracted from the #! time series and compared with the theoretical density computed with exact #! spectrum, S1, and the estimated spectrum, Sest. -clf() +plt.clf() import wafo.misc as wm dtyex = S1.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3) dtyest = Sest.to_t_pdf(pdef='Tt', paramt=(0, 10, 51), nit=3) @@ -69,29 +68,29 @@ wm.plot_histgrm(T, bins=bins, normed=True) dtyex.plot() dtyest.plot('-.') -axis([0, 10, 0, 0.35]) -show() +plt.axis([0, 10, 0, 0.35]) +plt.show() #! Section 1.4.3 Directional spectra #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -#! Here are a few lines of code, which produce directional spectra +#! Here are a few lines of code, which produce directional spectra #! with frequency independent and frequency dependent spreading. -clf() +plt.clf() plotflag = 1 -Nt = 101; # number of angles -th0 = pi / 2; # primary direction of waves -Sp = 15; # spreading parameter +Nt = 101 # number of angles +th0 = np.pi / 2 # primary direction of waves +Sp = 15 # spreading parameter -D1 = wsm.Spreading(type='cos', theta0=th0, method=None) # frequency independent -D12 = wsm.Spreading(type='cos', theta0=0, method='mitsuyasu') # frequency dependent +D1 = wsm.Spreading(type='cos', theta0=th0, method=None) +D12 = wsm.Spreading(type='cos', theta0=0, method='mitsuyasu') SD1 = D1.tospecdata2d(S1) SD12 = D12.tospecdata2d(S1) SD1.plot() -SD12.plot()#linestyle='dashdot') -show() +SD12.plot() # linestyle='dashdot') +plt.show() -#! 3D Simulation of the sea surface +#! 3D Simulation of the sea surface #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #! The simulations show that frequency dependent spreading leads to #! much more irregular surface so the orientation of waves is less @@ -120,7 +119,7 @@ show() #! The figure is not shown in the Tutorial # # Nx = 3; Ny = 2; Nt = 2 ^ 12; dx = 10; dy = 10;dt = 0.5; -# F = seasim(SD12, Nx, Ny, Nt, dx, dy, dt, 1, 0); +# F = seasim(SD12, Nx, Ny, Nt, dx, dy, dt, 1, 0); # Z = permute(F.Z, [3 1 2]); # [X, Y] = meshgrid(F.x, F.y); # N = Nx * Ny; @@ -137,16 +136,16 @@ show() #! Section 1.4.4 Fatigue, Load cycles and Markov models. #! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -#! Switching Markow chain of turningpoints +#! Switching Markow chain of turningpoints #! In fatigue applications the exact sample path is not important, but #! only the tops and bottoms of the load, called the sequence of turning #! points (TP). From the turning points one can extract load cycles, from #! which damage calculations and fatigue life predictions can be #! performed. #! -#! The commands below computes the intensity of rainflowcycles for -#! the Gaussian model with spectrum S1 using the Markov approximation. -#! The rainflow cycles found in the simulated load signal are shown in the +#! The commands below computes the intensity of rainflowcycles for +#! the Gaussian model with spectrum S1 using the Markov approximation. +#! The rainflow cycles found in the simulated load signal are shown in the #! figure. #clf() @@ -164,30 +163,30 @@ show() #! Section 1.4.5 Extreme value statistics #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Plot of yura87 data -clf() +plt.clf() import wafo.data as wd xn = wd.yura87() -#xn = load('yura87.dat'); -subplot(211) -plot(xn[::30, 0] / 3600, xn[::30, 1], '.') -title('Water level') -ylabel('(m)') +#xn = load('yura87.dat'); +plt.subplot(211) +plt.plot(xn[::30, 0] / 3600, xn[::30, 1], '.') +plt.title('Water level') +plt.ylabel('(m)') #! Formation of 5 min maxima yura = xn[:85500, 1] yura = np.reshape(yura, (285, 300)).T maxyura = yura.max(axis=0) -subplot(212) -plot(xn[299:85500:300, 0] / 3600, maxyura, '.') -xlabel('Time (h)') -ylabel('(m)') -title('Maximum 5 min water level') -show() +plt.subplot(212) +plt.plot(xn[299:85500:300, 0] / 3600, maxyura, '.') +plt.xlabel('Time (h)') +plt.ylabel('(m)') +plt.title('Maximum 5 min water level') +plt.show() #! Estimation of GEV for yuramax -clf() +plt.clf() import wafo.stats as ws phat = ws.genextreme.fit2(maxyura, method='ml') phat.plotfitsummary() -show() +plt.show() #disp('Block = 11, Last block') diff --git a/wafo/doc/tutorial_scripts/chapter2.py b/wafo/doc/tutorial_scripts/chapter2.py index fe4616a..8210878 100644 --- a/wafo/doc/tutorial_scripts/chapter2.py +++ b/wafo/doc/tutorial_scripts/chapter2.py @@ -1,6 +1,5 @@ +import wafo.plotbackend.plotbackend as plt import numpy as np -from scipy import * -from pylab import * # pyreport -o chapter2.html chapter2.py @@ -10,15 +9,15 @@ from pylab import * #! Chapter2 contains the commands used in Chapter 2 of the tutorial and #! present some tools for analysis of random functions with #! respect to their correlation, spectral and distributional properties. -#! The presentation is divided into three examples: +#! The presentation is divided into three examples: #! #! Example1 is devoted to estimation of different parameters in the model. #! Example2 deals with spectral densities and #! Example3 presents the use of WAFO to simulate samples of a Gaussian #! process. #! -#! Some of the commands are edited for fast computation. -#! +#! Some of the commands are edited for fast computation. +#! #! Section 2.1 Introduction and preliminary analysis #!==================================================== #! Example 1: Sea data @@ -37,15 +36,15 @@ tp = ts.turning_points() cc = tp.cycle_pairs() lc = cc.level_crossings() lc.plot() -show() +plt.show() #! Average number of upcrossings per time unit #!---------------------------------------------- -#! Next we compute the mean frequency as the average number of upcrossings -#! per time unit of the mean level (= 0); this may require interpolation in the -#! crossing intensity curve, as follows. +#! Next we compute the mean frequency as the average number of upcrossings +#! per time unit of the mean level (= 0); this may require interpolation in the +#! crossing intensity curve, as follows. T = xx[:, 0].max() - xx[:, 0].min() -f0 = np.interp(0, lc.args, lc.data, 0) / T #! zero up-crossing frequency +f0 = np.interp(0, lc.args, lc.data, 0) / T # zero up-crossing frequency print('f0 = %g' % f0) #! Turningpoints and irregularity factor @@ -60,17 +59,17 @@ print('fm = %g, alpha = %g, ' % (fm, alfa)) #!------------------------ #! We finish this section with some remarks about the quality #! of the measured data. Especially sea surface measurements can be -#! of poor quality. We shall now check the quality of the dataset {\tt xx}. -#! It is always good practice to visually examine the data -#! before the analysis to get an impression of the quality, +#! of poor quality. We shall now check the quality of the dataset {\tt xx}. +#! It is always good practice to visually examine the data +#! before the analysis to get an impression of the quality, #! non-linearities and narrow-bandedness of the data. -#! First we shall plot the data and zoom in on a specific region. +#! First we shall plot the data and zoom in on a specific region. #! A part of sea data is visualized with the following commands -clf() +plt.clf() ts.plot_wave('k-', tp, '*', nfig=1, nsub=1) -axis([0, 2, -2, 2]) -show() +plt.axis([0, 2, -2, 2]) +plt.show() #! Finding possible spurious points #!------------------------------------ @@ -91,58 +90,59 @@ inds, indg = wm.findoutliers(ts.data, zcrit, dcrit, ddcrit, verbose=True) #!---------------------------------------------------- #! Periodogram: Raw spectrum #! -clf() +plt.clf() Lmax = 9500 S = ts.tospecdata(L=Lmax) S.plot() -axis([0, 5, 0, 0.7]) -show() +plt.axis([0, 5, 0, 0.7]) +plt.show() -#! Calculate moments +#! Calculate moments #!------------------- mom, text = S.moment(nr=4) -print('sigma = %g, m0 = %g' % (sa, sqrt(mom[0]))) +print('sigma = %g, m0 = %g' % (sa, np.sqrt(mom[0]))) #! Section 2.2.1 Random functions in Spectral Domain - Gaussian processes #!-------------------------------------------------------------------------- -#! Smoothing of spectral estimate +#! Smoothing of spectral estimate #!---------------------------------- #! By decreasing Lmax the spectrum estimate becomes smoother. -clf() -Lmax0 = 200; Lmax1 = 50 +plt.clf() +Lmax0 = 200 +Lmax1 = 50 S1 = ts.tospecdata(L=Lmax0) S2 = ts.tospecdata(L=Lmax1) S1.plot('-.') S2.plot() -show() +plt.show() #! Estimated autocovariance #!---------------------------- #! Obviously knowing the spectrum one can compute the covariance -#! function. The following code will compute the covariance for the -#! unimodal spectral density S1 and compare it with estimated +#! function. The following code will compute the covariance for the +#! unimodal spectral density S1 and compare it with estimated #! covariance of the signal xx. -clf() +plt.clf() Lmax = 85 -R1 = S1.tocovdata(nr=1) +R1 = S1.tocovdata(nr=1) Rest = ts.tocovdata(lag=Lmax) R1.plot('.') Rest.plot() -axis([0, 25, -0.1, 0.25]) -show() +plt.axis([0, 25, -0.1, 0.25]) +plt.show() -#! We can see in Figure below that the covariance function corresponding to -#! the spectral density S2 significantly differs from the one estimated -#! directly from data. -#! It can be seen in Figure above that the covariance corresponding to S1 +#! We can see in Figure below that the covariance function corresponding to +#! the spectral density S2 significantly differs from the one estimated +#! directly from data. +#! It can be seen in Figure above that the covariance corresponding to S1 #! agrees much better with the estimated covariance function -clf() +plt.clf() R2 = S2.tocovdata(nr=1) R2.plot('.') Rest.plot() -show() +plt.show() #! Section 2.2.2 Transformed Gaussian models #!------------------------------------------- @@ -154,68 +154,69 @@ rho3 = ws.skew(xx[:, 1]) rho4 = ws.kurtosis(xx[:, 1]) sk, ku = S1.stats_nl(moments='sk') - + #! Comparisons of 3 transformations -clf() +plt.clf() import wafo.transform.models as wtm gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku).trdata() -g = wtm.TrLinear(mean=me, sigma=sa).trdata() # Linear transformation +g = wtm.TrLinear(mean=me, sigma=sa).trdata() # Linear transformation glc, gemp = lc.trdata(mean=me, sigma=sa) -glc.plot('b-') #! Transf. estimated from level-crossings -gh.plot('b-.') #! Hermite Transf. estimated from moments +glc.plot('b-') # Transf. estimated from level-crossings +gh.plot('b-.') # Hermite Transf. estimated from moments g.plot('r') -grid('on') -show() - +plt.grid('on') +plt.show() + #! Test Gaussianity of a stochastic process #!------------------------------------------ -#! TESTGAUSSIAN simulates e(g(u)-u) = int (g(u)-u)^2 du for Gaussian processes +#! TESTGAUSSIAN simulates e(g(u)-u) = int (g(u)-u)^2 du for Gaussian processes #! given the spectral density, S. The result is plotted if test0 is given. #! This is useful for testing if the process X(t) is Gaussian. #! If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level. -#! -#! As we see from the figure below: none of the simulated values of test1 is -#! above 1.00. Thus the data significantly departs from a Gaussian distribution. -clf() +#! +#! As we see from the figure below: none of the simulated values of test1 is +#! above 1.00. Thus the data significantly departs from a Gaussian distribution. +plt.clf() test0 = glc.dist2gauss() #! the following test takes time N = len(xx) test1 = S1.testgaussian(ns=N, cases=50, test0=test0) -is_gaussian = sum(test1 > test0) > 5 +is_gaussian = sum(test1 > test0) > 5 print(is_gaussian) -show() +plt.show() #! Normalplot of data xx #!------------------------ #! indicates that the underlying distribution has a "heavy" upper tail and a -#! "light" lower tail. -clf() +#! "light" lower tail. +plt.clf() import pylab ws.probplot(ts.data.ravel(), dist='norm', plot=pylab) -show() +plt.show() #! Section 2.2.3 Spectral densities of sea data #!----------------------------------------------- #! Example 2: Different forms of spectra #! import wafo.spectrum.models as wsm -clf() -Hm0 = 7; Tp = 11; +plt.clf() +Hm0 = 7 +Tp = 11 spec = wsm.Jonswap(Hm0=Hm0, Tp=Tp).tospecdata() spec.plot() -show() +plt.show() #! Directional spectrum and Encountered directional spectrum #! Directional spectrum -clf() +plt.clf() D = wsm.Spreading('cos2s') Sd = D.tospecdata2d(spec) Sd.plot() -show() +plt.show() ##!Encountered directional spectrum -##!--------------------------------- +##!--------------------------------- #clf() #Se = spec2spec(Sd,'encdir',0,10); #plotspec(Se), hold on @@ -254,10 +255,10 @@ show() #disp('Block = 20'),pause(pstate) # ##!#! Section 2.3 Simulation of transformed Gaussian process -##!#! Example 3: Simulation of random sea +##!#! Example 3: Simulation of random sea ##! The reconstruct function replaces the spurious points of seasurface by ##! simulated data on the basis of the remaining data and a transformed Gaussian -##! process. As noted previously one must be careful using the criteria +##! process. As noted previously one must be careful using the criteria ##! for finding spurious points when reconstructing a dataset, because ##! these criteria might remove the highest and steepest waves as we can see ##! in this plot where the spurious points is indicated with a '+' sign: @@ -269,13 +270,13 @@ show() ##!wafostamp('','(ER)') #disp('Block = 21'),pause(pstate) # -##! Compare transformation (grec) from reconstructed (y) +##! Compare transformation (grec) from reconstructed (y) ##! with original (glc) from (xx) #clf #trplot(g), hold on #plot(gemp(:,1),gemp(:,2)) #plot(glc(:,1),glc(:,2),'-.') -#plot(grec(:,1),grec(:,2)), hold off +#plot(grec(:,1),grec(:,2)), hold off #disp('Block = 22'),pause(pstate) # ##!#! @@ -284,7 +285,7 @@ show() #x = dat2gaus(y,grec); #Sx = dat2spec(x,L); #disp('Block = 23'),pause(pstate) -# +# ##!#! #clf #dt = spec2dt(Sx) @@ -294,27 +295,28 @@ show() #waveplot(ysim,'-') ##!wafostamp('','(CR)') #disp('Block = 24'),pause(pstate) -# +# #! Estimated spectrum compared to Torsethaugen spectrum #!------------------------------------------------------- - -clf() -fp = 1.1;dw = 0.01 + +plt.clf() +fp = 1.1 +dw = 0.01 H0 = S1.characteristic('Hm0')[0] -St = wsm.Torsethaugen(Hm0=H0,Tp=2*pi/fp).tospecdata(np.arange(0,5+dw/2,dw)) +St = wsm.Torsethaugen(Hm0=H0,Tp=2*np.pi/fp).tospecdata(np.arange(0,5+dw/2,dw)) S1.plot() St.plot('-.') -axis([0, 6, 0, 0.4]) -show() - +plt.axis([0, 6, 0, 0.4]) +plt.show() + #! Transformed Gaussian model compared to Gaussian model #!-------------------------------------------------------- dt = St.sampling_period() -va, sk, ku = St.stats_nl(moments='vsk' ) +va, sk, ku = St.stats_nl(moments='vsk') #sa = sqrt(va) gh = wtm.TrHermite(mean=me, sigma=sa, skew=sk, kurt=ku, ysigma=sa) - + ysim_t = St.sim(ns=240, dt=0.5) xsim_t = ysim_t.copy() xsim_t[:,1] = gh.gauss2dat(ysim_t[:,1]) @@ -322,4 +324,4 @@ xsim_t[:,1] = gh.gauss2dat(ysim_t[:,1]) ts_y = wo.mat2timeseries(ysim_t) ts_x = wo.mat2timeseries(xsim_t) ts_y.plot_wave(sym1='r.', ts=ts_x, sym2='b', sigma=sa, nsub=5, nfig=1) -show() +plt.show() diff --git a/wafo/doc/tutorial_scripts/chapter3.py b/wafo/doc/tutorial_scripts/chapter3.py index ec94b60..852d88a 100644 --- a/wafo/doc/tutorial_scripts/chapter3.py +++ b/wafo/doc/tutorial_scripts/chapter3.py @@ -1,55 +1,57 @@ +from wafo.plotbackend import plotbackend as plt import numpy as np -from scipy import * -from pylab import * #! CHAPTER3 Demonstrates distributions of wave characteristics #!============================================================= #! #! Chapter3 contains the commands used in Chapter3 in the tutorial. -#! -#! Some of the commands are edited for fast computation. +#! +#! Some of the commands are edited for fast computation. #! #! Section 3.2 Estimation of wave characteristics from data #!---------------------------------------------------------- #! Example 1 -#!~~~~~~~~~~ - +#!~~~~~~~~~~ + speed = 'fast' #speed = 'slow' - +import scipy.signal as ss import wafo.data as wd import wafo.misc as wm import wafo.objects as wo -xx = wd.sea() -xx[:,1] = wm.detrend(xx[:,1]) +import wafo.stats as ws +import wafo.spectrum.models as wsm +xx = wd.sea() +xx[:, 1] = ss.detrend(xx[:, 1]) ts = wo.mat2timeseries(xx) Tcrcr, ix = ts.wave_periods(vh=0, pdef='c2c', wdef='tw', rate=8) Tc, ixc = ts.wave_periods(vh=0, pdef='u2d', wdef='tw', rate=8) - + #! Histogram of crestperiod compared to the kernel density estimate #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ import wafo.kdetools as wk -clf() +plt.clf() print(Tc.mean()) print(Tc.max()) -t = linspace(0.01,8,200); +t = np.linspace(0.01,8,200); ftc = wk.TKDE(Tc, L2=0, inc=128) -plot(t,ftc.eval_grid(t), t, ftc.eval_grid_fast(t),'-.') -wm.plot_histgrm(Tc,normed=True) -title('Kernel Density Estimates') -axis([0, 8, 0, 0.5]) -show() - +plt.plot(t,ftc.eval_grid(t), t, ftc.eval_grid_fast(t),'-.') +wm.plot_histgrm(Tc, normed=True) +plt.title('Kernel Density Estimates') +plt.xlabel('Tc [s]') +plt.axis([0, 8, 0, 0.5]) +plt.show() + #! Extreme waves - model check: the highest and steepest wave #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -clf() +plt.clf() S, H = ts.wave_height_steepness(method=0) indS = S.argmax() indH = H.argmax() ts.plot_sp_wave([indH, indS],'k.') -show() +plt.show() #! Does the highest wave contradict a transformed Gaussian model? #!---------------------------------------------------------------- @@ -68,9 +70,9 @@ show() #muLstd = tranproc(mu1o-lamb*mu1oStd,fliplr(grec1)); #muUstd = tranproc(mu1o+lamb*mu1oStd,fliplr(grec1)); #plot (y1(inds1,1), [muLstd muUstd],'b-') -#axis([1482 1498 -1 3]), +#axis([1482 1498 -1 3]), #wafostamp([],'(ER)') -#disp('Block = 6'), +#disp('Block = 6'), #pause(pstate) # ##!#! Expected value (solid) compared to data removed @@ -83,520 +85,530 @@ show() #! Crest height PDF #!------------------ #! Transform data so that kde works better -clf() +plt.clf() wave_data = ts.wave_parameters() Ac = wave_data['Ac'] -L2 = 0.6; -import pylab -ws.probplot(Ac**L2, dist='norm', plot=pylab) -show() +L2 = 0.6 + +ws.probplot(Ac**L2, dist='norm', plot=plt) +plt.show() #!#! -#! -fac = kde(Ac,{'L2',L2},linspace(0.01,3,200)); -pdfplot(fac) -wafostamp([],'(ER)') -simpson(fac.x{1},fac.f) -disp('Block = 8'), pause(pstate) +plt.clf()# +fac = wk.TKDE(Ac,L2=L2)(np.linspace(0.01,3,200), output='plot') +fac.plot() +# wafostamp([],'(ER)') +fac.integrate(a=0.01, b=3) +print('Block = 8'), +# pause(pstate) #!#! Empirical crest height CDF -clf -Fac = flipud(cumtrapz(fac.x{1},flipud(fac.f))); -Fac = [fac.x{1} 1-Fac]; -Femp = plotedf(Ac,Fac); -axis([0 2 0 1]) -wafostamp([],'(ER)') -disp('Block = 9'), pause(pstate) +plt.clf() +Fac = fac.to_cdf() +Femp = ws.edf(Ac) +Fac.plot() +Femp.plot() +plt.axis([0, 2, 0, 1]) + +#wafostamp([],'(ER)') +#disp('Block = 9'), pause(pstate) #!#! Empirical crest height CDF compared to a Transformed Rayleigh approximation -facr = trraylpdf(fac.x{1},'Ac',grec1); -Facr = cumtrapz(facr.x{1},facr.f); -hold on -plot(facr.x{1},Facr,'.') -axis([1.25 2.25 0.95 1]) -wafostamp([],'(ER)') -disp('Block = 10'), pause(pstate) + +# facr = trraylpdf(fac.x{1},'Ac',grec1); +# Facr = cumtrapz(facr.x{1},facr.f); +# hold on +# plot(facr.x{1},Facr,'.') +# axis([1.25 2.25 0.95 1]) +# wafostamp([],'(ER)') +# disp('Block = 10'), pause(pstate) #!#! Joint pdf of crest period and crest amplitude -clf -kopt2 = kdeoptset('L2',0.5,'inc',256); -Tc = Tcf+Tcb; -fTcAc = kdebin([Tc Ac],kopt2); -fTcAc.labx={'Tc [s]' 'Ac [m]'} -pdfplot(fTcAc) -hold on -plot(Tc,Ac,'k.') -hold off -wafostamp([],'(ER)') -disp('Block = 11'), pause(pstate) +plt.clf() +Tcf = wave_data['Tcf'] +Tcb = wave_data['Tcb'] +Tc = Tcf + Tcb +fTcAc = wk.TKDE([Tc, Ac],L2=0.5, inc=256).eval_grid_fast(output='plot') +fTcAc.labels.labx = 'Tc [s]' +fTcAc.labels.laby = 'Ac [m]' +fTcAc.plot() +plt.hold(True) +plt.plot(Tc, Ac,'k.') +plt.hold(False) +plt.show() +#wafostamp([],'(ER)') +#disp('Block = 11'), pause(pstate) -#!#! Example 4: Simple wave characteristics obtained from Jonswap spectrum -clf -S = jonswap([],[5 10]); -[m, mt]= spec2mom(S,4,[],0); -disp('Block = 12'), pause(pstate) +#!#! Example 4: Simple wave characteristics obtained from Jonswap spectrum +plt.clf() +S = wsm.Jonswap(Hm0=5, Tp=10).tospecdata() +m, mt = S.moment(nr=4, even=False) +print(m) +print(mt) +# disp('Block = 12'), pause(pstate) -clf -spec2bw(S) -[ch Sa2] = spec2char(S,[1 3]) -disp('Block = 13'), pause(pstate) +plt.clf() +S.bandwidth(['alpha']) +ch, Sa2 = S.characteristic(['Hm0', 'Tm02']) + +# disp('Block = 13'), pause(pstate) #!#! Section 3.3.2 Explicit form approximations of wave characteristic densities #!#! Longuett-Higgins model for Tc and Ac -clf -t = linspace(0,15,100); -h = linspace(0,6,100); -flh = lh83pdf(t,h,[m(1),m(2),m(3)]); -disp('Block = 14'), pause(pstate) - -#!#! Transformed Longuett-Higgins model for Tc and Ac -clf -[sk, ku ]=spec2skew(S); -sa = sqrt(m(1)); -gh = hermitetr([],[sa sk ku 0]); -flhg = lh83pdf(t,h,[m(1),m(2),m(3)],gh); -disp('Block = 15'), pause(pstate) +# plt.clf() +# t = np.linspace(0,15,100) +# h = np.linspace(0,6,100) +# flh = lh83pdf(t, h, [m[0],m[1], m[2]) +# #disp('Block = 14'), pause(pstate) +# +# #!#! Transformed Longuett-Higgins model for Tc and Ac +# clf +# [sk, ku ]=spec2skew(S); +# sa = sqrt(m(1)); +# gh = hermitetr([],[sa sk ku 0]); +# flhg = lh83pdf(t,h,[m(1),m(2),m(3)],gh); +# disp('Block = 15'), pause(pstate) #!#! Cavanie model for Tc and Ac -clf -t = linspace(0,10,100); -h = linspace(0,7,100); -fcav = cav76pdf(t,h,[m(1) m(2) m(3) m(5)],[]); -disp('Block = 16'), pause(pstate) - -#!#! Example 5 Transformed Rayleigh approximation of crest- vs trough- amplitude -clf -xx = load('sea.dat'); -x = xx; -x(:,2) = detrend(x(:,2)); -SS = dat2spec2(x); -[sk, ku, me, si ] = spec2skew(SS); -gh = hermitetr([],[si sk ku me]); -Hs = 4*si; -r = (0:0.05:1.1*Hs)'; -fac_h = trraylpdf(r,'Ac',gh); -fat_h = trraylpdf(r,'At',gh); -h = (0:0.05:1.7*Hs)'; -facat_h = trraylpdf(h,'AcAt',gh); -pdfplot(fac_h) -hold on -pdfplot(fat_h,'--') -hold off -wafostamp([],'(ER)') -disp('Block = 17'), pause(pstate) - -#!#! -clf -TC = dat2tc(xx, me); -tc = tp2mm(TC); -Ac = tc(:,2); -At = -tc(:,1); -AcAt = Ac+At; -disp('Block = 18'), pause(pstate) - -#!#! -clf -Fac_h = [fac_h.x{1} cumtrapz(fac_h.x{1},fac_h.f)]; -subplot(3,1,1) -Fac = plotedf(Ac,Fac_h); -hold on -plot(r,1-exp(-8*r.^2/Hs^2),'.') -axis([1. 2. 0.9 1]) -title('Ac CDF') - -Fat_h = [fat_h.x{1} cumtrapz(fat_h.x{1},fat_h.f)]; -subplot(3,1,2) -Fat = plotedf(At,Fat_h); -hold on -plot(r,1-exp(-8*r.^2/Hs^2),'.') -axis([1. 2. 0.9 1]) -title('At CDF') - -Facat_h = [facat_h.x{1} cumtrapz(facat_h.x{1},facat_h.f)]; -subplot(3,1,3) -Facat = plotedf(AcAt,Facat_h); -hold on -plot(r,1-exp(-2*r.^2/Hs^2),'.') -axis([1.5 3.5 0.9 1]) -title('At+Ac CDF') - -wafostamp([],'(ER)') -disp('Block = 19'), pause(pstate) - -#!#! Section 3.4 Exact wave distributions in transformed Gaussian Sea -#!#! Section 3.4.1 Density of crest period, crest length or encountered crest period -clf -S1 = torsethaugen([],[6 8],1); -D1 = spreading(101,'cos',pi/2,[15],[],0); -D12 = spreading(101,'cos',0,[15],S1.w,1); -SD1 = mkdspec(S1,D1); -SD12 = mkdspec(S1,D12); -disp('Block = 20'), pause(pstate) - -#!#! Crest period -clf -tic -f_tc = spec2tpdf(S1,[],'Tc',[0 11 56],[],4); -toc -pdfplot(f_tc) -wafostamp([],'(ER)') -simpson(f_tc.x{1},f_tc.f) -disp('Block = 21'), pause(pstate) - -#!#! Crest length - -if strncmpi(speed,'slow',1) - opt1 = rindoptset('speed',5,'method',3); - opt2 = rindoptset('speed',5,'nit',2,'method',0); -else - #! fast - opt1 = rindoptset('speed',7,'method',3); - opt2 = rindoptset('speed',7,'nit',2,'method',0); -end - - -clf -if strncmpi(speed,'slow',1) - NITa = 5; -else - disp('NIT=5 may take time, running with NIT=3 in the following') - NITa = 3; -end -#!f_Lc = spec2tpdf2(S1,[],'Lc',[0 200 81],opt1); #! Faster and more accurate -f_Lc = spec2tpdf(S1,[],'Lc',[0 200 81],[],NITa); -pdfplot(f_Lc,'-.') -wafostamp([],'(ER)') -disp('Block = 22'), pause(pstate) - - -f_Lc_1 = spec2tpdf(S1,[],'Lc',[0 200 81],1.5,NITa); -#!f_Lc_1 = spec2tpdf2(S1,[],'Lc',[0 200 81],1.5,opt1); - -hold on -pdfplot(f_Lc_1) -wafostamp([],'(ER)') - -disp('Block = 23'), pause(pstate) -#!#! -clf -simpson(f_Lc.x{1},f_Lc.f) -simpson(f_Lc_1.x{1},f_Lc_1.f) - -disp('Block = 24'), pause(pstate) -#!#! -clf -tic - -f_Lc_d1 = spec2tpdf(rotspec(SD1,pi/2),[],'Lc',[0 300 121],[],NITa); -f_Lc_d12 = spec2tpdf(SD12,[],'Lc',[0 200 81],[],NITa); -#! f_Lc_d1 = spec2tpdf2(rotspec(SD1,pi/2),[],'Lc',[0 300 121],opt1); -#! f_Lc_d12 = spec2tpdf2(SD12,[],'Lc',[0 200 81],opt1); -toc -pdfplot(f_Lc_d1,'-.'), hold on -pdfplot(f_Lc_d12), hold off -wafostamp([],'(ER)') - -disp('Block = 25'), pause(pstate) - -#!#! - - -clf -opt1 = rindoptset('speed',5,'method',3); -SD1r = rotspec(SD1,pi/2); -if strncmpi(speed,'slow',1) - f_Lc_d1_5 = spec2tpdf(SD1r,[], 'Lc',[0 300 121],[],5); - pdfplot(f_Lc_d1_5), hold on -else - #! fast - disp('Run the following example only if you want a check on computing time') - disp('Edit the command file and remove #!') -end -f_Lc_d1_3 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],3); -f_Lc_d1_2 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],2); -f_Lc_d1_0 = spec2tpdf(SD1r,[],'Lc',[0 300 121],[],0); -#!f_Lc_d1_n4 = spec2tpdf2(SD1r,[],'Lc',[0 400 161],opt1); - -pdfplot(f_Lc_d1_3), hold on -pdfplot(f_Lc_d1_2) -pdfplot(f_Lc_d1_0) -#!pdfplot(f_Lc_d1_n4) - -#!simpson(f_Lc_d1_n4.x{1},f_Lc_d1_n4.f) - -disp('Block = 26'), pause(pstate) - -#!#! Section 3.4.2 Density of wave period, wave length or encountered wave period -#!#! Example 7: Crest period and high crest waves -clf -tic -xx = load('sea.dat'); -x = xx; -x(:,2) = detrend(x(:,2)); -SS = dat2spec(x); -si = sqrt(spec2mom(SS,1)); -SS.tr = dat2tr(x); -Hs = 4*si -method = 0; -rate = 2; -[S, H, Ac, At, Tcf, Tcb, z_ind, yn] = dat2steep(x,rate,method); -Tc = Tcf+Tcb; -t = linspace(0.01,8,200); -ftc1 = kde(Tc,{'L2',0},t); -pdfplot(ftc1) -hold on -#! f_t = spec2tpdf(SS,[],'Tc',[0 8 81],0,4); -f_t = spec2tpdf(SS,[],'Tc',[0 8 81],0,2); -simpson(f_t.x{1},f_t.f) -pdfplot(f_t,'-.') -hold off -wafostamp([],'(ER)') -toc -disp('Block = 27'), pause(pstate) - -#!#! -clf -tic - -if strncmpi(speed,'slow',1) - NIT = 4; -else - NIT = 2; -end -#! f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],[Hs/2],4); -tic -f_t2 = spec2tpdf(SS,[],'Tc',[0 8 81],Hs/2,NIT); -toc - -Pemp = sum(Ac>Hs/2)/sum(Ac>0) -simpson(f_t2.x{1},f_t2.f) -index = find(Ac>Hs/2); -ftc1 = kde(Tc(index),{'L2',0},t); -ftc1.f = Pemp*ftc1.f; -pdfplot(ftc1) -hold on -pdfplot(f_t2,'-.') -hold off -wafostamp([],'(ER)') -toc -disp('Block = 28'), pause(pstate) - -#!#! Example 8: Wave period for high crest waves -#! clf - tic - f_tcc2 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],-1); -toc - simpson(f_tcc2.x{1},f_tcc2.f) - f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],3,5); -#! f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],1,5); - simpson(f_tcc3.x{1},f_tcc3.f) - pdfplot(f_tcc2,'-.') - hold on - pdfplot(f_tcc3) - hold off - toc -disp('Block = 29'), pause(pstate) - -#!#! -clf -[TC tc_ind v_ind] = dat2tc(yn,[],'dw'); -N = length(tc_ind); -t_ind = tc_ind(1:2:N); -c_ind = tc_ind(2:2:N); -Pemp = sum(yn(t_ind,2)<-Hs/2 & yn(c_ind,2)>Hs/2)/length(t_ind) -ind = find(yn(t_ind,2)<-Hs/2 & yn(c_ind,2)>Hs/2); -spwaveplot(yn,ind(2:4)) -wafostamp([],'(ER)') -disp('Block = 30'), pause(pstate) - -#!#! -clf -Tcc = yn(v_ind(1+2*ind),1)-yn(v_ind(1+2*(ind-1)),1); -t = linspace(0.01,14,200); -ftcc1 = kde(Tcc,{'kernel' 'epan','L2',0},t); -ftcc1.f = Pemp*ftcc1.f; -pdfplot(ftcc1,'-.') -wafostamp([],'(ER)') -disp('Block = 31'), pause(pstate) - -tic -f_tcc22_1 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[Hs/2],-1); -toc -simpson(f_tcc22_1.x{1},f_tcc22_1.f) -hold on -pdfplot(f_tcc22_1) -hold off -wafostamp([],'(ER)') -disp('Block = 32'), pause(pstate) - -disp('The rest of this chapter deals with joint densities.') -disp('Some calculations may take some time.') -disp('You could experiment with other NIT.') -#!return - -#!#! Section 3.4.3 Joint density of crest period and crest height -#!#! Example 9. Some preliminary analysis of the data -clf -tic -yy = load('gfaksr89.dat'); -SS = dat2spec(yy); -si = sqrt(spec2mom(SS,1)); -SS.tr = dat2tr(yy); -Hs = 4*si -v = gaus2dat([0 0],SS.tr); -v = v(2) -toc -disp('Block = 33'), pause(pstate) - -#!#! -clf -tic -[TC, tc_ind, v_ind] = dat2tc(yy,v,'dw'); -N = length(tc_ind); -t_ind = tc_ind(1:2:N); -c_ind = tc_ind(2:2:N); -v_ind_d = v_ind(1:2:N+1); -v_ind_u = v_ind(2:2:N+1); -T_d = ecross(yy(:,1),yy(:,2),v_ind_d,v); -T_u = ecross(yy(:,1),yy(:,2),v_ind_u,v); - -Tc = T_d(2:end)-T_u(1:end); -Tt = T_u(1:end)-T_d(1:end-1); -Tcf = yy(c_ind,1)-T_u; -Ac = yy(c_ind,2)-v; -At = v-yy(t_ind,2); -toc -disp('Block = 34'), pause(pstate) - -#!#! -clf -tic -t = linspace(0.01,15,200); -kopt3 = kdeoptset('hs',0.25,'L2',0); -ftc1 = kde(Tc,kopt3,t); -ftt1 = kde(Tt,kopt3,t); -pdfplot(ftt1,'k') -hold on -pdfplot(ftc1,'k-.') -f_tc4 = spec2tpdf(SS,[],'Tc',[0 12 81],0,4,5); -f_tc2 = spec2tpdf(SS,[],'Tc',[0 12 81],0,2,5); -f_tc = spec2tpdf(SS,[],'Tc',[0 12 81],0,-1); -pdfplot(f_tc,'b') -hold off -legend('kde(Tt)','kde(Tc)','f_{tc}') -wafostamp([],'(ER)') -toc -disp('Block = 35'), pause(pstate) - -#!#! Example 10: Joint characteristics of a half wave: -#!#! position and height of a crest for a wave with given period -clf -tic -ind = find(4.4Hs/2); -plot(Tc(ind), Ac(ind),'.'); -hold on -pdfplot(flh_g,'k-.') -pdfplot(f_tcac_s) -toc -wafostamp([],'(ER)') -disp('Block = 39'), pause(pstate) - -#!#! -clf -#! f_tcac = spec2thpdf(SS,[],'TcAc',[0 12 81],[0:0.2:8],opt1); -#! pdfplot(f_tcac) -disp('Block = 40'), pause(pstate) - -#!#! Section 3.4.4 Joint density of crest and trough height -#!#! Section 3.4.5 Min-to-max distributions � Markov method -#!#! Example 11. (min-max problems with Gullfaks data) -#!#! Joint density of maximum and the following minimum -clf -tic -tp = dat2tp(yy); -Mm = fliplr(tp2mm(tp)); -fmm = kde(Mm); -f_mM = spec2mmtpdf(SS,[],'mm',[],[-7 7 51],opt2); - -pdfplot(f_mM,'-.') -hold on -pdfplot(fmm,'k-') -hold off -wafostamp([],'(ER)') -toc -disp('Block = 41'), pause(pstate) - -#!#! The joint density of �still water separated� maxima and minima. -clf -tic -ind = find(Mm(:,1)>v & Mm(:,2)Hs/2)/sum(Ac>0) +# simpson(f_t2.x{1},f_t2.f) +# index = find(Ac>Hs/2); +# ftc1 = kde(Tc(index),{'L2',0},t); +# ftc1.f = Pemp*ftc1.f; +# pdfplot(ftc1) +# hold on +# pdfplot(f_t2,'-.') +# hold off +# wafostamp([],'(ER)') +# toc +# disp('Block = 28'), pause(pstate) +# +# #!#! Example 8: Wave period for high crest waves +# #! clf +# tic +# f_tcc2 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],-1); +# toc +# simpson(f_tcc2.x{1},f_tcc2.f) +# f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],3,5); +# #! f_tcc3 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[0],1,5); +# simpson(f_tcc3.x{1},f_tcc3.f) +# pdfplot(f_tcc2,'-.') +# hold on +# pdfplot(f_tcc3) +# hold off +# toc +# disp('Block = 29'), pause(pstate) +# +# #!#! +# clf +# [TC tc_ind v_ind] = dat2tc(yn,[],'dw'); +# N = length(tc_ind); +# t_ind = tc_ind(1:2:N); +# c_ind = tc_ind(2:2:N); +# Pemp = sum(yn(t_ind,2)<-Hs/2 & yn(c_ind,2)>Hs/2)/length(t_ind) +# ind = find(yn(t_ind,2)<-Hs/2 & yn(c_ind,2)>Hs/2); +# spwaveplot(yn,ind(2:4)) +# wafostamp([],'(ER)') +# disp('Block = 30'), pause(pstate) +# +# #!#! +# clf +# Tcc = yn(v_ind(1+2*ind),1)-yn(v_ind(1+2*(ind-1)),1); +# t = linspace(0.01,14,200); +# ftcc1 = kde(Tcc,{'kernel' 'epan','L2',0},t); +# ftcc1.f = Pemp*ftcc1.f; +# pdfplot(ftcc1,'-.') +# wafostamp([],'(ER)') +# disp('Block = 31'), pause(pstate) +# +# tic +# f_tcc22_1 = spec2tccpdf(SS,[],'t>',[0 12 61],[Hs/2],[Hs/2],-1); +# toc +# simpson(f_tcc22_1.x{1},f_tcc22_1.f) +# hold on +# pdfplot(f_tcc22_1) +# hold off +# wafostamp([],'(ER)') +# disp('Block = 32'), pause(pstate) +# +# disp('The rest of this chapter deals with joint densities.') +# disp('Some calculations may take some time.') +# disp('You could experiment with other NIT.') +# #!return +# +# #!#! Section 3.4.3 Joint density of crest period and crest height +# #!#! Example 9. Some preliminary analysis of the data +# clf +# tic +# yy = load('gfaksr89.dat'); +# SS = dat2spec(yy); +# si = sqrt(spec2mom(SS,1)); +# SS.tr = dat2tr(yy); +# Hs = 4*si +# v = gaus2dat([0 0],SS.tr); +# v = v(2) +# toc +# disp('Block = 33'), pause(pstate) +# +# #!#! +# clf +# tic +# [TC, tc_ind, v_ind] = dat2tc(yy,v,'dw'); +# N = length(tc_ind); +# t_ind = tc_ind(1:2:N); +# c_ind = tc_ind(2:2:N); +# v_ind_d = v_ind(1:2:N+1); +# v_ind_u = v_ind(2:2:N+1); +# T_d = ecross(yy(:,1),yy(:,2),v_ind_d,v); +# T_u = ecross(yy(:,1),yy(:,2),v_ind_u,v); +# +# Tc = T_d(2:end)-T_u(1:end); +# Tt = T_u(1:end)-T_d(1:end-1); +# Tcf = yy(c_ind,1)-T_u; +# Ac = yy(c_ind,2)-v; +# At = v-yy(t_ind,2); +# toc +# disp('Block = 34'), pause(pstate) +# +# #!#! +# clf +# tic +# t = linspace(0.01,15,200); +# kopt3 = kdeoptset('hs',0.25,'L2',0); +# ftc1 = kde(Tc,kopt3,t); +# ftt1 = kde(Tt,kopt3,t); +# pdfplot(ftt1,'k') +# hold on +# pdfplot(ftc1,'k-.') +# f_tc4 = spec2tpdf(SS,[],'Tc',[0 12 81],0,4,5); +# f_tc2 = spec2tpdf(SS,[],'Tc',[0 12 81],0,2,5); +# f_tc = spec2tpdf(SS,[],'Tc',[0 12 81],0,-1); +# pdfplot(f_tc,'b') +# hold off +# legend('kde(Tt)','kde(Tc)','f_{tc}') +# wafostamp([],'(ER)') +# toc +# disp('Block = 35'), pause(pstate) +# +# #!#! Example 10: Joint characteristics of a half wave: +# #!#! position and height of a crest for a wave with given period +# clf +# tic +# ind = find(4.4Hs/2); +# plot(Tc(ind), Ac(ind),'.'); +# hold on +# pdfplot(flh_g,'k-.') +# pdfplot(f_tcac_s) +# toc +# wafostamp([],'(ER)') +# disp('Block = 39'), pause(pstate) +# +# #!#! +# clf +# #! f_tcac = spec2thpdf(SS,[],'TcAc',[0 12 81],[0:0.2:8],opt1); +# #! pdfplot(f_tcac) +# disp('Block = 40'), pause(pstate) +# +# #!#! Section 3.4.4 Joint density of crest and trough height +# #!#! Section 3.4.5 Min-to-max distributions Markov method +# #!#! Example 11. (min-max problems with Gullfaks data) +# #!#! Joint density of maximum and the following minimum +# clf +# tic +# tp = dat2tp(yy); +# Mm = fliplr(tp2mm(tp)); +# fmm = kde(Mm); +# f_mM = spec2mmtpdf(SS,[],'mm',[],[-7 7 51],opt2); +# +# pdfplot(f_mM,'-.') +# hold on +# pdfplot(fmm,'k-') +# hold off +# wafostamp([],'(ER)') +# toc +# disp('Block = 41'), pause(pstate) +# +# #!#! The joint density of still water separated maxima and minima. +# clf +# tic +# ind = find(Mm(:,1)>v & Mm(:,2) easier to evaluate using cellmode evaluation. -% Created by GL July 13, 2000 -% from commands used in Chapter 5 -% - -%% Chapter 5 Extreme value analysis - -%% Section 5.1 Weibull and Gumbel papers +## CHAPTER5 contains the commands used in Chapter 5 of the tutorial +# +# CALL: Chapter5 +# +# Some of the commands are edited for fast computation. +# Each set of commands is followed by a 'pause' command. +# + +# Tested on Matlab 5.3 +# History +# Added Return values by GL August 2008 +# Revised pab sept2005 +# Added sections -> easier to evaluate using cellmode evaluation. +# Created by GL July 13, 2000 +# from commands used in Chapter 5 +# + +## Chapter 5 Extreme value analysis + +## Section 5.1 Weibull and Gumbel papers +from __future__ import division +import numpy as np +import scipy.interpolate as si +from wafo.plotbackend import plotbackend as plt +import wafo.data as wd +import wafo.objects as wo +import wafo.stats as ws +import wafo.kdetools as wk pstate = 'off' -% Significant wave-height data on Weibull paper, -clf -Hs = load('atlantic.dat'); -wei = plotweib(Hs) -wafostamp([],'(ER)') -disp('Block = 1'),pause(pstate) - -%% -% Significant wave-height data on Gumbel paper, -clf -gum=plotgumb(Hs) -wafostamp([],'(ER)') - disp('Block = 2'),pause(pstate) - -%% -% Significant wave-height data on Normal probability paper, -plotnorm(log(Hs),1,0); -wafostamp([],'(ER)') -disp('Block = 3'),pause(pstate) - -%% -% Return values in the Gumbel distribution -clf -T=1:100000; -sT=gum(2) - gum(1)*log(-log(1-1./T)); -semilogx(T,sT), hold -N=1:length(Hs); Nmax=max(N); -plot(Nmax./N,sort(Hs,'descend'),'.') -title('Return values in the Gumbel model') -xlabel('Return period') -ylabel('Return value') -wafostamp([],'(ER)') -disp('Block = 4'),pause(pstate) - -%% Section 5.2 Generalized Pareto and Extreme Value distributions -%% Section 5.2.1 Generalized Extreme Value distribution - -% Empirical distribution of significant wave-height with estimated -% Generalized Extreme Value distribution, -gev=fitgev(Hs,'plotflag',2) -wafostamp([],'(ER)') -disp('Block = 5a'),pause(pstate) - -clf -x = linspace(0,14,200); -plotkde(Hs,[x;pdfgev(x,gev)]') -disp('Block = 5b'),pause(pstate) - -% Analysis of yura87 wave data. -% Wave data interpolated (spline) and organized in 5-minute intervals -% Normalized to mean 0 and std = 1 to get stationary conditions. -% maximum level over each 5-minute interval analysed by GEV -xn = load('yura87.dat'); -XI = 0:0.25:length(xn); -N = length(XI); N = N-mod(N,4*60*5); -YI = interp1(xn(:,1),xn(:,2),XI(1:N),'spline'); -YI = reshape(YI,4*60*5,N/(4*60*5)); % Each column holds 5 minutes of - % interpolated data. -Y5 = (YI-ones(1200,1)*mean(YI))./(ones(1200,1)*std(YI)); -Y5M = max(Y5); -Y5gev = fitgev(Y5M,'method','mps','plotflag',2) -wafostamp([],'(ER)') -disp('Block = 6'),pause(pstate) - -%% Section 5.2.2 Generalized Pareto distribution - -% Exceedances of significant wave-height data over level 3, -gpd3 = fitgenpar(Hs(Hs>3)-3,'plotflag',1); -wafostamp([],'(ER)') - -%% -figure -% Exceedances of significant wave-height data over level 7, -gpd7 = fitgenpar(Hs(Hs>7),'fixpar',[nan,nan,7],'plotflag',1); -wafostamp([],'(ER)') -disp('Block = 6'),pause(pstate) - -%% -%Simulates 100 values from the GEV distribution with parameters (0.3, 1, 2), then estimates the -%parameters using two different methods and plots the estimated distribution functions together -%with the empirical distribution. -Rgev = rndgev(0.3,1,2,1,100); -gp = fitgev(Rgev,'method','pwm'); -gm = fitgev(Rgev,'method','ml','start',gp.params,'plotflag',0); - -x=sort(Rgev); -plotedf(Rgev,gp,{'-','r-'}); -hold on -plot(x,cdfgev(x,gm),'--') -hold off -wafostamp([],'(ER)') -disp('Block =7'),pause(pstate) - -%% -% ; -Rgpd = rndgenpar(0.4,1,0,1,100); -plotedf(Rgpd); -hold on -gp = fitgenpar(Rgpd,'method','pkd','plotflag',0); -x=sort(Rgpd); -plot(x,cdfgenpar(x,gp)) -% gm = fitgenpar(Rgpd,'method','mom','plotflag',0); -% plot(x,cdfgenpar(x,gm),'g--') -gw = fitgenpar(Rgpd,'method','pwm','plotflag',0); -plot(x,cdfgenpar(x,gw),'g:') -gml = fitgenpar(Rgpd,'method','ml','plotflag',0); -plot(x,cdfgenpar(x,gml),'--') -gmps = fitgenpar(Rgpd,'method','mps','plotflag',0); -plot(x,cdfgenpar(x,gmps),'r-.') -hold off -wafostamp([],'(ER)') -disp('Block = 8'),pause(pstate) - -%% -% Return values for the GEV distribution -T = logspace(1,5,10); -[sT, sTlo, sTup] = invgev(1./T,Y5gev,'lowertail',false,'proflog',true); - -%T = 2:100000; -%k=Y5gev.params(1); mu=Y5gev.params(3); sigma=Y5gev.params(2); -%sT1 = invgev(1./T,Y5gev,'lowertail',false); -%sT=mu + sigma/k*(1-(-log(1-1./T)).^k); -clf -semilogx(T,sT,T,sTlo,'r',T,sTup,'r'), hold -N=1:length(Y5M); Nmax=max(N); -plot(Nmax./N,sort(Y5M,'descend'),'.') -title('Return values in the GEV model') -xlabel('Return priod') -ylabel('Return value') -grid on -disp('Block = 9'),pause(pstate) - -%% Section 5.3 POT-analysis - -% Estimated expected exceedance over level u as function of u. -clf -plotreslife(Hs,'umin',2,'umax',10,'Nu',200); -wafostamp([],'(ER)') -disp('Block = 10'),pause(pstate) - -%% -% Estimated distribution functions of monthly maxima -%with the POT method (solid), -% fitting a GEV (dashed) and the empirical distribution. - -% POT- method -gpd7 = fitgenpar(Hs(Hs>7)-7,'method','pwm','plotflag',0); -khat = gpd7.params(1); -sigmahat = gpd7.params(2); -muhat = length(Hs(Hs>7))/(7*3*2); -bhat = sigmahat/muhat^khat; -ahat = 7-(bhat-sigmahat)/khat; -x = linspace(5,15,200); -plot(x,cdfgev(x,khat,bhat,ahat)) -disp('Block = 11'),pause(pstate) - -%% -% Since we have data to compute the monthly maxima mm over -%42 months we can also try to fit a -% GEV distribution directly: -mm = zeros(1,41); -for i=1:41 - mm(i)=max(Hs(((i-1)*14+1):i*14)); -end - -gev=fitgev(mm); -hold on -plotedf(mm) -plot(x,cdfgev(x,gev),'--') -hold off -wafostamp([],'(ER)') -disp('Block = 12, Last block'),pause(pstate) +# Significant wave-height data on Weibull paper, + +fig = plt.figure() +ax = fig.add_subplot(111) +Hs = wd.atlantic() +wei = ws.weibull_min.fit(Hs) +tmp = ws.probplot(Hs, wei, ws.weibull_min, plot=ax) +plt.show() +#wafostamp([],'(ER)') +#disp('Block = 1'),pause(pstate) + +## +# Significant wave-height data on Gumbel paper, +plt.clf() +ax = fig.add_subplot(111) +gum = ws.gumbel_r.fit(Hs) +tmp1 = ws.probplot(Hs, gum, ws.gumbel_r, plot=ax) +#wafostamp([],'(ER)') +plt.show() +#disp('Block = 2'),pause(pstate) + +## +# Significant wave-height data on Normal probability paper, +plt.clf() +ax = fig.add_subplot(111) +phat = ws.norm.fit2(np.log(Hs)) +phat.plotresq() +#tmp2 = ws.probplot(np.log(Hs), phat, ws.norm, plot=ax) + +#wafostamp([],'(ER)') +plt.show() +#disp('Block = 3'),pause(pstate) + +## +# Return values in the Gumbel distribution +plt.clf() +T = np.r_[1:100000] +sT = gum[0] - gum[1] * np.log(-np.log1p(-1./T)) +plt.semilogx(T, sT) +plt.hold(True) +# ws.edf(Hs).plot() +Nmax = len(Hs) +N = np.r_[1:Nmax+1] + +plt.plot(Nmax/N, sorted(Hs, reverse=True), '.') +plt.title('Return values in the Gumbel model') +plt.xlabel('Return period') +plt.ylabel('Return value') +#wafostamp([],'(ER)') +plt.show() +#disp('Block = 4'),pause(pstate) + +## Section 5.2 Generalized Pareto and Extreme Value distributions +## Section 5.2.1 Generalized Extreme Value distribution + +# Empirical distribution of significant wave-height with estimated +# Generalized Extreme Value distribution, +gev = ws.genextreme.fit2(Hs) +gev.plotfitsummary() +# wafostamp([],'(ER)') +# disp('Block = 5a'),pause(pstate) + +plt.clf() +x = np.linspace(0,14,200) +kde = wk.TKDE(Hs, L2=0.5)(x, output='plot') +kde.plot() +plt.hold(True) +plt.plot(x, gev.pdf(x),'--') +# disp('Block = 5b'),pause(pstate) + +# Analysis of yura87 wave data. +# Wave data interpolated (spline) and organized in 5-minute intervals +# Normalized to mean 0 and std = 1 to get stationary conditions. +# maximum level over each 5-minute interval analysed by GEV +xn = wd.yura87() +XI = np.r_[1:len(xn):0.25] - .99 +N = len(XI) +N = N - np.mod(N, 4*60*5) + +YI = si.interp1d(xn[:, 0],xn[:, 1], kind='linear')(XI) +YI = YI.reshape(4*60*5, N/(4*60*5)) # Each column holds 5 minutes of + # interpolated data. +Y5 = (YI - YI.mean(axis=0)) / YI.std(axis=0) +Y5M = Y5.maximum(axis=0) +Y5gev = ws.genextreme.fit2(Y5M,method='mps') +Y5gev.plotfitsummary() +#wafostamp([],'(ER)') +#disp('Block = 6'),pause(pstate) + +## Section 5.2.2 Generalized Pareto distribution + +# Exceedances of significant wave-height data over level 3, +gpd3 = ws.genpareto.fit2(Hs[Hs>3]-3, floc=0) +gpd3.plotfitsummary() +#wafostamp([],'(ER)') + +## +plt.figure() +# Exceedances of significant wave-height data over level 7, +gpd7 = ws.genpareto.fit2(Hs(Hs>7), floc=7) +gpd7.plotfitsummary() +# wafostamp([],'(ER)') +# disp('Block = 6'),pause(pstate) + +## +#Simulates 100 values from the GEV distribution with parameters (0.3, 1, 2), +# then estimates the parameters using two different methods and plots the +# estimated distribution functions together with the empirical distribution. +Rgev = ws.genextreme.rvs(0.3,1,2,size=100) +gp = ws.genextreme.fit2(Rgev, method='mps'); +gm = ws.genextreme.fit2(Rgev, *gp.par.tolist(), method='ml') +gm.plotfitsummary() + +gp.plotecdf() +plt.hold(True) +plt.plot(x, gm.cdf(x), '--') +plt.hold(False) +#wafostamp([],'(ER)') +#disp('Block =7'),pause(pstate) + +## +# ; +Rgpd = ws.genpareto.rvs(0.4,0, 1,size=100) +gp = ws.genpareto.fit2(Rgpd, method='mps') +gml = ws.genpareto.fit2(Rgpd, method='ml') + +gp.plotecdf() +x = sorted(Rgpd) +plt.hold(True) +plt.plot(x, gml.cdf(x)) +# gm = fitgenpar(Rgpd,'method','mom','plotflag',0); +# plot(x,cdfgenpar(x,gm),'g--') +#gw = fitgenpar(Rgpd,'method','pwm','plotflag',0); +#plot(x,cdfgenpar(x,gw),'g:') +#gml = fitgenpar(Rgpd,'method','ml','plotflag',0); +#plot(x,cdfgenpar(x,gml),'--') +#gmps = fitgenpar(Rgpd,'method','mps','plotflag',0); +#plot(x,cdfgenpar(x,gmps),'r-.') +plt.hold(False) +#wafostamp([],'(ER)') +#disp('Block = 8'),pause(pstate) + +## +# Return values for the GEV distribution +T = np.logspace(1, 5, 10); +#[sT, sTlo, sTup] = invgev(1./T,Y5gev,'lowertail',false,'proflog',true); + +#T = 2:100000; +#k=Y5gev.params(1); mu=Y5gev.params(3); sigma=Y5gev.params(2); +#sT1 = invgev(1./T,Y5gev,'lowertail',false); +#sT=mu + sigma/k*(1-(-log(1-1./T)).^k); +plt.clf() +#plt.semilogx(T,sT,T,sTlo,'r',T,sTup,'r') +#plt.hold(True) +#N = np.r_[1:len(Y5M)] +#Nmax = max(N); +#plot(Nmax./N, sorted(Y5M,reverse=True), '.') +#plt.title('Return values in the GEV model') +#plt.xlabel('Return priod') +#plt.ylabel('Return value') +#plt.grid(True) +#disp('Block = 9'),pause(pstate) + +## Section 5.3 POT-analysis + +# Estimated expected exceedance over level u as function of u. +plt.clf() + +mrl = ws.reslife(Hs,'umin',2,'umax',10,'Nu',200); +mrl.plot() +#wafostamp([],'(ER)') +#disp('Block = 10'),pause(pstate) + +## +# Estimated distribution functions of monthly maxima +#with the POT method (solid), +# fitting a GEV (dashed) and the empirical distribution. + +# POT- method +gpd7 = ws.genpareto.fit2(Hs(Hs>7)-7, method='mps', floc=0) +khat, loc, sigmahat = gpd7.par + +muhat = len(Hs[Hs>7])/(7*3*2) +bhat = sigmahat/muhat**khat +ahat = 7-(bhat-sigmahat)/khat +x = np.linspace(5,15,200); +plt.plot(x,ws.genextreme.cdf(x, khat,bhat,ahat)) +# disp('Block = 11'),pause(pstate) + +## +# Since we have data to compute the monthly maxima mm over +#42 months we can also try to fit a +# GEV distribution directly: +mm = np.zeros((1,41)) +for i in range(41): + mm[i] = max(Hs[((i-1)*14+1):i*14]) + + +gev = ws.genextreme.fit2(mm) + + +plt.hold(True) +gev.plotecdf() + +plt.hold(False) +#wafostamp([],'(ER)') +#disp('Block = 12, Last block'),pause(pstate) diff --git a/wafo/graphutil.py b/wafo/graphutil.py index 0417456..b288dd6 100644 --- a/wafo/graphutil.py +++ b/wafo/graphutil.py @@ -45,7 +45,10 @@ def delete_text_object(gidtxt, figure=None, axis=None, verbose=False): def _delete_gid_objects(handle, gidtxt, verbose): objs = handle.findobj(lmatchfun) - name = handle.__name__ + try: + name = handle.__class__.__name__ + except AttributeError: + name = 'unknown object' msg = "Tried to delete a non-existing {0} from {1}".format(gidtxt, name) for obj in objs: @@ -127,7 +130,10 @@ def cltext(levels, percent=False, n=4, xs=0.036, ys=0.94, zs=0, figure=None, yss = yint[0] + ys * (yint[1] - yint[0]) # delete cltext object if it exists - delete_text_object(_CLTEXT_GID, axis=axis) + try: + delete_text_object(_CLTEXT_GID, axis=axis) + except Exception: + pass charHeight = 1.0 / 33.0 delta_y = charHeight diff --git a/wafo/spectrum/core.py b/wafo/spectrum/core.py index 87f847b..4ecf5dd 100644 --- a/wafo/spectrum/core.py +++ b/wafo/spectrum/core.py @@ -3447,8 +3447,10 @@ class SpecData1D(PlotData): >>> S.bandwidth([0,'eps2',2,3]) array([ 0.73062845, 0.34476034, 0.68277527, 2.90817052]) ''' - m, unused_mtxt = self.moment(nr=4, even=False) + m, unused_mtxt = self.moment(nr=4, even=False) + if isinstance(factors, str): + factors = [factors] fact_dict = dict(alpha=0, eps2=1, eps4=3, qp=3, Qp=3) fact = array([fact_dict.get(idx, idx) for idx in list(factors)], dtype=int)