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

Small enhancements

Browse Source 
-Added trdata method to TimeSeries class
-Added TrLinear class to transform/models.py
-Added more doctest examples to spectrum/core.py
master
Per.Andreas.Brodtkorb 15 years ago
parent da96ab184a
commit acab68c0a1
8 changed files with 464 additions and 194 deletions
Show Stats Download Patch File Download Diff File

43
pywafo/src/wafo/misc.py
Unescape Escape View File

@ -596,15 +596,15 @@ def rfcfilter(x, h, method=0):
if cmpfun1(yi, fmi): if cmpfun1(yi, fmi):
z1 = -1 z1 = -1
elif cmpfun1(fpi, yi): elif cmpfun1(fpi, yi):
z1 = + 1 z1 = +1
else: else:
z1 = 0 z1 = 0
t1, y1 = (t0, y0) if z1 == 0 else (ti, yi) t1, y1 = (t0, y0) if z1 == 0 else (ti, yi)
else: 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 z1 = -1
elif (((z0 == + 1) & cmpfun2(fmi, yi)) | ((z0 == -1) & cmpfun1(fpi, yi))): elif (((z0 == +1) & cmpfun2(fmi, yi)) | ((z0 == -1) & cmpfun1(fpi, yi))):
z1 = + 1 z1 = +1
else: else:
warnings.warn('Something wrong, i=%d' % tim1) warnings.warn('Something wrong, i=%d' % tim1)
@ -614,7 +614,7 @@ def rfcfilter(x, h, method=0):
elif z1 == -1: elif z1 == -1:
#% y1 = min([y0 xi]) #% y1 = min([y0 xi])
t1, y1 = (t0, y0) if y0 < yi else (ti, yi) t1, y1 = (t0, y0) if y0 < yi else (ti, yi)
elif z1 == + 1: elif z1 == +1:
#% y1 = max([y0 xi]) #% y1 = max([y0 xi])
t1, y1 = (t0, y0) if y0 > yi else (ti, yi) 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,:] >>> tph = x1[itph,:]
>>> a = pylab.plot(x1[:,0],x1[:,1],tp[:,0],tp[:,1],'ro',tph[:,1],tph[:,1],'k.') >>> a = pylab.plot(x1[:,0],x1[:,1],tp[:,0],tp[:,1],'ro',tph[:,1],tph[:,1],'k.')
>>> pylab.close('all') >>> 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 See also
--------- ---------
@ -1031,7 +1038,7 @@ def common_shape(*args, ** kwds):
shape = kwds.get('shape') shape = kwds.get('shape')
if shape is not None: if shape is not None:
if not isinstance(shape, (list, tuple)): if not isinstance(shape, (list, tuple)):
shape = (shape, ) shape = (shape,)
shapes.append(tuple(shape)) shapes.append(tuple(shape))
if len(set(shapes)) == 1: if len(set(shapes)) == 1:
# Common case where nothing needs to be broadcasted. # Common case where nothing needs to be broadcasted.
@ -1055,7 +1062,7 @@ def common_shape(*args, ** kwds):
if len(unique) > 2: if len(unique) > 2:
# There must be at least two non-1 lengths for this axis. # There must be at least two non-1 lengths for this axis.
raise ValueError("shape mismatch: two or more arrays have " 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: elif len(unique) == 2:
# There is exactly one non-1 length. The common shape will take this # There is exactly one non-1 length. The common shape will take this
# value. # value.
@ -1109,7 +1116,7 @@ def argsreduce(condition, * args):
""" """
newargs = atleast_1d(*args) newargs = atleast_1d(*args)
if not isinstance(newargs, list): if not isinstance(newargs, list):
newargs = [newargs,] newargs = [newargs, ]
expand_arr = (condition == condition) expand_arr = (condition == condition)
return [extract(condition, arr1 * expand_arr) for arr1 in newargs] return [extract(condition, arr1 * expand_arr) for arr1 in newargs]
@ -1541,15 +1548,15 @@ def meshgrid(*xi, ** kwargs):
ndim = len(args) ndim = len(args)
s0 = (1, ) * ndim s0 = (1,) * ndim
output = [x.reshape(s0[:i] + (-1, ) + s0[i + 1::]) for i, x in enumerate(args)] output = [x.reshape(s0[:i] + (-1,) + s0[i + 1::]) for i, x in enumerate(args)]
shape = [x.size for x in output] shape = [x.size for x in output]
if indexing == 'xy': if indexing == 'xy':
# switch first and second axis # switch first and second axis
output[0].shape = (1, -1) + (1, ) * (ndim - 2) output[0].shape = (1, -1) + (1,) * (ndim - 2)
output[1].shape = (-1, 1) + (1, ) * (ndim - 2) output[1].shape = (-1, 1) + (1,) * (ndim - 2)
shape[0], shape[1] = shape[1], shape[0] shape[0], shape[1] = shape[1], shape[0]
if sparse: if sparse:
@ -1722,7 +1729,7 @@ def tranproc(x, f, x0, *xi):
xo, fo, x0 = atleast_1d(x, f, x0) xo, fo, x0 = atleast_1d(x, f, x0)
xi = atleast_1d(*xi) xi = atleast_1d(*xi)
if not isinstance(xi, list): if not isinstance(xi, list):
xi = [xi,] xi = [xi, ]
N = len(xi) # N = number of derivatives N = len(xi) # N = number of derivatives
nmax = ceil((xo.ptp()) * 10 ** (7. / max(N, 1))) 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) 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. #% Transform X with the derivatives of f.
fxder = zeros((N, x0.size)) 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. 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 = [(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) fxder[k] = tranproc(fder[0], fder[1], x0)
#% Calculate the transforms of the derivatives of X. # Calculate the transforms of the derivatives of X.
#% First time derivative of y: y1 = f'(x)*x1 # First time derivative of y: y1 = f'(x)*x1
y1 = fxder[0] * xi[0] y1 = fxder[0] * xi[0]
y.append(y1) y.append(y1)

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

@ -14,14 +14,15 @@
from __future__ import division from __future__ import division
from wafo.transform.core import TrData 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.interpolate.interpolate import interp1d
from scipy.integrate.quadrature import cumtrapz from scipy.integrate.quadrature import cumtrapz
from wafo.interpolate import SmoothSpline
import warnings import warnings
import numpy as np import numpy as np
from numpy import (inf, pi, zeros, ones, sqrt, where, log, exp, sin, arcsin, mod, finfo, interp, #@UnresolvedImport 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 polyfit, r_, nonzero, cumsum, ravel, size, isnan, nan, floor, ceil, diff, array) #@UnresolvedImport
from numpy.fft import fft from numpy.fft import fft
from numpy.random import randn from numpy.random import randn
@ -359,13 +360,13 @@ class LevelCrossings(WafoData):
if ng == 1: if ng == 1:
gvar = opt.gvar * ones(ncr) gvar = opt.gvar * ones(ncr)
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)) ng = len(atleast_1d(opt.cvar))
if ng == 1: if ng == 1:
cvar = opt.cvar * ones(ncr) cvar = opt.cvar * ones(ncr)
else: 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 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
@ -401,9 +402,9 @@ class LevelCrossings(WafoData):
#u0 = interp1q(cros(:,2),cros(:,1),.5) #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 # NB! the smooth function does not always extrapolate well outside the edges
# causing poor estimate of g # causing poor estimate of g
# We may alleviate this problem by: forcing the extrapolation # 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 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) 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: if opt.chkder:
for ix in range(5): for ix in range(5):
@ -603,8 +604,6 @@ class TurningPoints(WafoData):
---------------- ----------------
data : array_like data : array_like
args : vector for 1D args : vector for 1D
''' '''
def __init__(self, *args, **kwds): def __init__(self, *args, **kwds):
super(TurningPoints, self).__init__(*args, **kwds) super(TurningPoints, self).__init__(*args, **kwds)
@ -857,7 +856,207 @@ class TimeSeries(WafoData):
fact = 2.0 * pi fact = 2.0 * pi
w = fact * f w = fact * f
return _wafospec.SpecData1D(S / fact, w) 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): def turning_points(self, h=0.0, wavetype=None):
''' '''
Return turning points (tp) from data, optionally rainflowfiltered. Return turning points (tp) from data, optionally rainflowfiltered.

