danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Fixed a bug in LevelCrossings.trdata

Browse Source 
Added plotflag to Plotter_1d and WafoData.
Otherwise small cosmetic fixes
master
per.andreas.brodtkorb 14 years ago
parent 5704928b44
commit 843ddb90c4
7 changed files with 236 additions and 126 deletions
Show Stats Download Patch File Download Diff File

1
pywafo/src/wafo/__init__.py
Unescape Escape View File

@ -2,6 +2,7 @@
from info import __doc__ from info import __doc__
import misc import misc
import data import data
import demos
import kdetools import kdetools
import objects import objects
import spectrum import spectrum

88
pywafo/src/wafo/objects.py
Unescape Escape View File

@ -179,7 +179,8 @@ class LevelCrossings(WafoData):
L0 = vstack((L0, r1 * L0 + sqrt(1 - r2 ** 2) * randn(1))) L0 = vstack((L0, r1 * L0 + sqrt(1 - r2 ** 2) * randn(1)))
#%Simulate the process, starting in L0 #%Simulate the process, starting in L0
lfilter = scipy.signal.lfilter lfilter = scipy.signal.lfilter
L = lfilter(1, [1, a1, a2], e, lfilter([1, a1, a2], 1, L0)) # TODO: lfilter crashes the example
L = lfilter(ones(1), [1, a1, a2], e, lfilter([1, a1, a2], ones(1), L0))
epsilon = 1.01 epsilon = 1.01
min_L = min(L) min_L = min(L)
@ -322,20 +323,20 @@ class LevelCrossings(WafoData):
>>> int(S.tr.dist2gauss()*100) >>> int(S.tr.dist2gauss()*100)
141 141
>>> int(g0emp.dist2gauss()*100) >>> int(g0emp.dist2gauss()*100)
1544 380995
>>> int(g0.dist2gauss()*100) >>> int(g0.dist2gauss()*100)
340 143
>>> int(g1.dist2gauss()*100) >>> int(g1.dist2gauss()*100)
352 162
>>> int(g2.dist2gauss()*100) >>> int(g2.dist2gauss()*100)
365 120
hold on, trplot(g1,g) # Check the fit g0.plot() # Check the fit.
trplot(g2)
See also troptset, dat2tr, trplot, findcross, smooth See also
troptset, dat2tr, trplot, findcross, smooth
NB! the transformated data will be N(0,1) NB! the transformated data will be N(0,1)
Reference Reference
--------- ---------
@ -346,21 +347,13 @@ class LevelCrossings(WafoData):
Vol 10, pp 13--47 Vol 10, pp 13--47
''' '''
# Tested on: Matlab 5.3, 5.2, 5.1
# History:
# by pab 29.12.2000
# based on lc2tr, but the inversion is faster.
# by IR and PJ
if mean is None: if mean is None:
mean = self.mean mean = self.mean
if sigma is None: if sigma is None:
sigma = self.stdev sigma = self.stdev
opt = DotDict(chkder=True, plotflag=False, csm=0.9, gsm=.05, opt = DotDict(chkder=True, plotflag=False, csm=0.9, gsm=.05,
param=(-5, 5, 513), delay=2, linextrap=True, ntr=inf, ne=7, cvar=1, gvar=1) param=(-5, 5, 513), delay=2, linextrap=True, ntr=10000, ne=7, gvar=1)
# If just 'defaults' passed in, return the default options in g
opt.update(options) opt.update(options)
param = opt.param param = opt.param
Ne = opt.ne Ne = opt.ne
@ -381,20 +374,10 @@ class LevelCrossings(WafoData):
else: else:
gvar = interp1d(linspace(0, 1, ng) , opt.gvar, kind='linear')(linspace(0, 1, ncr)) gvar = interp1d(linspace(0, 1, ng) , opt.gvar, kind='linear')(linspace(0, 1, ncr))
ng = len(atleast_1d(opt.cvar))
if ng == 1:
cvar = opt.cvar * ones(ncr)
else:
cvar = interp1d(linspace(0, 1, ng), opt.cvar, kind='linear')(linspace(0, 1, ncr))
uu = linspace(*param) uu = linspace(*param)
g1 = sigma * uu + mean g1 = sigma * uu + mean
g22 = lc2.copy()
if Ner > 0: # Compute correction factors if Ner > 0: # Compute correction factors
cor1 = trapz(lc2[0:Ner + 1], lc1[0:Ner + 1]) cor1 = trapz(lc2[0:Ner + 1], lc1[0:Ner + 1])
cor2 = trapz(lc2[-Ner - 1::], lc1[-Ner - 1::]) cor2 = trapz(lc2[-Ner - 1::], lc1[-Ner - 1::])
@ -402,27 +385,11 @@ class LevelCrossings(WafoData):
cor1 = 0 cor1 = 0
cor2 = 0 cor2 = 0
lc22 = hstack((0, cumtrapz(lc2, lc1) + cor1)) lc22 = hstack((0, cumtrapz(lc2, lc1) + cor1))
lc22 = (lc22 + 0.5) / (lc22[-1] + cor2 + 1) lc22 = (lc22 + 0.5) / (lc22[-1] + cor2 + 1)
lc11 = (lc1 - mean) / sigma lc11 = (lc1 - mean) / sigma
# find the mode lc22 = invnorm(lc22) #- ymean
imin = abs(lc22 - 0.15).argmin()
imax = abs(lc22 - 0.85).argmin()
inde = slice(imin, imax + 1)
lc222 = SmoothSpline(lc11[inde], g22[inde], opt.csm, opt.linextrap, cvar[inde])(lc11[inde])
#tmp = smooth(cros(inde,1),g2(inde,2),opt.csm,cros(inde,1),def,cvar(inde));
imax = lc222.argmax()
u0 = lc22[inde][imax]
#u0 = interp1q(cros(:,2),cros(:,1),.5)
#lc22 = ndtri(lc22) - u0 #
lc22 = invnorm(lc22) - u0
g2 = TrData(lc22.copy(), lc1.copy(), mean=mean, sigma=sigma) g2 = TrData(lc22.copy(), lc1.copy(), mean=mean, sigma=sigma)
g2.setplotter('step') g2.setplotter('step')
@ -432,9 +399,9 @@ class LevelCrossings(WafoData):
# to be linear outside the edges or choosing a lower value for csm2. # to be linear outside the edges or choosing a lower value for csm2.
inds = slice(Ne, ncr - Ne) # indices to points we are smoothing over inds = slice(Ne, ncr - Ne) # indices to points we are smoothing over
scros2 = SmoothSpline(lc11[inds], lc22[inds], opt.gsm, opt.linextrap, gvar[inds])(uu) slc22 = SmoothSpline(lc11[inds], lc22[inds], opt.gsm, opt.linextrap, gvar[inds])(uu)
g = TrData(scros2, g1, mean=mean, sigma=sigma) #*sa; #multiply with stdev g = TrData(slc22.copy(), g1.copy(), mean=mean, sigma=sigma) #*sa; #multiply with stdev
if opt.chkder: if opt.chkder:
for ix in range(5): for ix in range(5):
@ -739,7 +706,7 @@ class TimeSeries(WafoData):
>>> h = rf.plot() >>> h = rf.plot()
>>> S = ts.tospecdata() >>> S = ts.tospecdata()
The default L is set to 68 The default L is set to 325
>>> tp = ts.turning_points() >>> tp = ts.turning_points()
>>> mm = tp.cycle_pairs() >>> mm = tp.cycle_pairs()
@ -1215,29 +1182,20 @@ class TimeSeries(WafoData):
>>> g1, g1emp = ts.trdata(method='m', gvar=0.5 ) # Equal weight on all points >>> g1, g1emp = ts.trdata(method='m', gvar=0.5 ) # Equal weight on all points
>>> g2, g2emp = ts.trdata(method='n', gvar=[3.5, 0.5, 3.5]) # Less weight on the ends >>> g2, g2emp = ts.trdata(method='n', gvar=[3.5, 0.5, 3.5]) # Less weight on the ends
>>> int(S.tr.dist2gauss()*100) >>> int(S.tr.dist2gauss()*100)
593 141
>>> int(g0emp.dist2gauss()*100) >>> int(g0emp.dist2gauss()*100)
439 217949
>>> int(g0.dist2gauss()*100) >>> int(g0.dist2gauss()*100)
432 93
>>> int(g1.dist2gauss()*100) >>> int(g1.dist2gauss()*100)
234 66
>>> int(g2.dist2gauss()*100) >>> int(g2.dist2gauss()*100)
437 84
Hm0 = 7;
S = jonswap([],Hm0); g=ochitr([],[Hm0/4]);
S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4;
xs = spec2sdat(S,2^13);
g0 = dat2tr(xs,[],'plot','iter'); % Monitor the development
g1 = dat2tr(xs,'mnon','gvar', .5 ); % More weight on all points
g2 = dat2tr(xs,'nonl','gvar', [3.5 .5 3.5]); % Less weight on the ends
hold on, trplot(g1,g) % Check the fit
trplot(g2)
See also See also
-------- --------
troptset, lc2tr, cdf2tr, trplot LevelCrossings.trdata
wafo.transform.models
References References
---------- ----------
@ -1259,7 +1217,7 @@ class TimeSeries(WafoData):
# 'cvar',1,'gvar',1,'multip',0); # 'cvar',1,'gvar',1,'multip',0);
opt = DotDict(chkder=True, plotflag=False, csm=.95, gsm=.05, opt = DotDict(chkder=True, plotflag=False, csm=.95, gsm=.05,
param=[-5, 5, 513], delay=2, ntr=inf, linextrap=True, ne=7, cvar=1, gvar=1, param=[-5, 5, 513], delay=2, ntr=1000, linextrap=True, ne=7, cvar=1, gvar=1,
multip=False, crossdef='uM') multip=False, crossdef='uM')
opt.update(**options) opt.update(**options)

