From 49e5e9d6e40b4b52498345d42eac6596360f793f Mon Sep 17 00:00:00 2001 From: "per.andreas.brodtkorb" Date: Tue, 1 Feb 2011 23:55:49 +0000 Subject: [PATCH] Updated chapter scripts + small fixes elsewhere --- .../src/wafo/doc/tutorial_scripts/chapter2.py | 80 +++++++----------- .../src/wafo/doc/tutorial_scripts/chapter3.py | 84 ++++++++++--------- .../src/wafo/doc/tutorial_scripts/chapter4.py | 4 +- pywafo/src/wafo/objects.py | 15 ++-- pywafo/src/wafo/transform/core.py | 12 ++- 5 files changed, 90 insertions(+), 105 deletions(-) diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py index 28e1456..fe4616a 100644 --- a/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py +++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter2.py @@ -171,9 +171,9 @@ show() #! Test Gaussianity of a stochastic process #!------------------------------------------ #! 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. +#! 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. @@ -182,7 +182,7 @@ 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() @@ -295,49 +295,31 @@ show() ##!wafostamp('','(CR)') #disp('Block = 24'),pause(pstate) # -##!#! Estimated spectrum compared to Torsethaugen spectrum -#clf -#Tp = 1.1; -#H0 = 4*sqrt(spec2mom(S1,1)) -#St = torsethaugen([0:0.01:5],[H0 2*pi/Tp]); -#plotspec(S1) -#hold on -#plotspec(St,'-.') -#axis([0 6 0 0.4]) -##!wafostamp('','(ER)') -#disp('Block = 25'),pause(pstate) -# -##!#! -#clf -#Snorm = St; -#Snorm.S = Snorm.S/sa^2; -#dt = spec2dt(Snorm) -#disp('Block = 26'),pause(pstate) -# -##!#! -#clf -#[Sk Su] = spec2skew(St); -#sa = sqrt(spec2mom(St,1)); -#gh = hermitetr([],[sa sk ku me]); -#Snorm.tr = gh; -#disp('Block = 27'),pause(pstate) -# -##!#! Transformed Gaussian model compared to Gaussian model -#clf -#dt = 0.5; -#ysim_t = spec2sdat(Snorm,240,dt); -#xsim_t = dat2gaus(ysim_t,Snorm.tr); -#disp('Block = 28'),pause(pstate) -# -##!#! Compare -##! In order to compare the Gaussian and non-Gaussian models we need to scale -##! \verb+xsim_t+ #!{\tt xsim$_t$} -##! to have the same first spectral moment as -##! \verb+ysim_t+, #!{\tt ysim$_t$}, Since the process xsim_t has variance one -##! which will be done by the following commands. -#clf -#xsim_t(:,2) = sa*xsim_t(:,2); -#waveplot(xsim_t,ysim_t,5,1,sa,4.5,'r.','b') -##!wafostamp('','(CR)') -#disp('Block = 29, Last block'),pause(pstate) +#! Estimated spectrum compared to Torsethaugen spectrum +#!------------------------------------------------------- + +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)) +S1.plot() +St.plot('-.') +axis([0, 6, 0, 0.4]) +show() + +#! Transformed Gaussian model compared to Gaussian model +#!-------------------------------------------------------- +dt = St.sampling_period() +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]) + +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() diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py index 930ec3f..ec94b60 100644 --- a/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py +++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter3.py @@ -45,49 +45,51 @@ show() #! Extreme waves - model check: the highest and steepest wave #!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ clf() -method = 0; -rate = 8; -[S, H, Ac, At, Tcf, Tcb, z_ind, yn] = ... - dat2steep(xx,rate,method); -disp('Block = 4'), pause(pstate) - -clf -[Smax indS]=max(S) -[Amax indA]=max(Ac) -spwaveplot(yn,[indA indS],'k.') -wafostamp([],'(ER)') -disp('Block = 5'), pause(pstate) - -#!#! Does the highest wave contradict a transformed Gaussian model? -clf -inds1 = (5965:5974)'; #! points to remove -Nsim = 10; -[y1, grec1, g2, test, tobs, mu1o, mu1oStd] = ... - reconstruct(xx,inds1,Nsim); -spwaveplot(y1,indA-10) -hold on -plot(xx(inds1,1),xx(inds1,2),'+') -lamb = 2.; -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]), -wafostamp([],'(ER)') -disp('Block = 6'), -pause(pstate) - -#!#! Expected value (solid) compared to data removed -clf -plot(xx(inds1,1),xx(inds1,2),'+'), hold on -mu = tranproc(mu1o,fliplr(grec1)); -plot(y1(inds1,1), mu), hold off -disp('Block = 7'), pause(pstate) - -#!#! Crest height PDF +S, H = ts.wave_height_steepness(method=0) +indS = S.argmax() +indH = H.argmax() +ts.plot_sp_wave([indH, indS],'k.') +show() + +#! Does the highest wave contradict a transformed Gaussian model? +#!---------------------------------------------------------------- + +# TODO: Fix this + +#clf +#inds1 = (5965:5974)'; #! points to remove +#Nsim = 10; +#[y1, grec1, g2, test, tobs, mu1o, mu1oStd] = ... +# reconstruct(xx,inds1,Nsim); +#spwaveplot(y1,indA-10) +#hold on +#plot(xx(inds1,1),xx(inds1,2),'+') +#lamb = 2.; +#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]), +#wafostamp([],'(ER)') +#disp('Block = 6'), +#pause(pstate) +# +##!#! Expected value (solid) compared to data removed +#clf +#plot(xx(inds1,1),xx(inds1,2),'+'), hold on +#mu = tranproc(mu1o,fliplr(grec1)); +#plot(y1(inds1,1), mu), hold off +#disp('Block = 7'), pause(pstate) + +#! Crest height PDF +#!------------------ #! Transform data so that kde works better -clf +clf() +wave_data = ts.wave_parameters() +Ac = wave_data['Ac'] L2 = 0.6; -plotnorm(Ac.^L2) +import pylab +ws.probplot(Ac**L2, dist='norm', plot=pylab) +show() #!#! #! diff --git a/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py b/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py index 1209b03..a4017f9 100644 --- a/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py +++ b/pywafo/src/wafo/doc/tutorial_scripts/chapter4.py @@ -65,10 +65,10 @@ mM_rfc = tp.cycle_pairs(h=0.3) clf() subplot(122), mM.plot() -title('min-max cycle count') +title('min-max cycle pairs') subplot(121), mM_rfc.plot() -title('Rainflow cycle count') +title('Rainflow filtered cycles') show() #! Min-max and rainflow cycle distributions diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index 9f63c66..5813629 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -63,7 +63,7 @@ class LevelCrossings(WafoData): Member variables ---------------- data : array-like - number of upcrossings + number of upcrossings or upcrossingintensity args : array-like crossing levels @@ -1858,7 +1858,7 @@ class TimeSeries(WafoData): plot = plotbackend.plot subplot = plotbackend.subplot figs = [] - for iz in xrange(nfig): + for unused_iz in xrange(nfig): figs.append(plotbackend.figure()) plotbackend.title('Surface elevation from mean water level (MWL).') for ix in xrange(nsub): @@ -1883,10 +1883,12 @@ class TimeSeries(WafoData): return figs - def plot_sp_wave(self, wave_idx_, tz_idx=None, *args, **kwds): + def plot_sp_wave(self, wave_idx_, *args, **kwds): """ Plot specified wave(s) from timeseries - + + Parameters + ---------- wave_idx : integer vector of indices to waves we want to plot, i.e., wave numbers. tz_idx : integer vector @@ -1907,8 +1909,9 @@ class TimeSeries(WafoData): plot_wave, findtc """ wave_idx = atleast_1d(wave_idx_).flatten() + tz_idx = kwds.pop('tz_idx', None) if tz_idx is None: - tc_ind, tz_idx = findtc(self.data, 0, 'tw') # finding trough to trough waves + unused_tc_ind, tz_idx = findtc(self.data, 0, 'tw') # finding trough to trough waves dw = nonzero(abs(diff(wave_idx)) > 1)[0] Nsub = dw.size + 1 @@ -1932,7 +1935,7 @@ class TimeSeries(WafoData): Nfig = int(ceil(Nsub / 6)) Nsub = min(6, int(ceil(Nsub / Nfig))) figs = [] - for iy in range(Nfig): + for unused_iy in range(Nfig): figs.append(plotbackend.figure()) for ix in range(Nsub): plotbackend.subplot(Nsub, 1, mod(ix, Nsub) + 1) diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py index d05a98b..abebfa6 100644 --- a/pywafo/src/wafo/transform/core.py +++ b/pywafo/src/wafo/transform/core.py @@ -166,19 +166,17 @@ class TrData(WafoData, TrCommon): >>> g.dist2gauss() < 1e-16 True """ - def __init__(self, *args, **kwds): - self.ymean = kwds.pop('ymean', 0e0) - self.ysigma = kwds.pop('ysigma', 1e0) - self.mean = kwds.pop('mean', None) - self.sigma = kwds.pop('sigma', None) - + def __init__(self, *args, **kwds): options = dict(title='Transform', xlab='x', ylab='g(x)', plot_args=['r'], plot_args_children=['g--'],) options.update(**kwds) super(TrData, self).__init__(*args, **options) - + self.ymean = kwds.get('ymean', 0e0) + self.ysigma = kwds.get('ysigma', 1e0) + self.mean = kwds.get('mean', None) + self.sigma = kwds.get('sigma', None) if self.mean is None: #self.mean = np.mean(self.args) #