149
pywafo/src/wafo/spectrum/core.py
Unescape Escape View File

@ -1,5 +1,6 @@
from __future__ import division from __future__ import division
from scipy.misc.ppimport import ppimport from scipy.misc.ppimport import ppimport
from wafo.objects import mat2timeseries, TimeSeries
import warnings import warnings
import numpy as np import numpy as np
from numpy import (pi, inf, meshgrid, zeros, ones, where, nonzero, #@UnresolvedImport 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 dispersion_relation import w2k #, k2w
from wafo.wafodata import WafoData, now 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: try:
from wafo.gaussian import Rind from wafo.gaussian import Rind
except ImportError: except ImportError:
@ -28,7 +30,8 @@ except ImportError:
warnings.warn('Compile the c_library.pyd again!') warnings.warn('Compile the c_library.pyd again!')
c_library = None c_library = None
from wafo.transform import TrData #from wafo.transform import TrData
from wafo.transform.models import TrLinear
from wafo.plotbackend import plotbackend from wafo.plotbackend import plotbackend
@ -190,7 +193,7 @@ class SpecData1D(WafoData):
self.freqtype = 'w' self.freqtype = 'w'
self.angletype = '' self.angletype = ''
self.h = inf self.h = inf
self.tr = None self.tr = None #TrLinear()
self.phi = 0.0 self.phi = 0.0
self.v = 0.0 self.v = 0.0
self.norm = False self.norm = False
@ -572,7 +575,9 @@ class SpecData1D(WafoData):
#Fs = 2*freq(end)+eps; % sampling frequency #Fs = 2*freq(end)+eps; % sampling frequency
for ix in xrange(max_sim): 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); #%x2(:,2:end) = x2(:,2:end) -x1(:,2:end);
S2 = dat2spec(x2, L) S2 = dat2spec(x2, L)
S1 = dat2spec(x1, L) S1 = dat2spec(x1, L)
@ -727,9 +732,10 @@ class SpecData1D(WafoData):
if self.tr is None: if self.tr is None:
y = linspace(-5, 5, 513) g = TrLinear(var=m[0])
#y = linspace(-5, 5, 513)
#g = _wafotransform. #g = _wafotransform.
g = TrData(y, sqrt(m[0]) * y) #g = TrData(y, sqrt(m[0]) * y)
else: else:
g = self.tr g = self.tr
@ -1053,7 +1059,7 @@ class SpecData1D(WafoData):
xder[:, 0] = x[:, 0] xder[:, 0] = x[:, 0]
if spec.tr is not None: if spec.tr is not None:
print(' Transforming data.') #print(' Transforming data.')
g = spec.tr g = spec.tr
if derivative: if derivative:
for i in range(cases): 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) # 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', 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) Simulates a Randomized 2nd order non-linear wave X(t)
@ -1142,6 +1148,21 @@ class SpecData1D(WafoData):
Example 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)<sa)
np =100; dt = .2 np =100; dt = .2
[x1, x2] = spec2nlsdat(jonswap,np,dt) [x1, x2] = spec2nlsdat(jonswap,np,dt)
waveplot(x1,'r',x2,'g',1,1) waveplot(x1,'r',x2,'g',1,1)
@ -1281,8 +1302,8 @@ class SpecData1D(WafoData):
#if isempty(nmin),nmin = 2end % Must always be greater than 1 #if isempty(nmin),nmin = 2end % Must always be greater than 1
f_limit_up = df * nmax f_limit_up = df * nmax
f_limit_lo = df * nmin f_limit_lo = df * nmin
if verbose:
print('2nd order frequency Limits = %g,%g' % (f_limit_lo, f_limit_up)) print('2nd order frequency Limits = %g,%g' % (f_limit_lo, f_limit_up))
@ -1298,7 +1319,7 @@ class SpecData1D(WafoData):
## % 1'st order + 2'nd order component. ## % 1'st order + 2'nd order component.
## x2(:,2:end) =x(:,2:end)+ real(x2s(1:np,:))+real(x2d(1:np,:)) ## x2(:,2:end) =x(:,2:end)+ real(x2s(1:np,:))+real(x2d(1:np,:))
## else ## else
amp = amp.T amp = np.array(amp.T).ravel()
rvec, ivec = c_library.disufq(amp.real, amp.imag, w, kw, water_depth, g, nmin, nmax, cases, ns) rvec, ivec = c_library.disufq(amp.real, amp.imag, w, kw, water_depth, g, nmin, nmax, cases, ns)
svec = rvec + 1J * ivec svec = rvec + 1J * ivec
@ -1325,8 +1346,8 @@ class SpecData1D(WafoData):
composed of letters ['mvsk'] specifying which moments to compute: composed of letters ['mvsk'] specifying which moments to compute:
'm' = mean, 'm' = mean,
'v' = variance, 'v' = variance,
's' = (Fisher's) skew, 's' = skewness,
'k' = (Fisher's) kurtosis. 'k' = (Pearson's) kurtosis.
method : string method : string
'approximate' method due to Marthinsen & Winterstein (default) 'approximate' method due to Marthinsen & Winterstein (default)
'eigenvalue' method due to Kac and Siegert 'eigenvalue' method due to Kac and Siegert
@ -1358,16 +1379,15 @@ class SpecData1D(WafoData):
-------- --------
#Simulate a Transformed Gaussian process: #Simulate a Transformed Gaussian process:
>>> import wafo.spectrum.models as sm >>> 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() >>> S = Sj.tospecdata()
>>> me, va, sk, ku = S.stats_nl(moments='mvsk') >>> me, va, sk, ku = S.stats_nl(moments='mvsk')
>>> g = wtm.TrHermite(mean=me, sigma=Hs/4, skew=sk, kurt=ku, ysigma=Hs/4)
>>> ys = S.sim(15000) # Simulated in the Gaussian world
Hm0=7;Tp=11 >>> xs = g.gauss2dat(ys[:,1]) # Transformed to the real world
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
See also See also
--------- ---------
@ -1483,7 +1503,7 @@ class SpecData1D(WafoData):
## skew = sum((6*C2+8*E2).*E)/sa^3 % skewness ## skew = sum((6*C2+8*E2).*E)/sa^3 % skewness
## kurt = 3+48*sum((C2+E2).*E2)/sa^4 % kurtosis ## kurt = 3+48*sum((C2+E2).*E2)/sa^4 % kurtosis
return output 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. 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. 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. If 95% of TEST1 is less than TEST0 then X(t) is not Gaussian at a 5% level.
Example: Example:
Hm0 = 7; -------
S0 = jonswap([],Hm0); g=ochitr([],[Hm0/4]); S=S0; >>> import wafo.spectrum.models as sm
S.tr=g;S.tr(:,2)=g(:,2)*Hm0/4; >>> import wafo.transform.models as wtm
xs = spec2sdat(S,2^13); >>> import wafo.objects as wo
[g0 t0] = dat2tr(xs); >>> Hs = 7
t1 = testgaussian(S0,[2^13 50],t0); >>> 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 See also cov2sdat, dat2tr, troptset
''' '''
@ -1542,49 +1573,43 @@ class SpecData1D(WafoData):
# #
# opt = troptset(opt,'multip',1) # opt = troptset(opt,'multip',1)
if test0 is None: plotflag=0 if test0 is None else 1
plotflag=0
else:
plotflag=1
if cases>50: if cases>50:
print(' ... be patient this may take a while') 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() acf = self.tocovdata()
#R = spec2cov(S); #R = spec2cov(S);
test1 = []
for ix in range(rep): for ix in range(rep):
xs = acf.sim(ns=ns, cases=Nstep) 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]); #xs = cov2sdat(R,[ns Nstep]);
[g, tmp] = dat2tr(xs,method, **opt); #[g, tmp] = dat2tr(xs,method, **opt);
#test1 = [test1; tmp(:)] #test1 = [test1; tmp(:)]
print('finished %d of %d ' % (ix,rep) ) print('finished %d of %d ' % (ix+1,rep) )
if rep>1: if rep>1:
xs = acf.sim(ns=ns, cases=rem(cases,rep)) xs = acf.sim(ns=ns, cases=np.remainder(cases,rep))
[g, tmp] = dat2tr(xs,method,**opt); for iy in range(1, xs.shape[-1]):
#test1 = [test1; tmp(:)]; ts = TimeSeries(xs[:, iy], xs[:, 0].ravel())
g, tmp = ts.trdata(method, **opt)
test1.append(g.dist2gauss())
# if (nargout>0 || plotflag==0),
# test2=test1; if plotflag:
# end plotbackend.plot(test1,'o')
# plotbackend.plot([1, cases], [test0, test0],'--')
#
# if plotflag plotbackend.ylabel('e(g(u)-u)')
# plot(test1,'o'),hold on plotbackend.xlabel('Simulation number')
# if 1 return test1
# plot([1 cases],test0*[1 1],'--'),
# end
# hold off
# ylabel('e(g(u)-u)')
# xlabel('Simulation number')
# end
def moment(self, nr=2, even=True, j=0): def moment(self, nr=2, even=True, j=0):
''' Calculates spectral moments from spectrum ''' Calculates spectral moments from spectrum
@ -1647,8 +1672,8 @@ class SpecData1D(WafoData):
if self.freqtype in ['f', 'w']: if self.freqtype in ['f', 'w']:
vari = 't' vari = 't'
if self.freqtype == 'f': if self.freqtype == 'f':
f = 2. * pi * f f = 2. * pi * f
S = S / (2. * pi) S = S / (2. * pi)
else: else:
vari = 'x' vari = 'x'
S1 = abs(S) ** (j + 1.) S1 = abs(S) ** (j + 1.)
@ -2326,7 +2351,7 @@ class SpecData2D(WafoData):
##% By es 27.08.1999 ##% By es 27.08.1999
pi = pi
two_dim_spectra = ['dir', 'encdir', 'k2d'] two_dim_spectra = ['dir', 'encdir', 'k2d']
if self.type not in two_dim_spectra: if self.type not in two_dim_spectra:
raise ValueError('Unknown 2D spectrum type!') raise ValueError('Unknown 2D spectrum type!')

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