65
pywafo/src/wafo/spectrum/models.py
Unescape Escape View File

@ -97,9 +97,10 @@ def _gengamspec(wn, N=5, M=4):
Examples Examples
-------- --------
>>> import wafo.spectrum.models as wsm
>>> import numpy as np >>> import numpy as np
>>> wn = np.linspace(0,4,5) >>> wn = np.linspace(0,4,5)
>>> _gengamspec(wn, N=6, M=2) >>> wsm._gengamspec(wn, N=6, M=2)
array([ 0. , 1.16765216, 0.17309961, 0.02305179, 0.00474686]) array([ 0. , 1.16765216, 0.17309961, 0.02305179, 0.00474686])
See also See also
@ -232,7 +233,8 @@ class Bretschneider(ModelSpectrum):
Examples Examples
-------- --------
>>> S = Bretschneider(Hm0=6.5,Tp=10) >>> import wafo.spectrum.models as wsm
>>> S = wsm.Bretschneider(Hm0=6.5,Tp=10)
>>> S((0,1,2,3)) >>> S((0,1,2,3))
array([ 0. , 1.69350993, 0.06352698, 0.00844783]) array([ 0. , 1.69350993, 0.06352698, 0.00844783])
@ -289,21 +291,22 @@ def jonswap_peakfact(Hm0, Tp):
Examples Examples
-------- --------
>>> import wafo.spectrum.models as wsm
>>> import pylab as plb >>> import pylab as plb
>>> Tp,Hs = plb.meshgrid(range(4,8),range(2,6)) >>> Tp,Hs = plb.meshgrid(range(4,8),range(2,6))
>>> gam = jonswap_peakfact(Hs,Tp) >>> gam = wsm.jonswap_peakfact(Hs,Tp)
>>> Hm0 = plb.linspace(1,20) >>> Hm0 = plb.linspace(1,20)
>>> Tp = Hm0 >>> Tp = Hm0
>>> [T,H] = plb.meshgrid(Tp,Hm0) >>> [T,H] = plb.meshgrid(Tp,Hm0)
>>> gam = jonswap_peakfact(H,T) >>> gam = wsm.jonswap_peakfact(H,T)
>>> v = plb.arange(0,8) >>> v = plb.arange(0,8)
>>> h = plb.contourf(Tp,Hm0,gam,v);h=plb.colorbar() >>> h = plb.contourf(Tp,Hm0,gam,v);h=plb.colorbar()
>>> Hm0 = plb.arange(1,11) >>> Hm0 = plb.arange(1,11)
>>> Tp = plb.linspace(2,16) >>> Tp = plb.linspace(2,16)
>>> T,H = plb.meshgrid(Tp,Hm0) >>> T,H = plb.meshgrid(Tp,Hm0)
>>> gam = jonswap_peakfact(H,T) >>> gam = wsm.jonswap_peakfact(H,T)
>>> h = plb.plot(Tp,gam.T) >>> h = plb.plot(Tp,gam.T)
>>> h = plb.xlabel('Tp [s]') >>> h = plb.xlabel('Tp [s]')
>>> h = plb.ylabel('Peakedness parameter') >>> h = plb.ylabel('Peakedness parameter')
@ -356,8 +359,9 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'):
Example Example
-------- --------
>>> import wafo.spectrum.models as wsm
>>> fetch = 10000; u10 = 10 >>> fetch = 10000; u10 = 10
>>> ss = jonswap_seastate(u10, fetch, output='dict') >>> ss = wsm.jonswap_seastate(u10, fetch, output='dict')
>>> ss >>> ss
{'Ag': 0.016257903375341734, {'Ag': 0.016257903375341734,
'Hm0': 0.51083679198275533, 'Hm0': 0.51083679198275533,
@ -365,13 +369,13 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'):
'gamma': 2.4824142635861119, 'gamma': 2.4824142635861119,
'sigmaA': 0.075317331395172021, 'sigmaA': 0.075317331395172021,
'sigmaB': 0.091912084512251344} 'sigmaB': 0.091912084512251344}
>>> S = Jonswap(**ss) >>> S = wsm.Jonswap(**ss)
>>> S.Hm0 >>> S.Hm0
0.51083679198275533 0.51083679198275533
# Alternatively # Alternatively
>>> ss1 = jonswap_seastate(u10, fetch, output='list') >>> ss1 = wsm.jonswap_seastate(u10, fetch, output='list')
>>> S1 = Jonswap(*ss1) >>> S1 = wsm.Jonswap(*ss1)
>>> S1.Hm0 >>> S1.Hm0
0.51083679198275533 0.51083679198275533
@ -484,12 +488,14 @@ class Jonswap(ModelSpectrum):
Examples Examples
--------- ---------
>>> import pylab as plb >>> import pylab as plb
>>> S = Jonswap(Hm0=7, Tp=11,gamma=1) >>> import wafo.spectrum.models as wsm
>>> S = wsm.Jonswap(Hm0=7, Tp=11,gamma=1)
>>> w = plb.linspace(0,5) >>> w = plb.linspace(0,5)
>>> h = plb.plot(w,S(w)) >>> h = plb.plot(w,S(w))
>>> S2 = Bretschneider(Hm0=7, Tp=11) >>> S2 = wsm.Bretschneider(Hm0=7, Tp=11)
>>> assert(all(abs(S(w)-S2(w))<1.e-7),'JONSWAP with gamma=1 differs from Bretscneider, should be equal!') >>> all(abs(S(w)-S2(w))<1.e-7)
True
>>> plb.close('all') >>> plb.close('all')
See also See also
@ -652,10 +658,10 @@ def phi1(wi, h, g=9.81):
Example: Example:
------- -------
Transform a JONSWAP spectrum to a spectrum for waterdepth = 30 m Transform a JONSWAP spectrum to a spectrum for waterdepth = 30 m
>>> import wafo.spectrum.models as wsm
>>> S = Jonswap() >>> S = wsm.Jonswap()
>>> w = range(3.0) >>> w = range(3.0)
>>> S(w)*phi1(w,30.0) >>> S(w)*wsm.phi1(w,30.0)
array([ 0. , 1.0358056 , 0.03796281]) array([ 0. , 1.0358056 , 0.03796281])
@ -723,9 +729,10 @@ class Tmaspec(Jonswap):
Example Example
-------- --------
>>> import wafo.spectrum.models as wsm
>>> import pylab as plb >>> import pylab as plb
>>> w = plb.linspace(0,2.5) >>> w = plb.linspace(0,2.5)
>>> S = Tmaspec(h=10,gamma=1) # Bretschneider spectrum Hm0=7, Tp=11 >>> S = wsm.Tmaspec(h=10,gamma=1) # Bretschneider spectrum Hm0=7, Tp=11
>>> o=plb.plot(w,S(w)) >>> o=plb.plot(w,S(w))
>>> o=plb.plot(w,S(w,h=21)) >>> o=plb.plot(w,S(w,h=21))
>>> o=plb.plot(w,S(w,h=42)) >>> o=plb.plot(w,S(w,h=42))
@ -812,9 +819,10 @@ class Torsethaugen(ModelSpectrum):
Example Example
------- -------
>>> import wafo.spectrum.models as wsm
>>> import pylab as plb >>> import pylab as plb
>>> w = plb.linspace(0,4) >>> w = plb.linspace(0,4)
>>> S = Torsethaugen(Hm0=6, Tp=8) >>> S = wsm.Torsethaugen(Hm0=6, Tp=8)
>>> h=plb.plot(w,S(w),w,S.wind(w),w,S.swell(w)) >>> h=plb.plot(w,S(w),w,S.wind(w),w,S.swell(w))
See also See also
@ -1034,7 +1042,8 @@ class McCormick(Bretschneider):
Example: Example:
-------- --------
>>> S = McCormick(Hm0=6.5,Tp=10) >>> import wafo.spectrum.models as wsm
>>> S = wsm.McCormick(Hm0=6.5,Tp=10)
>>> S(range(4)) >>> S(range(4))
array([ 0. , 1.87865908, 0.15050447, 0.02994663]) array([ 0. , 1.87865908, 0.15050447, 0.02994663])
@ -1107,7 +1116,8 @@ class OchiHubble(ModelSpectrum):
Examples Examples
-------- --------
>>> S = OchiHubble(par=2) >>> import wafo.spectrum.models as wsm
>>> S = wsm.OchiHubble(par=2)
>>> S(range(4)) >>> S(range(4))
array([ 0. , 0.90155636, 0.04185445, 0.00583207]) array([ 0. , 0.90155636, 0.04185445, 0.00583207])
@ -1239,7 +1249,8 @@ class Wallop(Bretschneider):
Example: Example:
-------- --------
>>> S = Wallop(Hm0=6.5, Tp=10) >>> import wafo.spectrum.models as wsm
>>> S = wsm.Wallop(Hm0=6.5, Tp=10)
>>> S(range(4)) >>> S(range(4))
array([ 0.00000000e+00, 9.36921871e-01, 2.76991078e-03, array([ 0.00000000e+00, 9.36921871e-01, 2.76991078e-03,
7.72996150e-05]) 7.72996150e-05])
@ -1360,8 +1371,9 @@ class Spreading(object):
Examples Examples
-------- --------
>>> import wafo.spectrum.models as wsm
>>> import pylab as plb >>> import pylab as plb
>>> D = Spreading('cos2s',s_a=10.0) >>> D = wsm.Spreading('cos2s',s_a=10.0)
>>> w = plb.linspace(0,3,257) >>> w = plb.linspace(0,3,257)
>>> theta = plb.linspace(-pi,pi,129) >>> theta = plb.linspace(-pi,pi,129)
@ -1369,13 +1381,13 @@ class Spreading(object):
# Make frequency dependent direction # Make frequency dependent direction
>>> theta0 = lambda w: w*plb.pi/6.0 >>> theta0 = lambda w: w*plb.pi/6.0
>>> D2 = Spreading('cos2s',theta0=theta0) >>> D2 = wsm.Spreading('cos2s',theta0=theta0)
>>> t = plb.contour(D2(theta,w)[0]) >>> t = plb.contour(D2(theta,w)[0])
# Plot all spreading functions # Plot all spreading functions
>>> alltypes = ('cos2s','box','mises','poisson','sech2','wrap_norm') >>> alltypes = ('cos2s','box','mises','poisson','sech2','wrap_norm')
>>> for ix in range(len(alltypes)): >>> for ix in range(len(alltypes)):
... D3 = Spreading(alltypes[ix]) ... D3 = wsm.Spreading(alltypes[ix])
... t = plb.figure(ix) ... t = plb.figure(ix)
... t = plb.contour(D3(theta,w)[0]) ... t = plb.contour(D3(theta,w)[0])
... t = plb.title(alltypes[ix]) ... t = plb.title(alltypes[ix])
@ -1927,10 +1939,11 @@ class Spreading(object):
Example Example
------- -------
>>> S = Jonswap().tospecdata() >>> import wafo.spectrum.models as wsm
>>> D = Spreading('cos2s') >>> S = wsm.Jonswap().tospecdata()
>>> D = wsm.Spreading('cos2s')
>>> SD = D.tospecdata2d(S) >>> SD = D.tospecdata2d(S)
>>> SD.plot() >>> h = SD.plot()
See also spreading, rotspec, jonswap, torsethaugen See also spreading, rotspec, jonswap, torsethaugen
''' '''

