From 843ddb90c46509766d10b9eb188c5dd9d36eb69a Mon Sep 17 00:00:00 2001 From: "per.andreas.brodtkorb" Date: Fri, 28 Jan 2011 15:04:29 +0000 Subject: [PATCH] Fixed a bug in LevelCrossings.trdata Added plotflag to Plotter_1d and WafoData. Otherwise small cosmetic fixes --- pywafo/src/wafo/__init__.py | 3 +- pywafo/src/wafo/objects.py | 96 +++++------------ pywafo/src/wafo/spectrum/models.py | 65 +++++++----- pywafo/src/wafo/stats/estimation.py | 5 +- pywafo/src/wafo/transform/core.py | 5 +- pywafo/src/wafo/transform/models.py | 30 ++++-- pywafo/src/wafo/wafodata.py | 158 +++++++++++++++++++++++++--- 7 files changed, 236 insertions(+), 126 deletions(-) diff --git a/pywafo/src/wafo/__init__.py b/pywafo/src/wafo/__init__.py index 876ea6f..bb654e1 100644 --- a/pywafo/src/wafo/__init__.py +++ b/pywafo/src/wafo/__init__.py @@ -1,7 +1,8 @@ from info import __doc__ import misc -import data +import data +import demos import kdetools import objects import spectrum diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index dc40766..85d1153 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -179,7 +179,8 @@ class LevelCrossings(WafoData): L0 = vstack((L0, r1 * L0 + sqrt(1 - r2 ** 2) * randn(1))) #%Simulate the process, starting in L0 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 min_L = min(L) @@ -322,20 +323,20 @@ class LevelCrossings(WafoData): >>> int(S.tr.dist2gauss()*100) 141 >>> int(g0emp.dist2gauss()*100) - 1544 + 380995 >>> int(g0.dist2gauss()*100) - 340 + 143 >>> int(g1.dist2gauss()*100) - 352 + 162 >>> int(g2.dist2gauss()*100) - 365 + 120 - hold on, trplot(g1,g) # Check the fit - trplot(g2) - - See also troptset, dat2tr, trplot, findcross, smooth + g0.plot() # Check the fit. + + 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 --------- @@ -345,22 +346,14 @@ class LevelCrossings(WafoData): Marine structures, Design, Construction and Safety, 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: mean = self.mean if sigma is None: sigma = self.stdev 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) - # If just 'defaults' passed in, return the default options in g - + param=(-5, 5, 513), delay=2, linextrap=True, ntr=10000, ne=7, gvar=1) opt.update(options) param = opt.param Ne = opt.ne @@ -380,21 +373,11 @@ class LevelCrossings(WafoData): gvar = opt.gvar * ones(ncr) else: 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) - g1 = sigma * uu + mean - g22 = lc2.copy() - if Ner > 0: # Compute correction factors cor1 = trapz(lc2[0:Ner + 1], lc1[0:Ner + 1]) cor2 = trapz(lc2[-Ner - 1::], lc1[-Ner - 1::]) @@ -402,27 +385,11 @@ class LevelCrossings(WafoData): cor1 = 0 cor2 = 0 - lc22 = hstack((0, cumtrapz(lc2, lc1) + cor1)) lc22 = (lc22 + 0.5) / (lc22[-1] + cor2 + 1) lc11 = (lc1 - mean) / sigma - # find the mode - 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 + lc22 = invnorm(lc22) #- ymean g2 = TrData(lc22.copy(), lc1.copy(), mean=mean, sigma=sigma) g2.setplotter('step') @@ -432,9 +399,9 @@ class LevelCrossings(WafoData): # 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 - 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: for ix in range(5): @@ -739,7 +706,7 @@ class TimeSeries(WafoData): >>> h = rf.plot() >>> S = ts.tospecdata() - The default L is set to 68 + The default L is set to 325 >>> tp = ts.turning_points() >>> 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 >>> g2, g2emp = ts.trdata(method='n', gvar=[3.5, 0.5, 3.5]) # Less weight on the ends >>> int(S.tr.dist2gauss()*100) - 593 + 141 >>> int(g0emp.dist2gauss()*100) - 439 + 217949 >>> int(g0.dist2gauss()*100) - 432 + 93 >>> int(g1.dist2gauss()*100) - 234 + 66 >>> int(g2.dist2gauss()*100) - 437 - - 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) - + 84 + See also -------- - troptset, lc2tr, cdf2tr, trplot + LevelCrossings.trdata + wafo.transform.models References ---------- @@ -1259,7 +1217,7 @@ class TimeSeries(WafoData): # 'cvar',1,'gvar',1,'multip',0); 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') opt.update(**options) diff --git a/pywafo/src/wafo/spectrum/models.py b/pywafo/src/wafo/spectrum/models.py index 5bfa945..75e0858 100644 --- a/pywafo/src/wafo/spectrum/models.py +++ b/pywafo/src/wafo/spectrum/models.py @@ -97,9 +97,10 @@ def _gengamspec(wn, N=5, M=4): Examples -------- + >>> import wafo.spectrum.models as wsm >>> import numpy as np >>> 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]) See also @@ -232,7 +233,8 @@ class Bretschneider(ModelSpectrum): 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)) array([ 0. , 1.69350993, 0.06352698, 0.00844783]) @@ -289,21 +291,22 @@ def jonswap_peakfact(Hm0, Tp): Examples -------- + >>> import wafo.spectrum.models as wsm >>> import pylab as plb >>> 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) >>> Tp = Hm0 >>> [T,H] = plb.meshgrid(Tp,Hm0) - >>> gam = jonswap_peakfact(H,T) + >>> gam = wsm.jonswap_peakfact(H,T) >>> v = plb.arange(0,8) >>> h = plb.contourf(Tp,Hm0,gam,v);h=plb.colorbar() >>> Hm0 = plb.arange(1,11) >>> Tp = plb.linspace(2,16) >>> 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.xlabel('Tp [s]') >>> h = plb.ylabel('Peakedness parameter') @@ -356,8 +359,9 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'): Example -------- + >>> import wafo.spectrum.models as wsm >>> fetch = 10000; u10 = 10 - >>> ss = jonswap_seastate(u10, fetch, output='dict') + >>> ss = wsm.jonswap_seastate(u10, fetch, output='dict') >>> ss {'Ag': 0.016257903375341734, 'Hm0': 0.51083679198275533, @@ -365,13 +369,13 @@ def jonswap_seastate(u10, fetch=150000., method='lewis', g=9.81, output='dict'): 'gamma': 2.4824142635861119, 'sigmaA': 0.075317331395172021, 'sigmaB': 0.091912084512251344} - >>> S = Jonswap(**ss) + >>> S = wsm.Jonswap(**ss) >>> S.Hm0 0.51083679198275533 # Alternatively - >>> ss1 = jonswap_seastate(u10, fetch, output='list') - >>> S1 = Jonswap(*ss1) + >>> ss1 = wsm.jonswap_seastate(u10, fetch, output='list') + >>> S1 = wsm.Jonswap(*ss1) >>> S1.Hm0 0.51083679198275533 @@ -484,12 +488,14 @@ class Jonswap(ModelSpectrum): Examples --------- >>> 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) >>> h = plb.plot(w,S(w)) - >>> S2 = 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!') + >>> S2 = wsm.Bretschneider(Hm0=7, Tp=11) + >>> all(abs(S(w)-S2(w))<1.e-7) + True >>> plb.close('all') See also @@ -652,10 +658,10 @@ def phi1(wi, h, g=9.81): Example: ------- Transform a JONSWAP spectrum to a spectrum for waterdepth = 30 m - - >>> S = Jonswap() + >>> import wafo.spectrum.models as wsm + >>> S = wsm.Jonswap() >>> w = range(3.0) - >>> S(w)*phi1(w,30.0) + >>> S(w)*wsm.phi1(w,30.0) array([ 0. , 1.0358056 , 0.03796281]) @@ -723,9 +729,10 @@ class Tmaspec(Jonswap): Example -------- + >>> import wafo.spectrum.models as wsm >>> import pylab as plb >>> 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,h=21)) >>> o=plb.plot(w,S(w,h=42)) @@ -812,9 +819,10 @@ class Torsethaugen(ModelSpectrum): Example ------- + >>> import wafo.spectrum.models as wsm >>> import pylab as plb >>> 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)) See also @@ -1034,7 +1042,8 @@ class McCormick(Bretschneider): 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)) array([ 0. , 1.87865908, 0.15050447, 0.02994663]) @@ -1107,7 +1116,8 @@ class OchiHubble(ModelSpectrum): Examples -------- - >>> S = OchiHubble(par=2) + >>> import wafo.spectrum.models as wsm + >>> S = wsm.OchiHubble(par=2) >>> S(range(4)) array([ 0. , 0.90155636, 0.04185445, 0.00583207]) @@ -1239,7 +1249,8 @@ class Wallop(Bretschneider): 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)) array([ 0.00000000e+00, 9.36921871e-01, 2.76991078e-03, 7.72996150e-05]) @@ -1360,8 +1371,9 @@ class Spreading(object): Examples -------- + >>> import wafo.spectrum.models as wsm >>> 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) >>> theta = plb.linspace(-pi,pi,129) @@ -1369,13 +1381,13 @@ class Spreading(object): # Make frequency dependent direction >>> 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]) # Plot all spreading functions >>> alltypes = ('cos2s','box','mises','poisson','sech2','wrap_norm') >>> for ix in range(len(alltypes)): - ... D3 = Spreading(alltypes[ix]) + ... D3 = wsm.Spreading(alltypes[ix]) ... t = plb.figure(ix) ... t = plb.contour(D3(theta,w)[0]) ... t = plb.title(alltypes[ix]) @@ -1927,10 +1939,11 @@ class Spreading(object): Example ------- - >>> S = Jonswap().tospecdata() - >>> D = Spreading('cos2s') + >>> import wafo.spectrum.models as wsm + >>> S = wsm.Jonswap().tospecdata() + >>> D = wsm.Spreading('cos2s') >>> SD = D.tospecdata2d(S) - >>> SD.plot() + >>> h = SD.plot() See also spreading, rotspec, jonswap, torsethaugen ''' diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py index 4b4988d..55faf42 100644 --- a/pywafo/src/wafo/stats/estimation.py +++ b/pywafo/src/wafo/stats/estimation.py @@ -837,9 +837,10 @@ class FitDistribution(rv_frozen): plotbackend.subplot(2, 2, 2) self.plotepdf() plotbackend.subplot(2, 2, 3) - self.plotresprb() - plotbackend.subplot(2, 2, 4) self.plotresq() + plotbackend.subplot(2, 2, 4) + self.plotresprb() + fixstr = '' if not self.par_fix == None: numfix = len(self.i_fixed) diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py index 392856e..d05a98b 100644 --- a/pywafo/src/wafo/transform/core.py +++ b/pywafo/src/wafo/transform/core.py @@ -145,15 +145,16 @@ class TrData(WafoData, TrCommon): ------- Construct a linear transformation model >>> import numpy as np + >>> import wafo.transorm as wt >>> sigma = 5; mean = 1 >>> u = np.linspace(-5,5); x = sigma*u+mean; y = u - >>> g = TrData(y,x) + >>> g = wt.TrData(y,x) >>> g.mean array([ 1.]) >>> g.sigma array([ 5.]) - >>> g = TrData(y,x,mean=1,sigma=5) + >>> g = wt.TrData(y,x,mean=1,sigma=5) >>> g.mean 1 >>> g.sigma diff --git a/pywafo/src/wafo/transform/models.py b/pywafo/src/wafo/transform/models.py index 815f557..d59afd5 100644 --- a/pywafo/src/wafo/transform/models.py +++ b/pywafo/src/wafo/transform/models.py @@ -26,19 +26,20 @@ from core import TrCommon, TrData __all__ = ['TrHermite', 'TrLinear', 'TrOchi'] _example = ''' + >>> import numpy as np + >>> import wafo.spectrum.models as sm + >>> import wafo.transform.models as tm >>> std = 7./4 - >>> g = (sigma=std, ysigma=std) + >>> g = tm.(sigma=std, ysigma=std) Simulate a Transformed Gaussian process: - >>> import numpy as np - >>> import wafo.spectrum.models as sm >>> Sj = sm.Jonswap(Hm0=4*std, Tp=11) >>> 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 >>> me, va, sk, ku = S.stats_nl(moments='mvsk') - >>> g2 = (mean=me, var=va, skew=sk, kurt=ku, ysigma=std) + >>> g2 = tm.(mean=me, var=va, skew=sk, kurt=ku, ysigma=std) >>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world ''' class TrCommon2(TrCommon): @@ -111,7 +112,9 @@ class TrHermite(TrCommon2): -------- """ + _example.replace('', 'TrHermite') + """ >>> g.dist2gauss() - 3.9858776379926808 + 0.88230868748851499 + >>> g2.dist2gauss() + 1.1411663205144991 See also -------- @@ -337,6 +340,8 @@ class TrLinear(TrCommon2): """ + _example.replace('', 'TrLinear') + """ >>> g.dist2gauss() 0.0 + >>> g2.dist2gauss() + 3.8594770921678001e-31 See also -------- @@ -399,7 +404,9 @@ class TrOchi(TrCommon2): ------- """ + _example.replace('', 'TrOchi') + """ >>> g.dist2gauss() - 5.9322684525265501 + 1.410698801056657 + >>> g2.dist2gauss() + 1.988807188766706 See also -------- @@ -486,7 +493,8 @@ class TrOchi(TrCommon2): mean = self.mean sigma = self.sigma xn = atleast_1d(x) - xn = (xn - mean) / sigma + shape0 = xn.shape + xn = (xn.ravel() - mean) / sigma igp, = where(0 <= xn) igm, = where(xn < 0) @@ -498,6 +506,7 @@ class TrOchi(TrCommon2): if gb != 0: np.put(g, igm, (-expm1(-gb * xn[igm])) / gb) + g.shape = shape0 return (g - mean2) * self.ysigma / sigma2 + self.ymean def _gauss2dat(self, y, *yi): @@ -509,7 +518,7 @@ class TrOchi(TrCommon2): sigma = self.sigma yn = (atleast_1d(y) - self.ymean) / self.ysigma - xn = sigma2 * yn + mean2 + xn = sigma2 * yn.ravel() + mean2 igp, = where(0 <= xn) igm, = where(xn < 0) @@ -519,7 +528,8 @@ class TrOchi(TrCommon2): if gb != 0: np.put(xn, igm, -log1p(-gb * xn[igm]) / gb) - + + xn.shape = yn.shape return sigma * xn + mean def main(): diff --git a/pywafo/src/wafo/wafodata.py b/pywafo/src/wafo/wafodata.py index 443d2e0..b456d41 100644 --- a/pywafo/src/wafo/wafodata.py +++ b/pywafo/src/wafo/wafodata.py @@ -2,6 +2,7 @@ import warnings from plotbackend import plotbackend from time import gmtime, strftime import numpy as np +from scipy.integrate.quadrature import cumtrapz #@UnresolvedImport __all__ = ['WafoData', 'AxisLabels'] @@ -66,19 +67,21 @@ class WafoData(object): self.data = data self.args = args self.date = now() - self.plotter = None + self.plotter = kwds.pop('plotter', None) self.children = None - self.plot_args_children = kwds.get('plot_args_children',[]) - self.plot_kwds_children = kwds.get('plot_kwds_children',{}) - self.plot_args = kwds.get('plot_args',[]) - self.plot_kwds = kwds.get('plot_kwds',{}) + self.plot_args_children = kwds.pop('plot_args_children',[]) + self.plot_kwds_children = kwds.pop('plot_kwds_children',{}) + self.plot_args = kwds.pop('plot_args',[]) + self.plot_kwds = kwds.pop('plot_kwds',{}) self.labels = AxisLabels(**kwds) - self.setplotter() + if not self.plotter: + self.setplotter() def plot(self, *args, **kwds): tmp = None - if self.children != None: + plotflag = kwds.get('plotflag', None) + if not plotflag and self.children != None: plotbackend.hold('on') tmp = [] child_args = args if len(args) else tuple(self.plot_args_children) @@ -106,12 +109,10 @@ class WafoData(object): newcopy.__dict__.update(self.__dict__) return newcopy - def setplotter(self, plotmethod=None): ''' Set plotter based on the data type data_1d, data_2d, data_3d or data_nd ''' - if isinstance(self.args, (list, tuple)): # Multidimensional data ndim = len(self.args) if ndim < 2: @@ -140,10 +141,6 @@ class AxisLabels: newcopy = empty_copy(self) newcopy.__dict__.update(self.__dict__) return newcopy - #lbkwds = self.labels.__dict__.copy() - #labels = AxisLabels(**lbkwds) - #return labels - def labelfig(self): try: h1 = plotbackend.title(self.title) @@ -178,7 +175,6 @@ class Plotter_1d(object): plotmethod = 'plot' self.plotbackend = plotbackend try: - #self.plotfun = plotbackend.__dict__[plotmethod] self.plotfun = getattr(plotbackend, plotmethod) except: pass @@ -190,9 +186,136 @@ class Plotter_1d(object): args1 = tuple((wdata.args)) + (wdata.data,) + args else: 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() 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): """ @@ -201,6 +324,7 @@ class Plotter_2d(Plotter_1d): plotmethod : string defining type of plot. Options are: contour (default) + contourf mesh surf """ @@ -209,7 +333,9 @@ class Plotter_2d(Plotter_1d): if plotmethod is None: plotmethod = 'contour' super(Plotter_2d, self).__init__(plotmethod) - #self.plotfun = plotbackend.__dict__[plotmethod] + + def _plot(self,*args,**kwds): + pass def main():