@ -55,7 +55,7 @@ from numpy import (inf, atleast_1d, newaxis, any, minimum, maximum, array, #@Unr
from dispersion_relation import w2k from dispersion_relation import w2k
#ppimport.enable() #ppimport.enable()
#_wafospectrum = ppimport.ppimport('wafo.spectrum') #_wafospectrum = ppimport.ppimport('wafo.spectrum')
from core import SpecData1D from wafo.spectrum import SpecData1D
sech = lambda x: 1.0 / cosh(x) sech = lambda x: 1.0 / cosh(x)
eps = finfo(float).eps eps = finfo(float).eps
@ -542,9 +542,9 @@ class Jonswap(ModelSpectrum):
outsideJonswapRange = Tp > 5 * sqrt(Hm0) or Tp < 3.6 * sqrt(Hm0) outsideJonswapRange = Tp > 5 * sqrt(Hm0) or Tp < 3.6 * sqrt(Hm0)
if outsideJonswapRange: if outsideJonswapRange:
txt0 = ''' 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. The validity of the spectral density is questionable.
''' ''' % (Hm0, Tp)
warnings.warn(txt0) warnings.warn(txt0)
if gam < 1 or 7 < gam: if gam < 1 or 7 < gam:

4
pywafo/src/wafo/stats/core.py
Unescape Escape View File

@ -8,7 +8,7 @@ __all__ = ['edf']
def edf(x, method=2): def edf(x, method=2):
''' '''
Returns EDF Empirical Distribution Function. Returns Empirical Distribution Function (EDF).
Parameters Parameters
---------- ----------
@ -48,7 +48,7 @@ def edf(x, method=2):
def edfcnd(x, c=None, 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 Parameters
---------- ----------

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

@ -40,7 +40,7 @@ arr = asarray
all = alltrue all = alltrue
def chi2isf(p,df): def chi2isf(p, df):
return special.chdtri(df, p) return special.chdtri(df, p)
def chi2sf(x, df): def chi2sf(x, df):
return special.chdtrc(df, x) return special.chdtrc(df, x)
@ -92,7 +92,7 @@ class rv_frozen(object):
def __init__(self, dist, *args, **kwds): def __init__(self, dist, *args, **kwds):
self.dist = dist self.dist = dist
loc0, scale0 = map(kwds.get, ['loc', 'scale']) 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) args, loc0, scale0 = dist.fix_loc_scale(args, loc0, scale0)
self.par = args + (loc0, scale0) self.par = args + (loc0, scale0)
else: # rv_discrete else: # rv_discrete

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

@ -48,8 +48,8 @@ class TrCommon(object):
self.ymean = kwds.get('ymean', 0e0) self.ymean = kwds.get('ymean', 0e0)
self.ysigma = kwds.get('ysigma', 1e0) self.ysigma = kwds.get('ysigma', 1e0)
def __call__(self, x): def __call__(self, x, *xi):
return self._dat2gauss(x) return self._dat2gauss(x, *xi)
def dist2gauss(self, x=None, xnmin=-5, xnmax=5, n=513): def dist2gauss(self, x=None, xnmin=-5, xnmax=5, n=513):
""" """
@ -158,6 +158,7 @@ class TrData(WafoData, TrCommon):
1 1
>>> g.sigma >>> g.sigma
5 5
>>> g.dat2gauss(1,2,3)
Check that the departure from a Gaussian model is zero Check that the departure from a Gaussian model is zero

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

@ -15,20 +15,18 @@ TrOchi
# Licence: <your licence> # Licence: <your licence>
#------------------------------------------------------------------------------- #-------------------------------------------------------------------------------
#!/usr/bin/env python #!/usr/bin/env python
from __future__ import division
from scipy.optimize import brentq from scipy.optimize import brentq
from numpy import (sqrt, atleast_1d, abs, imag, sign, where, cos, arccos, ceil, #@UnresolvedImport from numpy import (sqrt, atleast_1d, abs, imag, sign, where, cos, arccos, ceil, #@UnresolvedImport
expm1, log1p, pi) #@UnresolvedImport expm1, log1p, pi) #@UnresolvedImport
import numpy as np import numpy as np
import warnings import warnings
from core import TrCommon from core import TrCommon
__all__=['TrHermite','TrOchi'] __all__ = ['TrHermite', 'TrLinear', 'TrOchi']
_example = ''' _example = '''
>>> std = 7./4 >>> std = 7./4
>>> g = <generic>(sigma=std, ysigma=std) >>> g = <generic>(sigma=std, ysigma=std)
>>> g.dist2gauss()
3.9858776379926808
Simulate a Transformed Gaussian process: Simulate a Transformed Gaussian process:
>>> import numpy as np >>> import numpy as np
@ -43,8 +41,9 @@ _example = '''
>>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world >>> xs = g2.gauss2dat(ys[:,1:]) # Transformed to the real world
''' '''
class TrHermite(TrCommon): class TrHermite(TrCommon):
__doc__ = TrCommon.__doc__.replace('<generic>','Hermite') + """ __doc__ = TrCommon.__doc__.replace('<generic>', 'Hermite') + """
pardef : scalar, integer pardef : scalar, integer
1 Winterstein et. al. (1994) parametrization [1]_ (default) 1 Winterstein et. al. (1994) parametrization [1]_ (default)
2 Winterstein (1988) parametrization [2]_ 2 Winterstein (1988) parametrization [2]_
@ -79,7 +78,9 @@ class TrHermite(TrCommon):
Example: Example:
-------- --------
""" + _example.replace('<generic>','TrHermite') + """ """ + _example.replace('<generic>', 'TrHermite') + """
>>> g.dist2gauss()
3.9858776379926808
See also See also
-------- --------
@ -96,9 +97,9 @@ class TrHermite(TrCommon):
'Nonlinear vibration models for extremes and fatigue.' 'Nonlinear vibration models for extremes and fatigue.'
J. Engng. Mech., ASCE, Vol 114, No 10, pp 1772-1790 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) super(TrHermite, self).__init__(*args, **kwds)
self.pardef = kwds.get('pardef',1) self.pardef = kwds.get('pardef', 1)
self._c3 = None self._c3 = None
self._c4 = None self._c4 = None
self._forward = None self._forward = None
@ -108,32 +109,32 @@ class TrHermite(TrCommon):
def _poly_par_from_stats(self): def _poly_par_from_stats(self):
skew = self.skew skew = self.skew
ga2 = self.kurt-3.0 ga2 = self.kurt - 3.0
if ga2 <= 0: if ga2 <= 0:
self._c4 = ga2/24. self._c4 = ga2 / 24.
self._c3 = skew/6. self._c3 = skew / 6.
elif self.pardef == 2: elif self.pardef == 2:
# Winterstein 1988 parametrization # 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') warnings.warn('Kurtosis too low compared to the skewness')
self._c4 = (sqrt(1.+1.5*ga2)-1.)/18. self._c4 = (sqrt(1. + 1.5 * ga2) - 1.) / 18.
self._c3 = skew/(6.*(1+6.*self._c4)) self._c3 = skew / (6. * (1 + 6. * self._c4))
else: else:
# Winterstein et. al. 1994 parametrization intended to # Winterstein et. al. 1994 parametrization intended to
# apply for the range: 0 <= ga2 < 12 and 0<= skew^2 < 2*ga2/3 # 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') warnings.warn('Kurtosis too low compared to the skewness')
if (ga2 < 0) or (12 < ga2): if (ga2 < 0) or (12 < ga2):
warnings.warn('Kurtosis must be between 0 and 12') 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.: if ga2 == 0.:
self._c4 = 0.0 self._c4 = 0.0
else: else:
c41 = (1.-1.43*skew**2./ga2)**(1.-0.1*(ga2+3.)**0.8) 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 self._c4 = 0.1 * ((1. + 1.25 * ga2) ** (1. / 3.) - 1.) * c41
if not np.isfinite(self._c3) or not np.isfinite(self._c4): if not np.isfinite(self._c3) or not np.isfinite(self._c4):
raise ValueError('Unable to calculate the polynomial') raise ValueError('Unable to calculate the polynomial')
@ -150,31 +151,31 @@ class TrHermite(TrCommon):
c4 = self._c4 c4 = self._c4
ma = self.mean ma = self.mean
sa = self.sigma sa = self.sigma
if abs(c4)<sqrt(eps): if abs(c4) < sqrt(eps):
c4 = 0.0 c4 = 0.0
#gdef = self.kurt-3.0 #gdef = self.kurt-3.0
if self.kurt<3.0: if self.kurt < 3.0:
p = np.poly1d([-c4, -c3, 1.+3.*c4, c3]) # forward, g p = np.poly1d([-c4, -c3, 1. + 3. * c4, c3]) # forward, g
self._forward = p self._forward = p
self._backward = None self._backward = None
else: else:
Km1 = np.sqrt(1.+2.*c3**2+6*c4**2) Km1 = np.sqrt(1. + 2. * c3 ** 2 + 6 * c4 ** 2)
p = np.poly1d(np.r_[c4, c3, 1.-3.*c4, -c3]/Km1) # backward G p = np.poly1d(np.r_[c4, c3, 1. - 3. * c4, -c3] / Km1) # backward G
self._forward = None self._forward = None
self._backward = p self._backward = p
#% Check if it is a strictly increasing function. #% Check if it is a strictly increasing function.
dp = p.deriv(m=1) #% Derivative dp = p.deriv(m=1) #% Derivative
r = dp.r #% Find roots of the derivative r = dp.r #% Find roots of the derivative
r = r[where(abs(imag(r))<eps)] # Keep only real roots r = r[where(abs(imag(r)) < eps)] # Keep only real roots
if r.size>0: if r.size > 0:
# Compute where it is possible to invert the polynomial # Compute where it is possible to invert the polynomial
if self.kurt<3.: if self.kurt < 3.:
self._x_limit = r self._x_limit = r
else: else:
self._x_limit = sa*p(r)+ma self._x_limit = sa * p(r) + ma
txt1 = ''' txt1 = '''
The polynomial is not a strictly increasing function. The polynomial is not a strictly increasing function.
The derivative of g(x) is infinite at x = %g''' % self._x_limit 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): def check_forward(self, x):
if not (self._x_limit is None): if not (self._x_limit is None):
x00 = self._x_limit 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 cdef = 1
else: 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 errtxt = 'Unable to invert the polynomial \n %s' % txt2
raise ValueError(errtxt) raise ValueError(errtxt)
np.disp('However, successfully inverted the polynomial\n %s' % txt2) np.disp('However, successfully inverted the polynomial\n %s' % txt2)
def _dat2gauss(self,x, *xi): def _dat2gauss(self, x, *xi):
if len(xi)>0: if len(xi) > 0:
raise ValueError('Transforming derivatives is not implemented!') raise ValueError('Transforming derivatives is not implemented!')
xn = atleast_1d(x) xn = atleast_1d(x)
self.check_forward(xn) self.check_forward(xn)
xn = (xn-self.mean)/self.sigma xn = (xn - self.mean) / self.sigma
if self._forward is None: if self._forward is None:
#Inverting the polynomial #Inverting the polynomial
yn = self._poly_inv(self._backward,xn) yn = self._poly_inv(self._backward, xn)
else: else:
yn = self._forward(xn) yn = self._forward(xn)
return yn*self.ysigma + self.ymean return yn * self.ysigma + self.ymean
def _gauss2dat(self, y, *yi): def _gauss2dat(self, y, *yi):
if len(yi)>0: if len(yi) > 0:
raise ValueError('Transforming derivatives is not implemented!') 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) #self.check_forward(y)
if self._backward is None: if self._backward is None:
#% Inverting the polynomial #% Inverting the polynomial
#%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
xn = self._poly_inv(self._forward,yn) xn = self._poly_inv(self._forward, yn)
else: else:
xn = self._backward(yn) xn = self._backward(yn)
return self.sigma*xn + self.mean return self.sigma * xn + self.mean
def _poly_inv(self, p, xn): def _poly_inv(self, p, xn):
''' '''
Invert polynomial Invert polynomial
''' '''
if p.order<2: if p.order < 2:
return xn return xn
elif p.order==2: elif p.order == 2:
# Quadratic: Solve a*u**2+b*u+c = xn # Quadratic: Solve a*u**2+b*u+c = xn
coefs = p.coeffs coefs = p.coeffs
a = coefs[0] a = coefs[0]
b = coefs[1] b = coefs[1]
c = coefs[2]-xn c = coefs[2] - xn
t = 0.5*(b+sign(b)*sqrt(b**2-4*a*c)) t = 0.5 * (b + sign(b) * sqrt(b ** 2 - 4 * a * c))
#so1 = t/a # largest solution #so1 = t/a # largest solution
so2 = -c/t # smallest solution so2 = -c / t # smallest solution
return so2 return so2
elif p.order==3: elif p.order == 3:
# Solve # Solve
# K*(c4*u^3+c3*u^2+(1-3*c4)*u-c3) = xn = (x-ma)/sa # 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 # -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] a = coefs[0]
b = coefs[1] 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 #% 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 #p1 = (b-a**2/3)/3
#q0 = (c + x0*(2.*x0/3.-b))/2. #q0 = (c + x0*(2.*x0/3.-b))/2.
#q0 = x0**3 -a*b/6 +c/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 ## # z^3+3*p1*z+2*q0=0
## c3 = self._c3 ## c3 = self._c3
@ -272,24 +273,59 @@ class TrHermite(TrCommon):
#q0 = x0**3-1.5*b1*(x0+xn) #q0 = x0**3-1.5*b1*(x0+xn)
if not (self._x_limit is None): # % Three real roots if not (self._x_limit is None): # % Three real roots
d = sqrt(-p1) d = sqrt(-p1)
theta1 = arccos(-q0/d**3)/3 theta1 = arccos(-q0 / d ** 3) / 3
th2 = np.r_[0, -2*pi/3, 2*pi/3] th2 = np.r_[0, -2 * pi / 3, 2 * pi / 3]
x1 = abs(2*d*cos(theta1[ceil(len(xn)/2)] + th2)-x0) x1 = abs(2 * d * cos(theta1[ceil(len(xn) / 2)] + th2) - x0)
ix = x1.argmin() # % choose the smallest solution 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 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 #% Find the real root of the monic polynomial
A0 = (q1-q0)**(1./3.) A0 = (q1 - q0) ** (1. / 3.)
B0 = -(q1+q0)**(1./3.) B0 = -(q1 + q0) ** (1. / 3.)
return A0+B0-x0 #% real root return A0 + B0 - x0 #% real root
#%% The other complex roots are given by #%% 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
#%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('<generic>', '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('<generic>', '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): class TrOchi(TrCommon):
__doc__ = TrCommon.__doc__.replace('<generic>','Ochi') + """ __doc__ = TrCommon.__doc__.replace('<generic>', 'Ochi') + """
Description Description
----------- -----------
@ -322,8 +358,10 @@ class TrOchi(TrCommon):
Example Example
------- -------
""" + _example.replace('<generic>','TrOchi') + """ """ + _example.replace('<generic>', 'TrOchi') + """
>>> g.dist2gauss()
5.9322684525265501
See also See also
-------- --------
spec2skew, hermitetr, lc2tr, dat2tr spec2skew, hermitetr, lc2tr, dat2tr
@ -347,13 +385,13 @@ class TrOchi(TrCommon):
def _par_from_stats(self): def _par_from_stats(self):
skew = self.skew skew = self.skew
if abs(skew)>2.82842712474619: if abs(skew) > 2.82842712474619:
raise ValueError('Skewness must be less than 2.82842') raise ValueError('Skewness must be less than 2.82842')
mean1 = self.mean mean1 = self.mean
sigma1 = self.sigma sigma1 = self.sigma
if skew==0: if skew == 0:
self._phat = [sigma1, mean1, 0, 0, 1, 0] self._phat = [sigma1, mean1, 0, 0, 1, 0]
return return
@ -367,22 +405,22 @@ class TrOchi(TrCommon):
# Set up the 2D non-linear equations for a and sig2^2: # 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 ]' # 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: # 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. a1 = 1. # % Start interval where sig2^2 is located.
a2 = 2. a2 = 2.
sig22 = brentq(g2,a1,a2) #% smallest solution for sig22 sig22 = brentq(g2, a1, a2) #% smallest solution for sig22
a = sign(skew)*sqrt(abs(sig22-1)/2)/sig22 a = sign(skew) * sqrt(abs(sig22 - 1) / 2) / sig22
gam_a = 1.28*a gam_a = 1.28 * a
gam_b = 3*a gam_b = 3 * a
sigma2 = sqrt(sig22) sigma2 = sqrt(sig22)
#% Solve the following 2nd order equation to obtain ma2 #% Solve the following 2nd order equation to obtain ma2
#% a*(sig2^2+ma2^2)+ma2 = 0 #% a*(sig2^2+ma2^2)+ma2 = 0
my2 = (-1.-sqrt(1.-4.*a**2*sig22))/a #% Largest mean my2 = (-1. - sqrt(1. - 4. * a ** 2 * sig22)) / a #% Largest mean
mean2 = a*sig22/my2 #% choose the smallest mean mean2 = a * sig22 / my2 #% choose the smallest mean
self._phat = [sigma1, mean1, gam_a, gam_b, sigma2, mean2] self._phat = [sigma1, mean1, gam_a, gam_b, sigma2, mean2]
return return
@ -391,8 +429,8 @@ class TrOchi(TrCommon):
''' '''
Returns ga, gb, sigma2, mean2 Returns ga, gb, sigma2, mean2
''' '''
if (self._phat is None or self.sigma!=self._phat[0] if (self._phat is None or self.sigma != self._phat[0]
or self.mean!=self._phat[1]): or self.mean != self._phat[1]):
self._par_from_stats() self._par_from_stats()
#sigma1 = self._phat[0] #sigma1 = self._phat[0]
#mean1 = self._phat[1] #mean1 = self._phat[1]
@ -402,63 +440,63 @@ class TrOchi(TrCommon):
mean2 = self._phat[5] mean2 = self._phat[5]
return ga, gb, sigma2, mean2 return ga, gb, sigma2, mean2
def _dat2gauss(self,x, *xi): def _dat2gauss(self, x, *xi):
if len(xi)>0: if len(xi) > 0:
raise ValueError('Transforming derivatives is not implemented!') raise ValueError('Transforming derivatives is not implemented!')
ga, gb, sigma2, mean2 = self._get_par() ga, gb, sigma2, mean2 = self._get_par()
mean = self.mean mean = self.mean
sigma = self.sigma sigma = self.sigma
xn = atleast_1d(x) xn = atleast_1d(x)
xn = (xn-mean)/sigma xn = (xn - mean) / sigma
igp, = where(0 <= xn) igp, = where(0 <= xn)
igm, = where(xn < 0) igm, = where(xn < 0)
g = xn.copy() g = xn.copy()
if ga!=0: if ga != 0:
np.put(g,igp,(-expm1(-ga*xn[igp]))/ga) np.put(g, igp, (-expm1(-ga * xn[igp])) / ga)
if gb!=0: if gb != 0:
np.put(g,igm,(-expm1(-gb*xn[igm]))/gb) 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): def _gauss2dat(self, y, *yi):
if len(yi)>0: if len(yi) > 0:
raise ValueError('Transforming derivatives is not implemented!') raise ValueError('Transforming derivatives is not implemented!')
ga, gb, sigma2, mean2 = self._get_par() ga, gb, sigma2, mean2 = self._get_par()
mean = self.mean mean = self.mean
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 + mean2
igp, = where(0 <= xn) igp, = where(0 <= xn)
igm, = where(xn < 0) igm, = where(xn < 0)
if ga!=0: if ga != 0:
np.put(xn, igp, -log1p(-ga*xn[igp])/ga) np.put(xn, igp, -log1p(-ga * xn[igp]) / ga)
if gb!=0: if gb != 0:
np.put(xn, igm, -log1p(-gb*xn[igm])/gb) np.put(xn, igm, -log1p(-gb * xn[igm]) / gb)
return sigma*xn + mean return sigma * xn + mean
def main(): def main():
import pylab import pylab
g = TrHermite(skew=0.1,kurt=3.01) g = TrHermite(skew=0.1, kurt=3.01)
g.dist2gauss() g.dist2gauss()
#g = TrOchi(skew=0.56) #g = TrOchi(skew=0.56)
x = np.linspace(-5,5) x = np.linspace(-5, 5)
y = g(x) 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.plot(x,y,x,x,':',g.gauss2dat(y),y,'r')
pylab.show() pylab.show()
np.disp('finito') np.disp('finito')
if __name__=='__main__': if __name__ == '__main__':
if True : # False: # if True : # False: #
import doctest import doctest
doctest.testmod() doctest.testmod()

Write Preview
Loading…
Cancel
Save