5
pywafo/src/wafo/stats/estimation.py
Unescape Escape View File

@ -837,9 +837,10 @@ class FitDistribution(rv_frozen):
plotbackend.subplot(2, 2, 2) plotbackend.subplot(2, 2, 2)
self.plotepdf() self.plotepdf()
plotbackend.subplot(2, 2, 3) plotbackend.subplot(2, 2, 3)
self.plotresprb()
plotbackend.subplot(2, 2, 4)
self.plotresq() self.plotresq()
plotbackend.subplot(2, 2, 4)
self.plotresprb()
fixstr = '' fixstr = ''
if not self.par_fix == None: if not self.par_fix == None:
numfix = len(self.i_fixed) numfix = len(self.i_fixed)

5
pywafo/src/wafo/transform/core.py
Unescape Escape View File

@ -145,15 +145,16 @@ class TrData(WafoData, TrCommon):
------- -------
Construct a linear transformation model Construct a linear transformation model
>>> import numpy as np >>> import numpy as np
>>> import wafo.transorm as wt
>>> sigma = 5; mean = 1 >>> sigma = 5; mean = 1
>>> u = np.linspace(-5,5); x = sigma*u+mean; y = u >>> u = np.linspace(-5,5); x = sigma*u+mean; y = u
>>> g = TrData(y,x) >>> g = wt.TrData(y,x)
>>> g.mean >>> g.mean
array([ 1.]) array([ 1.])
>>> g.sigma >>> g.sigma
array([ 5.]) array([ 5.])
>>> g = TrData(y,x,mean=1,sigma=5) >>> g = wt.TrData(y,x,mean=1,sigma=5)
>>> g.mean >>> g.mean
1 1
>>> g.sigma >>> g.sigma

