From acab68c0a1f1b788ea7e1ad031dc5ea63041a1a1 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Thu, 1 Apr 2010 09:07:13 +0000 Subject: [PATCH] Small enhancements -Added trdata method to TimeSeries class -Added TrLinear class to transform/models.py -Added more doctest examples to spectrum/core.py --- pywafo/src/wafo/misc.py | 43 +++--- pywafo/src/wafo/objects.py | 219 ++++++++++++++++++++++++-- pywafo/src/wafo/spectrum/core.py | 149 ++++++++++-------- pywafo/src/wafo/spectrum/models.py | 6 +- pywafo/src/wafo/stats/core.py | 4 +- pywafo/src/wafo/stats/estimation.py | 4 +- pywafo/src/wafo/transform/core.py | 5 +- pywafo/src/wafo/transform/models.py | 228 ++++++++++++++++------------ 8 files changed, 464 insertions(+), 194 deletions(-) diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index 20c9eb3..c4d230e 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -596,15 +596,15 @@ def rfcfilter(x, h, method=0): if cmpfun1(yi, fmi): z1 = -1 elif cmpfun1(fpi, yi): - z1 = + 1 + z1 = +1 else: z1 = 0 t1, y1 = (t0, y0) if z1 == 0 else (ti, yi) else: - if (((z0 == + 1) & cmpfun1(yi, fmi)) | ((z0 == -1) & cmpfun2(yi, fpi))): + if (((z0 == +1) & cmpfun1(yi, fmi)) | ((z0 == -1) & cmpfun2(yi, fpi))): z1 = -1 - elif (((z0 == + 1) & cmpfun2(fmi, yi)) | ((z0 == -1) & cmpfun1(fpi, yi))): - z1 = + 1 + elif (((z0 == +1) & cmpfun2(fmi, yi)) | ((z0 == -1) & cmpfun1(fpi, yi))): + z1 = +1 else: warnings.warn('Something wrong, i=%d' % tim1) @@ -614,7 +614,7 @@ def rfcfilter(x, h, method=0): elif z1 == -1: #% y1 = min([y0 xi]) t1, y1 = (t0, y0) if y0 < yi else (ti, yi) - elif z1 == + 1: + elif z1 == +1: #% y1 = max([y0 xi]) t1, y1 = (t0, y0) if y0 > yi else (ti, yi) @@ -672,6 +672,13 @@ def findtp(x, h=0.0, kind=None): >>> tph = x1[itph,:] >>> a = pylab.plot(x1[:,0],x1[:,1],tp[:,0],tp[:,1],'ro',tph[:,1],tph[:,1],'k.') >>> pylab.close('all') + >>> itp + array([ 11, 21, 22, 24, 26, 28, 31, 39, 43, 45, 47, 51, 56, + 64, 70, 78, 82, 84, 89, 94, 101, 108, 119, 131, 141, 148, + 149, 150, 159, 173, 184, 190, 199]) + >>> itph + array([ 11, 64, 28, 31, 47, 51, 39, 56, 70, 94, 78, 89, 101, + 108, 119, 148, 131, 141, 0, 159, 173, 184, 190]) See also --------- @@ -1031,7 +1038,7 @@ def common_shape(*args, ** kwds): shape = kwds.get('shape') if shape is not None: if not isinstance(shape, (list, tuple)): - shape = (shape, ) + shape = (shape,) shapes.append(tuple(shape)) if len(set(shapes)) == 1: # Common case where nothing needs to be broadcasted. @@ -1055,7 +1062,7 @@ def common_shape(*args, ** kwds): if len(unique) > 2: # There must be at least two non-1 lengths for this axis. raise ValueError("shape mismatch: two or more arrays have " - "incompatible dimensions on axis %r." % (axis, )) + "incompatible dimensions on axis %r." % (axis,)) elif len(unique) == 2: # There is exactly one non-1 length. The common shape will take this # value. @@ -1109,7 +1116,7 @@ def argsreduce(condition, * args): """ newargs = atleast_1d(*args) if not isinstance(newargs, list): - newargs = [newargs,] + newargs = [newargs, ] expand_arr = (condition == condition) return [extract(condition, arr1 * expand_arr) for arr1 in newargs] @@ -1541,15 +1548,15 @@ def meshgrid(*xi, ** kwargs): ndim = len(args) - s0 = (1, ) * ndim - output = [x.reshape(s0[:i] + (-1, ) + s0[i + 1::]) for i, x in enumerate(args)] + s0 = (1,) * ndim + output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::]) for i, x in enumerate(args)] shape = [x.size for x in output] if indexing == 'xy': # switch first and second axis - output[0].shape = (1, -1) + (1, ) * (ndim - 2) - output[1].shape = (-1, 1) + (1, ) * (ndim - 2) + output[0].shape = (1, -1) + (1,) * (ndim - 2) + output[1].shape = (-1, 1) + (1,) * (ndim - 2) shape[0], shape[1] = shape[1], shape[0] if sparse: @@ -1722,7 +1729,7 @@ def tranproc(x, f, x0, *xi): xo, fo, x0 = atleast_1d(x, f, x0) xi = atleast_1d(*xi) if not isinstance(xi, list): - xi = [xi,] + xi = [xi, ] N = len(xi) # N = number of derivatives nmax = ceil((xo.ptp()) * 10 ** (7. / max(N, 1))) xo, fo = trangood(xo, fo, min_x=min(x0), max_x=max(x0), max_n=nmax) @@ -1748,15 +1755,15 @@ def tranproc(x, f, x0, *xi): #% Transform X with the derivatives of f. fxder = zeros((N, x0.size)) - fder = vstack((xo, fo)).T + fder = vstack((xo, fo)) for k in range(N): #% Derivation of f(x) using a difference method. - n = fder.shape[0] + n = fder.shape[-1] #%fder = [(fder(1:n-1,1)+fder(2:n,1))/2 diff(fder(:,2))./diff(fder(:,1))] - fder = vstack([(fder[0:n - 1, 0] + fder[1:n, 0]) / 2, diff(fder[:, 1]) / hn]) + fder = vstack([(fder[0, 0:n - 1] + fder[0, 1:n]) / 2, diff(fder[1, :]) / hn]) fxder[k] = tranproc(fder[0], fder[1], x0) - #% Calculate the transforms of the derivatives of X. - #% First time derivative of y: y1 = f'(x)*x1 + # Calculate the transforms of the derivatives of X. + # First time derivative of y: y1 = f'(x)*x1 y1 = fxder[0] * xi[0] y.append(y1) diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index 6cac986..94fbe67 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -14,14 +14,15 @@ from __future__ import division from wafo.transform.core import TrData +from wafo.transform.models import TrHermite, TrOchi, TrLinear +from wafo.interpolate import SmoothSpline from scipy.interpolate.interpolate import interp1d from scipy.integrate.quadrature import cumtrapz -from wafo.interpolate import SmoothSpline import warnings import numpy as np from numpy import (inf, pi, zeros, ones, sqrt, where, log, exp, sin, arcsin, mod, finfo, interp, #@UnresolvedImport - newaxis, linspace, arange, sort, all, abs, linspace, vstack, hstack, atleast_1d, #@UnresolvedImport + newaxis, linspace, arange, sort, all, abs, vstack, hstack, atleast_1d, #@UnresolvedImport polyfit, r_, nonzero, cumsum, ravel, size, isnan, nan, floor, ceil, diff, array) #@UnresolvedImport from numpy.fft import fft from numpy.random import randn @@ -359,13 +360,13 @@ class LevelCrossings(WafoData): if ng == 1: gvar = opt.gvar * ones(ncr) 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)) + cvar = interp1d(linspace(0, 1, ng), opt.cvar, kind='linear')(linspace(0, 1, ncr)) @@ -383,7 +384,7 @@ class LevelCrossings(WafoData): cor2 = 0 - lc22 = hstack((0,cumtrapz(lc2, lc1) + cor1)) + lc22 = hstack((0, cumtrapz(lc2, lc1) + cor1)) lc22 = (lc22 + 0.5) / (lc22[-1] + cor2 + 1) lc11 = (lc1 - mean) / sigma @@ -401,9 +402,9 @@ class LevelCrossings(WafoData): #u0 = interp1q(cros(:,2),cros(:,1),.5) - lc22 = ndtri(lc22)-u0 #invnorm(lc22, -u0, 1); + lc22 = ndtri(lc22) - u0 #invnorm(lc22, -u0, 1); - g2 = TrData(lc22.copy(), lc1.copy(), mean, sigma**2) + g2 = TrData(lc22.copy(), lc1.copy(), mean, sigma ** 2) # NB! the smooth function does not always extrapolate well outside the edges # causing poor estimate of g # We may alleviate this problem by: forcing the extrapolation @@ -412,7 +413,7 @@ class LevelCrossings(WafoData): inds = slice(Ne, ncr - Ne) # indices to points we are smoothing over scros2 = SmoothSpline(lc11[inds], lc22[inds], opt.gsm, opt.lin_extrap, gvar[inds])(uu) - g = TrData(scros2, g1, mean, sigma**2) #*sa; #multiply with stdev + g = TrData(scros2, g1, mean, sigma ** 2) #*sa; #multiply with stdev if opt.chkder: for ix in range(5): @@ -603,8 +604,6 @@ class TurningPoints(WafoData): ---------------- data : array_like args : vector for 1D - - ''' def __init__(self, *args, **kwds): super(TurningPoints, self).__init__(*args, **kwds) @@ -857,7 +856,207 @@ class TimeSeries(WafoData): fact = 2.0 * pi w = fact * f return _wafospec.SpecData1D(S / fact, w) + def trdata(self, method='nonlinear', **options): + ''' + Estimate transformation, g, from data. + + CALL: [g test cmax irr g2] = dat2tr(x,def,options); + + g,g2 = the smoothed and empirical transformation, respectively. + A two column matrix if multip=0. + If multip=1 it is a 2*(m-1) column matrix where the + first and second column is the transform + for values in column 2 and third and fourth column is the + transform for values in column 3 ...... + + test = int (g(u)-u)^2 du where int. limits is given by param. This + is a measure of departure of the data from the Gaussian model. + + Parameters + ---------- + + method : string + 'nonlinear' : transform based on smoothed crossing intensity (default) + 'mnonlinear': transform based on smoothed marginal distribution + 'hermite' : transform based on cubic Hermite polynomial + 'ochi' : transform based on exponential function + 'linear' : identity. + + options = options structure with the following fields: + csm,gsm - defines the smoothing of the logarithm of crossing intensity + and the transformation g, respectively. Valid values must + be 0<=csm,gsm<=1. (default csm=0.9, gsm=0.05) + Smaller values gives smoother functions. + param - vector which defines the region of variation of the data x. + (default see lc2tr). + plotflag - 0 no plotting (Default) + 1 plots empirical and smoothed g(u) and the theoretical for + a Gaussian model. + 2 monitor the development of the estimation + linextrap - 0 use a regular smoothing spline + 1 use a smoothing spline with a constraint on the ends to + ensure linear extrapolation outside the range of the data. + (default) + gvar - Variances for the empirical transformation, g. (default 1) + ne - Number of extremes (maxima & minima) to remove from the + estimation of the transformation. This makes the + estimation more robust against outliers. (default 7) + ntr - Maximum length of empirical crossing intensity or CDF. + The empirical crossing intensity or CDF is interpolated + linearly before smoothing if their lengths exceeds Ntr. + A reasonable NTR will significantly speed up the + estimation for long time series without loosing any + accuracy. NTR should be chosen greater than + PARAM(3). (default 1000) + multip - 0 the data in columns belong to the same seastate (default). + 1 the data in columns are from separate seastates. + + DAT2TR estimates the transformation in a transformed Gaussian model. + Assumption: a Gaussian process, Y, is related to the + non-Gaussian process, X, by Y = g(X). + + The empirical crossing intensity is usually very irregular. + More than one local maximum of the empirical crossing intensity + may cause poor fit of the transformation. In such case one + should use a smaller value of CSM. In order to check the effect + of smoothing it is recomended to also plot g and g2 in the same plot or + plot the smoothed g against an interpolated version of g (when CSM=GSM=1). + If x is likely to cross levels higher than 5 standard deviations + then the vector param has to be modified. For example if x is + unlikely to cross a level of 7 standard deviations one can use + PARAM=[-7 7 513]. + + Example + ------- + >>> import wafo.spectrum.models as sm + >>> import wafo.transform.models as tm + >>> from wafo.objects import mat2timeseries + >>> Hs = 7.0 + >>> Sj = sm.Jonswap(Hm0=Hs) + >>> S = Sj.tospecdata() #Make spectrum object from numerical values + >>> S.tr = tm.TrOchi(mean=0, skew=0.16, kurt=0, sigma=Hs/4, ysigma=Hs/4) + >>> xs = S.sim(ns=2**16) + >>> ts = mat2timeseries(xs) + >>> g0, gemp = ts.trdata(monitor=True) # Monitor the development + >>> g1, gemp = ts.trdata(method='m', gvar=0.5 ) # Equal weight on all points + >>> g2, gemp = ts.trdata(method='n', gvar=[3.5, 0.5, 3.5]) # Less weight on the ends + >>> S.tr.dist2gauss() + 5.9322684525265501 + >>> np.round(gemp.dist2gauss()) + 6.0 + >>> np.round(g0.dist2gauss()) + 4.0 + >>> np.round(g1.dist2gauss()) + 4.0 + >>> np.round(g2.dist2gauss()) + 4.0 + + 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 + -------- + troptset, lc2tr, cdf2tr, trplot + References + ---------- + Rychlik, I. , Johannesson, P and Leadbetter, M. R. (1997) + "Modelling and statistical analysis of ocean wavedata using + transformed Gaussian process." + Marine structures, Design, Construction and Safety, Vol. 10, No. 1, pp 13--47 + + + Brodtkorb, P, Myrhaug, D, and Rue, H (1999) + "Joint distribution of wave height and crest velocity from + reconstructed data" + in Proceedings of 9th ISOPE Conference, Vol III, pp 66-73 + ''' +# Tested on: Matlab 5.3, 5.2, 5.1 +# History: +# revised pab Dec2004 +# -Fixed bug: string comparison for def at fault. +# revised pab Nov2004 +# -Fixed bug: linextrap was not accounted for +# revised pab july 2004 +# revised pab 3 april 2004 +# -fixed a bug in hermite estimation: excess changed to kurtosis +# revised pab 29.12.2000 +# - added example, hermite and ochi options +# - replaced optional arguments with a options struct +# - default param is now [-5 5 513] -> better to have the discretization +# represented with exact numbers, especially when calculating +# derivatives of the transformation numerically. +# revised pab 19.12.2000 +# - updated call edf(X,-inf,[],monitor) to edf(X,[],monitor) +# due to new calling syntax for edf +# modifed pab 24.09.2000 +# - changed call from norminv to wnorminv +# - also removed the 7 lowest and 7 highest points from +# the estimation using def='mnonlinear' +# (This is similar to what lc2tr does. lc2tr removes +# the 9 highest and 9 lowest TP from the estimation) +# modified pab 09.06.2000 +# - made all the *empirical options secret. +# - Added 'mnonlinear' and 'mempirical' +# - Fixed the problem of multip==1 and def=='empirical' by interpolating +# with spline to ensure that the length of g is fixed +# - Replaced the test statistic for def=='empirical' with the one +# obtained when csm1=csm2=1. (Previously only the smoothed test +# statistic where returned) +# modified pab 12.10.1999 +# fixed a bug +# added secret output of empirical estimate g2 +# modified by svi 29.09.1999 +# changed input def by adding new options. +# revised by pab 11.08.99 +# changed name from dat2tran to dat2tr +# modified by Per A. Brodtkorb 12.05.1999,15.08.98 +# added secret option: to accept multiple data, to monitor the steps +# of estimation of the transformation +# also removed some code and replaced it with a call to lc2tr (cross2tr) +# making the maintainance easier +# + + #opt = troptset('plotflag','off','csm',.95,'gsm',.05,.... + # 'param',[-5 5 513],'delay',2,'linextrap','on','ne',7,... + # 'cvar',1,'gvar',1,'multip',0); + + + opt = DotDict(chkder=True, plotflag=True, csm=.95, gsm=.05, + param=[-5, 5, 513], delay=2, ntr=inf, linextrap=True, ne=7, cvar=1, gvar=1, + multip=False, crossdef='uM') + opt.update(**options) + + ma = self.data.mean() + sa = self.data.std() + if method.startswith('lin'): + return TrLinear(mean=ma, sigma=sa) + + if method[0] == 'n': + tp = self.turning_points() + mM = tp.cycle_pairs() + lc = mM.level_crossings(opt.crossdef) + return lc.trdata() + elif method[0] == 'm': + return cdftr() + elif method[0] == 'h': + ga1 = np.skew(self.data) + ga2 = np.kurtosis(self.data, fisher=True) #kurt(xx(n+1:end))-3; + up = min(4 * (4 * ga1 / 3) ** 2, 13) + lo = (ga1 ** 2) * 3 / 2; + kurt1 = min(up, max(ga2, lo)) + 3 + return TrHermite(mean=ma, var=sa ** 2, skew=ga1, kurt=kurt1) + elif method[0] == 'o': + ga1 = np.skew(self.data) + return TrOchi(mean=ma, var=sa ** 2, skew=ga1) + def turning_points(self, h=0.0, wavetype=None): ''' Return turning points (tp) from data, optionally rainflowfiltered. diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py index 583d5c2..9c88a1b 100644 --- a/pywafo/src/wafo/spectrum/core.py +++ b/pywafo/src/wafo/spectrum/core.py @@ -1,5 +1,6 @@ from __future__ import division from scipy.misc.ppimport import ppimport +from wafo.objects import mat2timeseries, TimeSeries import warnings import numpy as np from numpy import (pi, inf, meshgrid, zeros, ones, where, nonzero, #@UnresolvedImport @@ -17,7 +18,8 @@ from pylab import stineman_interp from dispersion_relation import w2k #, k2w from wafo.wafodata import WafoData, now -from wafo.misc import sub_dict_select, nextpow2, discretize, JITImport, tranproc + +from wafo.misc import sub_dict_select, nextpow2, discretize, JITImport #, tranproc try: from wafo.gaussian import Rind except ImportError: @@ -28,7 +30,8 @@ except ImportError: warnings.warn('Compile the c_library.pyd again!') c_library = None -from wafo.transform import TrData +#from wafo.transform import TrData +from wafo.transform.models import TrLinear from wafo.plotbackend import plotbackend @@ -190,7 +193,7 @@ class SpecData1D(WafoData): self.freqtype = 'w' self.angletype = '' self.h = inf - self.tr = None + self.tr = None #TrLinear() self.phi = 0.0 self.v = 0.0 self.norm = False @@ -572,7 +575,9 @@ class SpecData1D(WafoData): #Fs = 2*freq(end)+eps; % sampling frequency for ix in xrange(max_sim): - [x2, x1] = spec2nlsdat(SL, [np, cases], [], iseed, method, fnLimit) + #[x2, x1] = spec2nlsdat(SL, [np, cases], [], iseed, method, fnLimit) + [x2, x1] = self.sim_nl(ns=np, cases=cases, dt=None, iseed=iseed, method=method, + fnlimit=fn_limit) #%x2(:,2:end) = x2(:,2:end) -x1(:,2:end); S2 = dat2spec(x2, L) S1 = dat2spec(x1, L) @@ -727,9 +732,10 @@ class SpecData1D(WafoData): if self.tr is None: - y = linspace(-5, 5, 513) + g = TrLinear(var=m[0]) + #y = linspace(-5, 5, 513) #g = _wafotransform. - g = TrData(y, sqrt(m[0]) * y) + #g = TrData(y, sqrt(m[0]) * y) else: g = self.tr @@ -1053,7 +1059,7 @@ class SpecData1D(WafoData): xder[:, 0] = x[:, 0] if spec.tr is not None: - print(' Transforming data.') + #print(' Transforming data.') g = spec.tr if derivative: for i in range(cases): @@ -1079,7 +1085,7 @@ class SpecData1D(WafoData): # function [x2,x,svec,dvec,amp]=spec2nlsdat(spec,np,dt,iseed,method,truncationLimit) def sim_nl(self, ns=None, cases=1, dt=None, iseed=None, method='random', - fnlimit=1.4142, reltol=1e-3, g=9.81): + fnlimit=1.4142, reltol=1e-3, g=9.81, verbose=False): """ Simulates a Randomized 2nd order non-linear wave X(t) @@ -1142,6 +1148,21 @@ class SpecData1D(WafoData): Example -------- + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap();S = Sj.tospecdata() + >>> ns =100; dt = .2 + >>> x1 = S.sim_nl(ns,dt=dt) + + >>> import numpy as np + >>> import scipy.stats as st + >>> x2 = S.sim_nl(ns=20000,cases=20) + >>> truth1 = [0,np.sqrt(S.moment(1)[0])] + S.stats_nl(moments='sk') + >>> funs = [np.mean,np.std,st.skew,st.kurtosis] + >>> for fun,trueval in zip(funs,truth1): + ... res = fun(x2[:,1::], axis=0) + ... m = res.mean() + ... sa = res.std() + ... assert(np.abs(m-trueval)>> import wafo.spectrum.models as sm - >>> Sj = sm.Jonswap() + >>> import wafo.transform.models as wtm + >>> Hs = 7. + >>> Sj = sm.Jonswap(Hm0=Hs, Tp=11) >>> S = Sj.tospecdata() >>> me, va, sk, ku = S.stats_nl(moments='mvsk') - - - Hm0=7;Tp=11 - S = jonswap([],[Hm0 Tp]); [sk, ku, me]=spec2skew(S) - g=hermitetr([],[Hm0/4 sk ku me]); g2=[g(:,1), g(:,2)*Hm0/4] - ys = spec2sdat(S,15000) % Simulated in the Gaussian world - xs = gaus2dat(ys,g2) % Transformed to the real world + >>> g = wtm.TrHermite(mean=me, sigma=Hs/4, skew=sk, kurt=ku, ysigma=Hs/4) + >>> ys = S.sim(15000) # Simulated in the Gaussian world + >>> xs = g.gauss2dat(ys[:,1]) # Transformed to the real world + See also --------- @@ -1483,7 +1503,7 @@ class SpecData1D(WafoData): ## skew = sum((6*C2+8*E2).*E)/sa^3 % skewness ## kurt = 3+48*sum((C2+E2).*E2)/sa^4 % kurtosis return output - def testgaussian(ns,test0=None, cases=100, method='nonlinear',**opt): + def testgaussian(self, ns,test0=None, cases=100, method='nonlinear',**opt): ''' TESTGAUSSIAN Test if a stochastic process is Gaussian. @@ -1507,13 +1527,24 @@ class SpecData1D(WafoData): 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. - Example: - Hm0 = 7; - S0 = jonswap([],Hm0); g=ochitr([],[Hm0/4]); S=S0; - S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4; - xs = spec2sdat(S,2^13); - [g0 t0] = dat2tr(xs); - t1 = testgaussian(S0,[2^13 50],t0); + Example: + ------- + >>> import wafo.spectrum.models as sm + >>> import wafo.transform.models as wtm + >>> import wafo.objects as wo + >>> Hs = 7 + >>> Sj = sm.Jonswap(Hm0=Hs) + >>> S0 = Sj.tospecdata() + >>> ns =100; dt = .2 + >>> x1 = S0.sim(ns, dt=dt) + + >>> S = S0.copy() + >>> me, va, sk, ku = S.stats_nl(moments='mvsk') + >>> S.tr = wtm.TrHermite(mean=me, sigma=Hs/4, skew=sk, kurt=ku, ysigma=Hs/4) + >>> ys = wo.mat2timeseries(S.sim(ns=2**13)) + >>> g0, gemp = ys.trdata() + >>> t0 = g0.dist2gauss() + >>> t1 = S0.testgaussian(ns=2**13, t0=t0, cases=50) See also cov2sdat, dat2tr, troptset ''' @@ -1542,49 +1573,43 @@ class SpecData1D(WafoData): # # opt = troptset(opt,'multip',1) - if test0 is None: - plotflag=0 - else: - plotflag=1 - + plotflag=0 if test0 is None else 1 if cases>50: print(' ... be patient this may take a while') - test1 = [] - rep = floor(ns*cases/maxsize)+1 + + rep = int(np.floor(ns*cases/maxsize)+1) - Nstep = floor(cases/rep); + Nstep = np.floor(cases/rep); acf = self.tocovdata() #R = spec2cov(S); - + test1 = [] for ix in range(rep): xs = acf.sim(ns=ns, cases=Nstep) + for iy in range(1, xs.shape[-1]): + ts = TimeSeries(xs[:, iy], xs[:, 0].ravel()) + g, tmp = ts.trdata(method, **opt) + test1.append(g.dist2gauss()) #xs = cov2sdat(R,[ns Nstep]); - [g, tmp] = dat2tr(xs,method, **opt); + #[g, tmp] = dat2tr(xs,method, **opt); #test1 = [test1; tmp(:)] - print('finished %d of %d ' % (ix,rep) ) + print('finished %d of %d ' % (ix+1,rep) ) if rep>1: - xs = acf.sim(ns=ns, cases=rem(cases,rep)) - [g, tmp] = dat2tr(xs,method,**opt); - #test1 = [test1; tmp(:)]; - - -# if (nargout>0 || plotflag==0), -# test2=test1; -# end -# -# -# if plotflag -# plot(test1,'o'),hold on -# if 1 -# plot([1 cases],test0*[1 1],'--'), -# end -# hold off -# ylabel('e(g(u)-u)') -# xlabel('Simulation number') -# end + xs = acf.sim(ns=ns, cases=np.remainder(cases,rep)) + for iy in range(1, xs.shape[-1]): + ts = TimeSeries(xs[:, iy], xs[:, 0].ravel()) + g, tmp = ts.trdata(method, **opt) + test1.append(g.dist2gauss()) + + if plotflag: + plotbackend.plot(test1,'o') + plotbackend.plot([1, cases], [test0, test0],'--') + + plotbackend.ylabel('e(g(u)-u)') + plotbackend.xlabel('Simulation number') + return test1 def moment(self, nr=2, even=True, j=0): ''' Calculates spectral moments from spectrum @@ -1647,8 +1672,8 @@ class SpecData1D(WafoData): if self.freqtype in ['f', 'w']: vari = 't' if self.freqtype == 'f': - f = 2. * pi * f - S = S / (2. * pi) + f = 2. * pi * f + S = S / (2. * pi) else: vari = 'x' S1 = abs(S) ** (j + 1.) @@ -2326,7 +2351,7 @@ class SpecData2D(WafoData): ##% By es 27.08.1999 - pi = pi + two_dim_spectra = ['dir', 'encdir', 'k2d'] if self.type not in two_dim_spectra: raise ValueError('Unknown 2D spectrum type!') diff --git a/pywafo/src/wafo/spectrum/models.py b/pywafo/src/wafo/spectrum/models.py index 1e02867..9f90f22 100644 --- a/pywafo/src/wafo/spectrum/models.py +++ b/pywafo/src/wafo/spectrum/models.py @@ -55,7 +55,7 @@ from numpy import (inf, atleast_1d, newaxis, any, minimum, maximum, array, #@Unr from dispersion_relation import w2k #ppimport.enable() #_wafospectrum = ppimport.ppimport('wafo.spectrum') -from core import SpecData1D +from wafo.spectrum import SpecData1D sech = lambda x: 1.0 / cosh(x) eps = finfo(float).eps @@ -542,9 +542,9 @@ class Jonswap(ModelSpectrum): outsideJonswapRange = Tp > 5 * sqrt(Hm0) or Tp < 3.6 * sqrt(Hm0) if outsideJonswapRange: txt0 = ''' - Hm0,Tp is outside the JONSWAP range. + Hm0=%g,Tp=%g is outside the JONSWAP range. The validity of the spectral density is questionable. - ''' + ''' % (Hm0, Tp) warnings.warn(txt0) if gam < 1 or 7 < gam: diff --git a/pywafo/src/wafo/stats/core.py b/pywafo/src/wafo/stats/core.py index 8380c43..cf6c43a 100644 --- a/pywafo/src/wafo/stats/core.py +++ b/pywafo/src/wafo/stats/core.py @@ -8,7 +8,7 @@ __all__ = ['edf'] def edf(x, method=2): ''' - Returns EDF Empirical Distribution Function. + Returns Empirical Distribution Function (EDF). Parameters ---------- @@ -48,7 +48,7 @@ def edf(x, method=2): def edfcnd(x, c=None, method=2): ''' - EDFCND Empirical Distribution Function CoNDitioned that X>=c. + Returns empirical Distribution Function CoNDitioned that X>=c (EDFCND). Parameters ---------- diff --git a/pywafo/src/wafo/stats/estimation.py b/pywafo/src/wafo/stats/estimation.py index b75db37..5bf0d2c 100644 --- a/pywafo/src/wafo/stats/estimation.py +++ b/pywafo/src/wafo/stats/estimation.py @@ -40,7 +40,7 @@ arr = asarray all = alltrue -def chi2isf(p,df): +def chi2isf(p, df): return special.chdtri(df, p) def chi2sf(x, df): return special.chdtrc(df, x) @@ -92,7 +92,7 @@ class rv_frozen(object): def __init__(self, dist, *args, **kwds): self.dist = dist loc0, scale0 = map(kwds.get, ['loc', 'scale']) - if hasattr(dist, 'fix_loc_scale'): #isinstance(dist, _WAFODIST.rv_continuous): + if hasattr(dist, 'fix_loc_scale'): #isinstance(dist, rv_continuous): args, loc0, scale0 = dist.fix_loc_scale(args, loc0, scale0) self.par = args + (loc0, scale0) else: # rv_discrete diff --git a/pywafo/src/wafo/transform/core.py b/pywafo/src/wafo/transform/core.py index 39b7a1d..8ab3de5 100644 --- a/pywafo/src/wafo/transform/core.py +++ b/pywafo/src/wafo/transform/core.py @@ -48,8 +48,8 @@ class TrCommon(object): self.ymean = kwds.get('ymean', 0e0) self.ysigma = kwds.get('ysigma', 1e0) - def __call__(self, x): - return self._dat2gauss(x) + def __call__(self, x, *xi): + return self._dat2gauss(x, *xi) def dist2gauss(self, x=None, xnmin=-5, xnmax=5, n=513): """ @@ -158,6 +158,7 @@ class TrData(WafoData, TrCommon): 1 >>> g.sigma 5 + >>> g.dat2gauss(1,2,3) Check that the departure from a Gaussian model is zero diff --git a/pywafo/src/wafo/transform/models.py b/pywafo/src/wafo/transform/models.py index c870e2a..c81c9db 100644 --- a/pywafo/src/wafo/transform/models.py +++ b/pywafo/src/wafo/transform/models.py @@ -15,20 +15,18 @@ TrOchi # Licence: #------------------------------------------------------------------------------- #!/usr/bin/env python - +from __future__ import division from scipy.optimize import brentq from numpy import (sqrt, atleast_1d, abs, imag, sign, where, cos, arccos, ceil, #@UnresolvedImport expm1, log1p, pi) #@UnresolvedImport import numpy as np import warnings from core import TrCommon -__all__=['TrHermite','TrOchi'] +__all__ = ['TrHermite', 'TrLinear', 'TrOchi'] _example = ''' >>> std = 7./4 >>> g = (sigma=std, ysigma=std) - >>> g.dist2gauss() - 3.9858776379926808 Simulate a Transformed Gaussian process: >>> import numpy as np @@ -43,8 +41,9 @@ _example = ''' >>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world ''' + class TrHermite(TrCommon): - __doc__ = TrCommon.__doc__.replace('','Hermite') + """ + __doc__ = TrCommon.__doc__.replace('', 'Hermite') + """ pardef : scalar, integer 1 Winterstein et. al. (1994) parametrization [1]_ (default) 2 Winterstein (1988) parametrization [2]_ @@ -79,7 +78,9 @@ class TrHermite(TrCommon): Example: -------- - """ + _example.replace('','TrHermite') + """ + """ + _example.replace('', 'TrHermite') + """ + >>> g.dist2gauss() + 3.9858776379926808 See also -------- @@ -96,9 +97,9 @@ class TrHermite(TrCommon): 'Nonlinear vibration models for extremes and fatigue.' J. Engng. Mech., ASCE, Vol 114, No 10, pp 1772-1790 """ - def __init__(self,*args, **kwds): + def __init__(self, *args, **kwds): super(TrHermite, self).__init__(*args, **kwds) - self.pardef = kwds.get('pardef',1) + self.pardef = kwds.get('pardef', 1) self._c3 = None self._c4 = None self._forward = None @@ -108,32 +109,32 @@ class TrHermite(TrCommon): def _poly_par_from_stats(self): skew = self.skew - ga2 = self.kurt-3.0 + ga2 = self.kurt - 3.0 if ga2 <= 0: - self._c4 = ga2/24. - self._c3 = skew/6. + self._c4 = ga2 / 24. + self._c3 = skew / 6. elif self.pardef == 2: # Winterstein 1988 parametrization - if skew**2>8*(ga2+3.)/9.: + if skew ** 2 > 8 * (ga2 + 3.) / 9.: warnings.warn('Kurtosis too low compared to the skewness') - self._c4 = (sqrt(1.+1.5*ga2)-1.)/18. - self._c3 = skew/(6.*(1+6.*self._c4)) + self._c4 = (sqrt(1. + 1.5 * ga2) - 1.) / 18. + self._c3 = skew / (6. * (1 + 6. * self._c4)) else: # Winterstein et. al. 1994 parametrization intended to # apply for the range: 0 <= ga2 < 12 and 0<= skew^2 < 2*ga2/3 - if skew**2 > 2*(ga2)/3: + if skew ** 2 > 2 * (ga2) / 3: warnings.warn('Kurtosis too low compared to the skewness') if (ga2 < 0) or (12 < ga2): warnings.warn('Kurtosis must be between 0 and 12') - self._c3 = skew/6*(1-0.015*abs(skew)+0.3*skew**2)/(1+0.2*ga2) + self._c3 = skew / 6 * (1 - 0.015 * abs(skew) + 0.3 * skew ** 2) / (1 + 0.2 * ga2) if ga2 == 0.: self._c4 = 0.0 else: - c41 = (1.-1.43*skew**2./ga2)**(1.-0.1*(ga2+3.)**0.8) - self._c4 = 0.1*((1.+1.25*ga2)**(1./3.)-1.)*c41 + c41 = (1. - 1.43 * skew ** 2. / ga2) ** (1. - 0.1 * (ga2 + 3.) ** 0.8) + self._c4 = 0.1 * ((1. + 1.25 * ga2) ** (1. / 3.) - 1.) * c41 if not np.isfinite(self._c3) or not np.isfinite(self._c4): raise ValueError('Unable to calculate the polynomial') @@ -150,31 +151,31 @@ class TrHermite(TrCommon): c4 = self._c4 ma = self.mean sa = self.sigma - if abs(c4)0: + if r.size > 0: # Compute where it is possible to invert the polynomial - if self.kurt<3.: + if self.kurt < 3.: self._x_limit = r else: - self._x_limit = sa*p(r)+ma + self._x_limit = sa * p(r) + ma txt1 = ''' The polynomial is not a strictly increasing function. The derivative of g(x) is infinite at x = %g''' % self._x_limit @@ -183,82 +184,82 @@ class TrHermite(TrCommon): def check_forward(self, x): if not (self._x_limit is None): x00 = self._x_limit - txt2 = 'for the given interval x = [%g, %g]' % (x[0],x[-1]) + txt2 = 'for the given interval x = [%g, %g]' % (x[0], x[-1]) - if any(np.logical_and(x[0]<= x00, x00 <= x[-1])): + if any(np.logical_and(x[0] <= x00, x00 <= x[-1])): cdef = 1 else: - cdef = sum( np.logical_xor(x00 <= x[0] , x00 <= x[-1])) + cdef = sum(np.logical_xor(x00 <= x[0] , x00 <= x[-1])) - if np.mod(cdef,2): + if np.mod(cdef, 2): errtxt = 'Unable to invert the polynomial \n %s' % txt2 raise ValueError(errtxt) np.disp('However, successfully inverted the polynomial\n %s' % txt2) - def _dat2gauss(self,x, *xi): - if len(xi)>0: + def _dat2gauss(self, x, *xi): + if len(xi) > 0: raise ValueError('Transforming derivatives is not implemented!') xn = atleast_1d(x) self.check_forward(xn) - xn = (xn-self.mean)/self.sigma + xn = (xn - self.mean) / self.sigma if self._forward is None: #Inverting the polynomial - yn = self._poly_inv(self._backward,xn) + yn = self._poly_inv(self._backward, xn) else: yn = self._forward(xn) - return yn*self.ysigma + self.ymean + return yn * self.ysigma + self.ymean def _gauss2dat(self, y, *yi): - if len(yi)>0: + if len(yi) > 0: raise ValueError('Transforming derivatives is not implemented!') - yn = (atleast_1d(y)-self.mean)/self.ysigma + yn = (atleast_1d(y) - self.ymean) / self.ysigma #self.check_forward(y) if self._backward is None: #% Inverting the polynomial #%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - xn = self._poly_inv(self._forward,yn) + xn = self._poly_inv(self._forward, yn) else: xn = self._backward(yn) - return self.sigma*xn + self.mean + return self.sigma * xn + self.mean def _poly_inv(self, p, xn): ''' Invert polynomial ''' - if p.order<2: + if p.order < 2: return xn - elif p.order==2: + elif p.order == 2: # Quadratic: Solve a*u**2+b*u+c = xn coefs = p.coeffs a = coefs[0] b = coefs[1] - c = coefs[2]-xn - t = 0.5*(b+sign(b)*sqrt(b**2-4*a*c)) + c = coefs[2] - xn + t = 0.5 * (b + sign(b) * sqrt(b ** 2 - 4 * a * c)) #so1 = t/a # largest solution - so2 = -c/t # smallest solution + so2 = -c / t # smallest solution return so2 - elif p.order==3: + elif p.order == 3: # Solve # K*(c4*u^3+c3*u^2+(1-3*c4)*u-c3) = xn = (x-ma)/sa # -c4*xn^3-c3*xn^2+(1+3*c4)*xn+c3 = u - coefs = p.coeffs[1::]/p.coeffs[0] + coefs = p.coeffs[1::] / p.coeffs[0] a = coefs[0] b = coefs[1] - c = coefs[2]-xn/p.coeffs[0] + c = coefs[2] - xn / p.coeffs[0] - x0 = a/3. + x0 = a / 3. #% substitue xn = z-x0 and divide by c4 => z^3 + 3*p1*z+2*q0 = 0 - p1 = b/3-x0**2 + p1 = b / 3 - x0 ** 2 #p1 = (b-a**2/3)/3 #q0 = (c + x0*(2.*x0/3.-b))/2. #q0 = x0**3 -a*b/6 +c/2 - q0 = x0*(x0**2-b/2)+c/2 + q0 = x0 * (x0 ** 2 - b / 2) + c / 2 ## # z^3+3*p1*z+2*q0=0 ## c3 = self._c3 @@ -272,24 +273,59 @@ class TrHermite(TrCommon): #q0 = x0**3-1.5*b1*(x0+xn) if not (self._x_limit is None): # % Three real roots d = sqrt(-p1) - theta1 = arccos(-q0/d**3)/3 - th2 = np.r_[0, -2*pi/3, 2*pi/3] - x1 = abs(2*d*cos(theta1[ceil(len(xn)/2)] + th2)-x0) + theta1 = arccos(-q0 / d ** 3) / 3 + th2 = np.r_[0, -2 * pi / 3, 2 * pi / 3] + x1 = abs(2 * d * cos(theta1[ceil(len(xn) / 2)] + th2) - x0) ix = x1.argmin() # % choose the smallest solution - return 2.*d*cos(theta1 + th2[ix])-x0 + return 2. * d * cos(theta1 + th2[ix]) - x0 else: # %Only one real root exist - q1 = sqrt((q0)**2+p1**3) + q1 = sqrt((q0) ** 2 + p1 ** 3) #% Find the real root of the monic polynomial - A0 = (q1-q0)**(1./3.) - B0 = -(q1+q0)**(1./3.) - return A0+B0-x0 #% real root + A0 = (q1 - q0) ** (1. / 3.) + B0 = -(q1 + q0) ** (1. / 3.) + return A0 + B0 - x0 #% real root #%% The other complex roots are given by #%x= -(A0+B0)/2+(A0-B0)*sqrt(3)/2-x0 #%x=-(A0+B0)/2+(A0-B0)*sqrt(-3)/2-x0 + +class TrLinear(TrCommon): + __doc__ = TrCommon.__doc__.replace('', 'Linear') + """ + Description + ----------- + The linear transformation model is monotonic linear polynomial, calibrated + such that the first 2 moments of the transformed model G(y)=g^-1(y) match + the moments of the true process. + + Example: + -------- + """ + _example.replace('', 'TrLinear') + """ + >>> g.dist2gauss() + 0.0 + + See also + -------- + spec2skew, ochitr, lc2tr, dat2tr + + """ + def _dat2gauss(self, x, *xi): + sratio = atleast_1d(self.ysigma / self.sigma) + y = (atleast_1d(x) - self.mean) * sratio + self.ymean + if len(xi) > 0: + y = [y, ] + [ ix * sratio for ix in xi] + return y + + def _gauss2dat(self, y, *yi): + sratio = atleast_1d(self.sigma / self.ysigma) + x = (atleast_1d(y) - self.ymean) * sratio + self.mean + if len(yi) > 0: + x = [x, ] + [iy * sratio for iy in yi] + return x + + class TrOchi(TrCommon): - __doc__ = TrCommon.__doc__.replace('','Ochi') + """ + __doc__ = TrCommon.__doc__.replace('', 'Ochi') + """ Description ----------- @@ -322,8 +358,10 @@ class TrOchi(TrCommon): Example ------- - """ + _example.replace('','TrOchi') + """ - + """ + _example.replace('', 'TrOchi') + """ + >>> g.dist2gauss() + 5.9322684525265501 + See also -------- spec2skew, hermitetr, lc2tr, dat2tr @@ -347,13 +385,13 @@ class TrOchi(TrCommon): def _par_from_stats(self): skew = self.skew - if abs(skew)>2.82842712474619: + if abs(skew) > 2.82842712474619: raise ValueError('Skewness must be less than 2.82842') mean1 = self.mean sigma1 = self.sigma - if skew==0: + if skew == 0: self._phat = [sigma1, mean1, 0, 0, 1, 0] return @@ -367,22 +405,22 @@ class TrOchi(TrCommon): # Set up the 2D non-linear equations for a and sig2^2: # g1='[x(2)-2.*x(1).^2.*x(2).^2-P1, 2.*x(1).*x(2).^2.*(3-8.*x(1).^2.*x(2))-P2 ]' # Or solve the following 1D non-linear equation for sig2^2: - g2 = lambda x: -sqrt(abs(x-1)*2)*(3.*x-4*abs(x-1))+abs(skew) + g2 = lambda x:-sqrt(abs(x - 1) * 2) * (3. * x - 4 * abs(x - 1)) + abs(skew) a1 = 1. # % Start interval where sig2^2 is located. a2 = 2. - sig22 = brentq(g2,a1,a2) #% smallest solution for sig22 - a = sign(skew)*sqrt(abs(sig22-1)/2)/sig22 - gam_a = 1.28*a - gam_b = 3*a + sig22 = brentq(g2, a1, a2) #% smallest solution for sig22 + a = sign(skew) * sqrt(abs(sig22 - 1) / 2) / sig22 + gam_a = 1.28 * a + gam_b = 3 * a sigma2 = sqrt(sig22) #% Solve the following 2nd order equation to obtain ma2 #% a*(sig2^2+ma2^2)+ma2 = 0 - my2 = (-1.-sqrt(1.-4.*a**2*sig22))/a #% Largest mean - mean2 = a*sig22/my2 #% choose the smallest mean + my2 = (-1. - sqrt(1. - 4. * a ** 2 * sig22)) / a #% Largest mean + mean2 = a * sig22 / my2 #% choose the smallest mean self._phat = [sigma1, mean1, gam_a, gam_b, sigma2, mean2] return @@ -391,8 +429,8 @@ class TrOchi(TrCommon): ''' Returns ga, gb, sigma2, mean2 ''' - if (self._phat is None or self.sigma!=self._phat[0] - or self.mean!=self._phat[1]): + if (self._phat is None or self.sigma != self._phat[0] + or self.mean != self._phat[1]): self._par_from_stats() #sigma1 = self._phat[0] #mean1 = self._phat[1] @@ -402,63 +440,63 @@ class TrOchi(TrCommon): mean2 = self._phat[5] return ga, gb, sigma2, mean2 - def _dat2gauss(self,x, *xi): - if len(xi)>0: + def _dat2gauss(self, x, *xi): + if len(xi) > 0: raise ValueError('Transforming derivatives is not implemented!') ga, gb, sigma2, mean2 = self._get_par() mean = self.mean sigma = self.sigma xn = atleast_1d(x) - xn = (xn-mean)/sigma + xn = (xn - mean) / sigma igp, = where(0 <= xn) igm, = where(xn < 0) g = xn.copy() - if ga!=0: - np.put(g,igp,(-expm1(-ga*xn[igp]))/ga) + if ga != 0: + np.put(g, igp, (-expm1(-ga * xn[igp])) / ga) - if gb!=0: - np.put(g,igm,(-expm1(-gb*xn[igm]))/gb) + if gb != 0: + np.put(g, igm, (-expm1(-gb * xn[igm])) / gb) - return (g-mean2)*self.ysigma/sigma2 + self.ymean + return (g - mean2) * self.ysigma / sigma2 + self.ymean - def _gauss2dat(self,y, *yi): - if len(yi)>0: + def _gauss2dat(self, y, *yi): + if len(yi) > 0: raise ValueError('Transforming derivatives is not implemented!') ga, gb, sigma2, mean2 = self._get_par() mean = self.mean sigma = self.sigma - yn = (atleast_1d(y)-self.ymean)/self.ysigma - xn = sigma2*yn+mean2 + yn = (atleast_1d(y) - self.ymean) / self.ysigma + xn = sigma2 * yn + mean2 igp, = where(0 <= xn) igm, = where(xn < 0) - if ga!=0: - np.put(xn, igp, -log1p(-ga*xn[igp])/ga) + if ga != 0: + np.put(xn, igp, -log1p(-ga * xn[igp]) / ga) - if gb!=0: - np.put(xn, igm, -log1p(-gb*xn[igm])/gb) + if gb != 0: + np.put(xn, igm, -log1p(-gb * xn[igm]) / gb) - return sigma*xn + mean + return sigma * xn + mean def main(): import pylab - g = TrHermite(skew=0.1,kurt=3.01) + g = TrHermite(skew=0.1, kurt=3.01) g.dist2gauss() #g = TrOchi(skew=0.56) - x = np.linspace(-5,5) + x = np.linspace(-5, 5) y = g(x) - pylab.plot(np.abs(x-g.gauss2dat(y))) + pylab.plot(np.abs(x - g.gauss2dat(y))) #pylab.plot(x,y,x,x,':',g.gauss2dat(y),y,'r') pylab.show() np.disp('finito') -if __name__=='__main__': +if __name__ == '__main__': if True : # False: # import doctest doctest.testmod()