28
pywafo/src/wafo/transform/models.py
Unescape Escape View File

@ -26,19 +26,20 @@ from core import TrCommon, TrData
__all__ = ['TrHermite', 'TrLinear', 'TrOchi'] __all__ = ['TrHermite', 'TrLinear', 'TrOchi']
_example = ''' _example = '''
>>> import numpy as np
>>> import wafo.spectrum.models as sm
>>> import wafo.transform.models as tm
>>> std = 7./4 >>> std = 7./4
>>> g = <generic>(sigma=std, ysigma=std) >>> g = tm.<generic>(sigma=std, ysigma=std)
Simulate a Transformed Gaussian process: Simulate a Transformed Gaussian process:
>>> import numpy as np
>>> import wafo.spectrum.models as sm
>>> Sj = sm.Jonswap(Hm0=4*std, Tp=11) >>> Sj = sm.Jonswap(Hm0=4*std, Tp=11)
>>> w = np.linspace(0,4,256) >>> w = np.linspace(0,4,256)
>>> S = Sj.tospecdata(w) #Make spectrum object from numerical values >>> S = Sj.tospecdata(w) # Make spectrum object from numerical values
>>> ys = S.sim(ns=15000) # Simulated in the Gaussian world >>> ys = S.sim(ns=15000) # Simulated in the Gaussian world
>>> me, va, sk, ku = S.stats_nl(moments='mvsk') >>> me, va, sk, ku = S.stats_nl(moments='mvsk')
>>> g2 = <generic>(mean=me, var=va, skew=sk, kurt=ku, ysigma=std) >>> g2 = tm.<generic>(mean=me, var=va, skew=sk, kurt=ku, ysigma=std)
>>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world >>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world
''' '''
class TrCommon2(TrCommon): class TrCommon2(TrCommon):
@ -111,7 +112,9 @@ class TrHermite(TrCommon2):
-------- --------
""" + _example.replace('<generic>', 'TrHermite') + """ """ + _example.replace('<generic>', 'TrHermite') + """
>>> g.dist2gauss() >>> g.dist2gauss()
3.9858776379926808 0.88230868748851499
>>> g2.dist2gauss()
1.1411663205144991
See also See also
-------- --------
@ -337,6 +340,8 @@ class TrLinear(TrCommon2):
""" + _example.replace('<generic>', 'TrLinear') + """ """ + _example.replace('<generic>', 'TrLinear') + """
>>> g.dist2gauss() >>> g.dist2gauss()
0.0 0.0
>>> g2.dist2gauss()
3.8594770921678001e-31
See also See also
-------- --------
@ -399,7 +404,9 @@ class TrOchi(TrCommon2):
------- -------
""" + _example.replace('<generic>', 'TrOchi') + """ """ + _example.replace('<generic>', 'TrOchi') + """
>>> g.dist2gauss() >>> g.dist2gauss()
5.9322684525265501 1.410698801056657
>>> g2.dist2gauss()
1.988807188766706
See also See also
-------- --------
@ -486,7 +493,8 @@ class TrOchi(TrCommon2):
mean = self.mean mean = self.mean
sigma = self.sigma sigma = self.sigma
xn = atleast_1d(x) xn = atleast_1d(x)
xn = (xn - mean) / sigma shape0 = xn.shape
xn = (xn.ravel() - mean) / sigma
igp, = where(0 <= xn) igp, = where(0 <= xn)
igm, = where(xn < 0) igm, = where(xn < 0)
@ -498,6 +506,7 @@ class TrOchi(TrCommon2):
if gb != 0: if gb != 0:
np.put(g, igm, (-expm1(-gb * xn[igm])) / gb) np.put(g, igm, (-expm1(-gb * xn[igm])) / gb)
g.shape = shape0
return (g - mean2) * self.ysigma / sigma2 + self.ymean return (g - mean2) * self.ysigma / sigma2 + self.ymean
def _gauss2dat(self, y, *yi): def _gauss2dat(self, y, *yi):
@ -509,7 +518,7 @@ class TrOchi(TrCommon2):
sigma = self.sigma sigma = self.sigma
yn = (atleast_1d(y) - self.ymean) / self.ysigma yn = (atleast_1d(y) - self.ymean) / self.ysigma
xn = sigma2 * yn + mean2 xn = sigma2 * yn.ravel() + mean2
igp, = where(0 <= xn) igp, = where(0 <= xn)
igm, = where(xn < 0) igm, = where(xn < 0)
@ -520,6 +529,7 @@ class TrOchi(TrCommon2):
if gb != 0: if gb != 0:
np.put(xn, igm, -log1p(-gb * xn[igm]) / gb) np.put(xn, igm, -log1p(-gb * xn[igm]) / gb)
xn.shape = yn.shape
return sigma * xn + mean return sigma * xn + mean
def main(): def main():

158
pywafo/src/wafo/wafodata.py
Unescape Escape View File

@ -2,6 +2,7 @@ import warnings
from plotbackend import plotbackend from plotbackend import plotbackend
from time import gmtime, strftime from time import gmtime, strftime
import numpy as np import numpy as np
from scipy.integrate.quadrature import cumtrapz #@UnresolvedImport
__all__ = ['WafoData', 'AxisLabels'] __all__ = ['WafoData', 'AxisLabels']
@ -66,19 +67,21 @@ class WafoData(object):
self.data = data self.data = data
self.args = args self.args = args
self.date = now() self.date = now()
self.plotter = None self.plotter = kwds.pop('plotter', None)
self.children = None self.children = None
self.plot_args_children = kwds.get('plot_args_children',[]) self.plot_args_children = kwds.pop('plot_args_children',[])
self.plot_kwds_children = kwds.get('plot_kwds_children',{}) self.plot_kwds_children = kwds.pop('plot_kwds_children',{})
self.plot_args = kwds.get('plot_args',[]) self.plot_args = kwds.pop('plot_args',[])
self.plot_kwds = kwds.get('plot_kwds',{}) self.plot_kwds = kwds.pop('plot_kwds',{})
self.labels = AxisLabels(**kwds) self.labels = AxisLabels(**kwds)
self.setplotter() if not self.plotter:
self.setplotter()
def plot(self, *args, **kwds): def plot(self, *args, **kwds):
tmp = None tmp = None
if self.children != None: plotflag = kwds.get('plotflag', None)
if not plotflag and self.children != None:
plotbackend.hold('on') plotbackend.hold('on')
tmp = [] tmp = []
child_args = args if len(args) else tuple(self.plot_args_children) child_args = args if len(args) else tuple(self.plot_args_children)
@ -106,12 +109,10 @@ class WafoData(object):
newcopy.__dict__.update(self.__dict__) newcopy.__dict__.update(self.__dict__)
return newcopy return newcopy
def setplotter(self, plotmethod=None): def setplotter(self, plotmethod=None):
''' '''
Set plotter based on the data type data_1d, data_2d, data_3d or data_nd Set plotter based on the data type data_1d, data_2d, data_3d or data_nd
''' '''
if isinstance(self.args, (list, tuple)): # Multidimensional data if isinstance(self.args, (list, tuple)): # Multidimensional data
ndim = len(self.args) ndim = len(self.args)
if ndim < 2: if ndim < 2:
@ -140,10 +141,6 @@ class AxisLabels:
newcopy = empty_copy(self) newcopy = empty_copy(self)
newcopy.__dict__.update(self.__dict__) newcopy.__dict__.update(self.__dict__)
return newcopy return newcopy
#lbkwds = self.labels.__dict__.copy()
#labels = AxisLabels(**lbkwds)
#return labels
def labelfig(self): def labelfig(self):
try: try:
h1 = plotbackend.title(self.title) h1 = plotbackend.title(self.title)
@ -178,7 +175,6 @@ class Plotter_1d(object):
plotmethod = 'plot' plotmethod = 'plot'
self.plotbackend = plotbackend self.plotbackend = plotbackend
try: try:
#self.plotfun = plotbackend.__dict__[plotmethod]
self.plotfun = getattr(plotbackend, plotmethod) self.plotfun = getattr(plotbackend, plotmethod)
except: except:
pass pass
@ -190,10 +186,137 @@ class Plotter_1d(object):
args1 = tuple((wdata.args)) + (wdata.data,) + args args1 = tuple((wdata.args)) + (wdata.data,) + args
else: else:
args1 = tuple((wdata.args,)) + (wdata.data,) + args args1 = tuple((wdata.args,)) + (wdata.data,) + args
h1 = self.plotfun(*args1, **kwds) plotflag = kwds.pop('plotflag', None)
if plotflag:
h1 = self._plot(plotflag, *args1, **kwds)
else:
h1 = self.plotfun(*args1, **kwds)
h2 = wdata.labels.labelfig() h2 = wdata.labels.labelfig()
return h1, h2 return h1, h2
def _plot(self, plotflag, *args1, **kwds):
x = args1[0]
data = transformdata(x,args1[1], plotflag)
dataCI = ()
h1 = plot1d(x,data, dataCI, plotflag, *args1[2:], **kwds)
return h1
def plot1d(args,data,dataCI,plotflag,*varargin, **kwds):
plottype = np.mod(plotflag,10)
if plottype==0: # % No plotting
return []
elif plottype==1:
H = plotbackend.plot(args,data,*varargin, **kwds)
elif plottype==2:
H = plotbackend.step(args,data,*varargin, **kwds)
elif plottype==3:
H = plotbackend.stem(args,data,*varargin, **kwds)
elif plottype==4:
H = plotbackend.errorbar(args,data,dataCI[:,0]-data,dataCI[:,1]-data,*varargin, **kwds)
elif plottype==5:
H = plotbackend.bar(args,data,*varargin, **kwds)
elif plottype==6:
level = 0
if np.isfinite(level):
H = plotbackend.fill_between(args,data,level,*varargin, **kwds);
else:
H = plotbackend.fill_between(args,data,*varargin, **kwds);
scale = plotscale(plotflag);
logXscale = any(scale=='x');
logYscale = any(scale=='y');
logZscale = any(scale=='z');
ax = plotbackend.gca()
if logXscale:
plotbackend.setp(ax,xscale='log')
if logYscale:
plotbackend.setp(ax,yscale='log')
if logZscale:
plotbackend.setp(ax,zscale='log')
transFlag = np.mod(plotflag//10,10)
logScale = logXscale or logYscale or logZscale
if logScale or (transFlag ==5 and not logScale):
ax = list(plotbackend.axis())
fmax1 = data.max()
if transFlag==5 and not logScale:
ax[3] = 11*np.log10(fmax1)
ax[2] = ax[3]-40
else:
ax[3] = 1.15*fmax1;
ax[2] = ax[3]*1e-4;
plotbackend.axis(ax)
if dataCI and plottype < 3:
plotbackend.hold('on')
plot1d(args,dataCI,(),plotflag,'r--');
return H
def plotscale(plotflag):
'''
Return plotscale from plotflag
CALL scale = plotscale(plotflag)
plotflag = integer defining plotscale.
Let scaleId = floor(plotflag/100).
If scaleId < 8 then:
0 'linear' : Linear scale on all axes.
1 'xlog' : Log scale on x-axis.
2 'ylog' : Log scale on y-axis.
3 'xylog' : Log scale on xy-axis.
4 'zlog' : Log scale on z-axis.
5 'xzlog' : Log scale on xz-axis.
6 'yzlog' : Log scale on yz-axis.
7 'xyzlog' : Log scale on xyz-axis.
otherwise
if (mod(scaleId,10)>0) : Log scale on x-axis.
if (mod(floor(scaleId/10),10)>0) : Log scale on y-axis.
if (mod(floor(scaleId/100),10)>0) : Log scale on z-axis.
scale = string defining plotscale valid options are:
'linear', 'xlog', 'ylog', 'xylog', 'zlog', 'xzlog',
'yzlog', 'xyzlog'
Example
plotscale(100) % xlog
plotscale(200) % xlog
plotscale(1000) % ylog
See also plotscale
'''
scaleId = plotflag//100
if scaleId>7:
logXscaleId = np.mod(scaleId,10)>0
logYscaleId = (np.mod(scaleId//10, 10)>0)*2
logZscaleId = (np.mod(scaleId//100, 10)>0)*4
scaleId = logYscaleId +logXscaleId+logZscaleId
scales = ['linear','xlog','ylog','xylog','zlog','xzlog','yzlog','xyzlog']
return scales[scaleId]
def transformdata(x,f,plotflag):
transFlag = np.mod(plotflag//10,10)
if transFlag==0:
data = f
elif transFlag==1:
data = 1-f
elif transFlag==2:
data = cumtrapz(f,x)
elif transFlag==3:
data = 1-cumtrapz(f,x)
if transFlag in (4,5):
if transFlag==4:
data = -np.log1p(-cumtrapz(f,x))
else:
if any(f<0):
raise ValueError('Invalid plotflag: Data or dataCI is negative, but must be positive')
data = 10*np.log10(f)
return data
class Plotter_2d(Plotter_1d): class Plotter_2d(Plotter_1d):
""" """
Parameters Parameters
@ -201,6 +324,7 @@ class Plotter_2d(Plotter_1d):
plotmethod : string plotmethod : string
defining type of plot. Options are: defining type of plot. Options are:
contour (default) contour (default)
contourf
mesh mesh
surf surf
""" """
@ -209,7 +333,9 @@ class Plotter_2d(Plotter_1d):
if plotmethod is None: if plotmethod is None:
plotmethod = 'contour' plotmethod = 'contour'
super(Plotter_2d, self).__init__(plotmethod) super(Plotter_2d, self).__init__(plotmethod)
#self.plotfun = plotbackend.__dict__[plotmethod]
def _plot(self,*args,**kwds):
pass
def main(): def main():

Write Preview
Loading…
Cancel
Save