Small bugfixes

master
Per.Andreas.Brodtkorb 15 years ago
parent d0b1e917ce
commit d0c62a08fa

@ -1,6 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<projectDescription>
<name>PYWAFO</name>
<name>google_pywafo</name>
<comment></comment>
<projects>
</projects>

@ -3,7 +3,7 @@
<pydev_project>
<pydev_pathproperty name="org.python.pydev.PROJECT_SOURCE_PATH">
<path>/PYWAFO/src</path>
<path>/google_pywafo/src</path>
</pydev_pathproperty>
<pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.6</pydev_property>
<pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Default</pydev_property>

@ -17,6 +17,7 @@ python setup.py sdist bdist_wininst upload --show-response
import os, sys
sys.argv.append("develop")
#sys.argv.append("install")
DISTUTILS_DEBUG = True
pkg_name = 'wafo'
root_dir = os.path.join('src',pkg_name)
@ -44,7 +45,9 @@ testscripts = [os.path.join(subtst, f) for subtst in subtests
f.endswith('.old') or f.endswith('.bak'))]
datadir = 'data'
datafiles = [os.path.join(datadir, f) for f in os.listdir(os.path.join(root_dir, datadir))
if not (f.endswith('.py') or f.endswith('test') )]
if not (f.startswith('.') or f.endswith('~') or
f.endswith('.old') or f.endswith('.bak') or
f.endswith('.py') or f.endswith('test') )]
libs = [f for f in os.listdir(os.path.join(root_dir)) if f.endswith('.pyd') ]
packagedata = testscripts + datafiles + libs #['c_library.pyd'] #,'disufq1.c','diffsumfunq.pyd','diffsumfunq.pyf','findrfc.c','rfc.pyd','rfc.pyf']

@ -1,7 +1,12 @@
Metadata-Version: 1.0
Name: WAFO
Name: wafo
Version: 0.11
Summary:
Summary: UNKNOWN
Home-page: http://www.maths.lth.se/matstat/wafo/
Author: WAFO-group
Author-email: wafo@maths.lth.se
License: GPL
Description:
WAFO
=====
WAFO is a toolbox Python routines for statistical analysis and simulation of random waves and random loads.
@ -84,9 +89,4 @@ Statistics
- Reported bugs.
- List of publications related to WAFO.
Home-page: http://www.maths.lth.se/matstat/wafo/
Author: WAFO-group
Author-email: wafo@maths.lth.se
License: GPL
Description: UNKNOWN
Platform: UNKNOWN

@ -1,5 +1 @@
wafo\spectrum
wafo\covariance
wafo
wafo\data
wafo\transform

Binary file not shown.

@ -19,7 +19,7 @@ from numpy.lib.shape_base import vstack
from numpy.lib.function_base import linspace
import polynomial as pl
class PPform1(object):
class PPform(object):
"""The ppform of the piecewise polynomials is given in terms of coefficients
and breaks. The polynomial in the ith interval is
x_{i} <= x < x_{i+1}

@ -9,7 +9,7 @@
# Licence: LGPL
#-------------------------------------------------------------------------------
#!/usr/bin/env python
#import numpy as np
import numpy as np
from scipy.special import gamma
from numpy import pi, atleast_2d #@UnresolvedImport
from misc import tranproc, trangood
@ -112,29 +112,29 @@ class kde(object):
if d != self.d:
if d == 1 and m == self.d:
# points was passed in as a row vector
points = reshape(points, (self.d, 1))
points = np.reshape(points, (self.d, 1))
m = 1
else:
msg = "points have dimension %s, dataset has dimension %s" % (d,
self.d)
raise ValueError(msg)
result = zeros((m,), points.dtype)
result = np.zeros((m,), points.dtype)
if m >= self.n:
# there are more points than data, so loop over data
for i in range(self.n):
diff = self.dataset[:,i,newaxis] - points
tdiff = dot(self.inv_cov, diff)
energy = sum(diff*tdiff,axis=0)/2.0
result += exp(-energy)
diff = self.dataset[:,i, np.newaxis] - points
tdiff = np.dot(self.inv_cov, diff)
energy = np.sum(diff*tdiff,axis=0)/2.0
result += np.exp(-energy)
else:
# loop over points
for i in range(m):
diff = self.dataset - points[:,i,newaxis]
tdiff = dot(self.inv_cov, diff)
diff = self.dataset - points[:,i,np.newaxis]
tdiff = np.dot(self.inv_cov, diff)
energy = sum(diff*tdiff,axis=0)/2.0
result[i] = sum(exp(-energy),axis=0)
result[i] = np.sum(np.exp(-energy),axis=0)
result /= self._norm_factor

@ -1779,7 +1779,7 @@ def tranproc(x, f, x0, * xi):
y.append(y4)
if N > 4:
warnings.warn('Transformation of derivatives of order>4 not supported.')
return y #0,y1,y2,y3,y4
return y #y0,y1,y2,y3,y4
def test_common_shape():

@ -13,10 +13,14 @@
from __future__ import division
from wafo.transform.core import TrData
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, #@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
polyfit, r_, nonzero, cumsum, ravel, size, isnan, nan, floor, ceil, diff, array) #@UnresolvedImport
from numpy.fft import fft
@ -26,10 +30,11 @@ from pylab import stineman_interp
from matplotlib.mlab import psd
import scipy.signal
from scipy.special import erf
from scipy.special import erf, ndtri
from wafo.misc import (nextpow2, findtp, findtc, findcross, sub_dict_select,
ecross, JITImport)
ecross, JITImport, DotDict)
from wafodata import WafoData
from plotbackend import plotbackend
import matplotlib
@ -177,18 +182,23 @@ class LevelCrossings(WafoData):
Z = ((u >= 0) * 2 - 1) * sqrt(-2 * log(G))
## sumcr = trapz(lc(:,1),lc(:,2))
## lc(:,2) = lc(:,2)/sumcr
## mcr = trapz(lc(:,1),lc(:,1).*lc(:,2))
## scr = trapz(lc(:,1),lc(:,1).^2.*lc(:,2))
## scr = sqrt(scr-mcr^2)
## g = lc2tr(lc,mcr,scr)
##
## f = [u u]
## f(:,2) = tranproc(Z,fliplr(g))
##
## process = tranproc(L,f)
## process = [(1:length(process)) process]
sumcr = trapz(self.data, self.args)
lc = self.data / sumcr
lc1 = self.args
mcr = trapz(lc1 * lc, lc1)
scr = trapz(lc1 ** 2 * lc, lc1)
scr = sqrt(scr - mcr ** 2)
lc2 = LevelCrossings(lc, lc1, mean=mcr, stdev=scr)
g = lc2.trdata()
f = [u, u]
f = g.dat2gauss(Z)
G = TrData(f, u)
process = G.dat2gauss(L)
return np.vstack((arange(len(process)), process)).T
##
##
## %Check the result without reference to getrfc:
@ -224,36 +234,40 @@ class LevelCrossings(WafoData):
Parameters
----------
options = structure with the fields:
csm, gsm - defines the smoothing of the crossing intensity and the
transformation g. Valid values must be
0<=csm,gsm<=1. (default csm = 0.9 gsm=0.05)
mean, sigma : real scalars
mean and standard deviation of the process
**options :
csm, gsm : real scalars
defines the smoothing of the crossing intensity and the transformation g.
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.
param :
vector which defines the region of variation of the data X.
(default [-5 5 513]).
monitor : bool
if true monitor development of estimation
linextrap : bool
if true use a smoothing spline with a constraint on the ends to
ensure linear extrapolation outside the range of the data. (default)
otherwise use a regular smoothing spline
cvar, gvar : real scalars
Variances for the crossing intensity and the empirical transformation, g. (default 1)
ne : scalar integer
Number of extremes (maxima & minima) to remove from the estimation
of the transformation. This makes the estimation more robust against
outliers. (default 7)
ntr : scalar integer
Maximum length of empirical crossing intensity. The empirical
crossing intensity is interpolated linearly before smoothing if the
length 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)
monitor monitor development of 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)
cvar - Variances for the crossing intensity. (default 1)
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.
The empirical crossing intensity is interpolated
linearly before smoothing if the length 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)
Returns
-------
gs, ge : TrData objects
smoothed and empirical estimate of the transformation g.
ma,sa = mean and standard deviation of the process
Notes
-----
@ -268,15 +282,21 @@ class LevelCrossings(WafoData):
Example
-------
Hm0 = 7
S = jonswap([],Hm0); g=ochitr([],[Hm0/4]);
S.tr = g; S.tr(:,2)=g(:,2)*Hm0/4;
xs = spec2sdat(S,2^13);
lc = dat2lc(xs)
g0 = lc2tr2(lc,0,Hm0/4,'plot','iter'); % Monitor the development
g1 = lc2tr2(lc,0,Hm0/4,troptset('gvar', .5 )); % Equal weight on all points
g2 = lc2tr2(lc,0,Hm0/4,'gvar', [3.5 .5 3.5]); % Less weight on the ends
hold on, trplot(g1,g) % Check the fit
>>> import wafo.spectrum.models as sm
>>> import wafo.transform.models as tm
>>> Hs = 7.0
>>> Sj = sm.Jonswap(Hm0=Hs)
>>> S = Sj.tospecdata() #Make spectrum object from numerical values
>>> S.tr = tm.TrOchi(mean=0, skew=sk, kurt=0, sigma=Hs/4, ysigma=Hs/4)
>>> xs = S.sim(ns=2**13)
>>> ts = mat2timeseries(xs)
>>> tp = ts.turning_points()
>>> mm = tp.cycle_pairs()
>>> lc = mm.level_crossings()
>>> g0, gemp = lc.trdata(monitor=True) # Monitor the development
>>> g1, gemp = lc.trdata(gvar=0.5 ) # Equal weight on all points
>>> g2, gemp = lc.trdata(gvar=[3.5, 0.5, 3.5]) # Less weight on the ends
hold on, trplot(g1,g) # Check the fit
trplot(g2)
See also troptset, dat2tr, trplot, findcross, smooth
@ -298,138 +318,108 @@ class LevelCrossings(WafoData):
# by pab 29.12.2000
# based on lc2tr, but the inversion is faster.
# by IR and PJ
pass
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, lin_extrap=True, ntr=1000, ne=7, cvar=1, gvar=1)
# If just 'defaults' passed in, return the default options in g
opt.update(options)
param = opt.param
Ne = opt.ne
ncr = len(self.data)
if ncr > opt.ntr and opt.ntr > 0:
x0 = linspace(self.args[Ne], self.args[-1 - Ne], opt.ntr)
lc1, lc2 = x0, interp(x0, self.args, self.data)
Ne = 0
Ner = opt.ne
ncr = opt.ntr
else:
Ner = 0
lc1, lc2 = self.args, self.data
ng = len(opt.gvar)
if ng == 1:
gvar = opt.gvar * ones(ncr)
else:
gvar = interp1d(linspace(0, 1, ng) , opt.gvar, kind='linear')( linspace(0, 1, ncr))
ng = len(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::])
else:
cor1 = 0
cor2 = 0
lc22 = 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.lin_extrap, cvar[inde])(lc11[inde])
# opt = troptset('chkder','on','plotflag','off','csm',.9,'gsm',.05,....
# 'param',[-5 5 513],'delay',2,'linextrap','on','ntr',1000,'ne',7,'cvar',1,'gvar',1);
# # If just 'defaults' passed in, return the default options in g
# if nargin==1 && nargout <= 1 && isequal(cross,'defaults')
# g = opt;
# return
# end
# error(nargchk(3,inf,nargin))
# if nargin>=4 , opt = troptset(opt,varargin{:}); end
# csm2 = opt.gsm;
# param = opt.param;
# ptime = opt.delay;
# Ne = opt.ne;
# switch opt.chkder;
# case 'off', chkder = 0;
# case 'on', chkder = 1;
# otherwise, chkder = opt.chkder;
# end
# switch opt.linextrap;
# case 'off', def = 0;
# case 'on', def = 1;
# otherwise, def = opt.linextrap;
# end
# switch opt.plotflag
# case {'none','off'}, plotflag = 0;
# case 'final', plotflag = 1;
# case 'iter', plotflag = 2;
# otherwise, plotflag = opt.plotflag;
# end
# ncr = length(cross);
# if ncr>opt.ntr && opt.ntr>0,
# x0 = linspace(cross(1+Ne,1),cross(end-Ne,1),opt.ntr)';
# cros = [ x0,interp1q(cross(:,1),cross(:,2),x0)];
# Ne = 0;
# Ner = opt.ne;
# ncr = opt.ntr;
# else
# Ner = 0;
# cros=cross;
# end
#
# ng = length(opt.gvar);
# if ng==1
# gvar = opt.gvar(ones(ncr,1));
# else
# gvar = interp1(linspace(0,1,ng)',opt.gvar(:),linspace(0,1,ncr)','*linear');
# end
# ng = length(opt.cvar);
# if ng==1
# cvar = opt.cvar(ones(ncr,1));
# else
# cvar = interp1(linspace(0,1,ng)',opt.cvar(:),linspace(0,1,ncr)','*linear');
# end
#
# g = zeros(param(3),2);
#
# uu = levels(param);
#
# g(:,1) = sa*uu' + ma;
#
# g2 = cros;
#
# if Ner>0, # Compute correction factors
# cor1 = trapz(cross(1:Ner+1,1),cross(1:Ner+1,2));
# cor2 = trapz(cross(end-Ner-1:end,1),cross(end-Ner-1:end,2));
# else
# cor1 = 0;
# cor2 = 0;
# end
# cros(:,2) = cumtrapz(cros(:,1),cros(:,2))+cor1;
# cros(:,2) = (cros(:,2)+.5)/(cros(end,2) + cor2 +1);
# cros(:,1) = (cros(:,1)-ma)/sa;
#
# # find the mode
# [tmp,imin]= min(abs(cros(:,2)-.15));
# [tmp,imax]= min(abs(cros(:,2)-.85));
# inde = imin:imax;
#tmp = smooth(cros(inde,1),g2(inde,2),opt.csm,cros(inde,1),def,cvar(inde));
#
# [tmp imax] = max(tmp);
# u0 = cros(inde(imax),1);
# #u0 = interp1q(cros(:,2),cros(:,1),.5)
#
#
# cros(:,2) = invnorm(cros(:,2),-u0,1);
#
# g2(:,2) = cros(:,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
# # to be linear outside the edges or choosing a lower value for csm2.
#
# inds = 1+Ne:ncr-Ne;# indices to points we are smoothing over
# scros2 = smooth(cros(inds,1),cros(inds,2),csm2,uu,def,gvar(inds));
#
# g(:,2) = scros2';#*sa; #multiply with stdev
#
# if chkder~=0
# for ix = 1:5
# dy = diff(g(:,2));
# if any(dy<=0)
# warning('WAFO:LCTR2','The empirical crossing spectrum is not sufficiently smoothed.')
# disp(' The estimated transfer function, g, is not ')
# disp(' a strictly increasing function.')
# dy(dy>0)=eps;
# gvar = -([dy;0]+[0;dy])/2+eps;
# g(:,2) = smooth(g(:,1),g(:,2),1,g(:,1),def,ix*gvar);
# else
# break
# end
# end
# end
# if 0, #either
# test = sqrt((param(2)-param(1))/(param(3)-1)*sum((uu-scros2).^2));
# else # or
# #test=sqrt(simpson(uu,(uu-scros2).^2));
# # or
# test=sqrt(trapz(uu,(uu-scros2).^2));
# end
#
#
# if plotflag>0,
# trplot(g ,g2,ma,sa)
# #legend(['Smoothed (T=' num2str(test) ')'],'g(u)=u','Not smoothed',0)
# #ylabel('Transfer function g(u)')
# #xlabel('Crossing level u')
#
# if plotflag>1,pause(ptime),end
# end
imax = lc222.argmax()
u0 = lc22[inde][imax]
#u0 = interp1q(cros(:,2),cros(:,1),.5)
lc22 = ndtri(lc22)-u0 #invnorm(lc22, -u0, 1);
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
# 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.lin_extrap, gvar[inds])(uu)
g = TrData(scros2, g1, mean, sigma**2) #*sa; #multiply with stdev
if opt.chkder:
for ix in range(5):
dy = diff(g.data)
if any(dy <= 0):
warnings.warn(
''' The empirical crossing spectrum is not sufficiently smoothed.
The estimated transfer function, g, is not a strictly increasing function.
''')
eps = finfo(float).eps
dy[dy > 0] = eps
gvar = -(hstack((dy, 0)) + hstack((0, dy))) / 2 + eps
g.data = SmoothSpline(g.args, g.data, 1, opt.lin_extrap, ix * gvar)(g.args)
else:
break
if opt.plotflag > 0:
g.plot()
g2.plot()
return g, g2
class CyclePairs(WafoData):
'''
@ -707,6 +697,7 @@ class TimeSeries(WafoData):
n = len(self.data)
self.args = range(0, n)
def sampling_period(self):
'''
Returns sampling interval
@ -1330,7 +1321,7 @@ def sensortypeid(*sensortypes):
>>> sensortypeid('W','v')
[11, 10]
>>> sensortypeid('rubbish')
[1.#QNAN]
[nan]
See also
--------

@ -5,10 +5,10 @@ import os
def compile_all():
# Install gfortran and run the following to build the module:
#compile_format = 'f2py %s %s -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71'
compile_format = 'f2py.py %s %s -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71'
# Install microsoft visual c++ .NET 2003 and run the following to build the module:
compile_format = 'f2py.py %s %s -c'
#compile_format = 'f2py.py %s %s -c'
pyfs = ('c_library.pyf',)
files =('c_functions.c',)

@ -16,16 +16,15 @@ 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
from wafo.misc import sub_dict_select, nextpow2, discretize, JITImport, tranproc
try:
from wafo.gaussian import Rind
except ImportError:
Rind = None
try:
from wafo import c_library
except ImportError:
warnings.warn('Compile the c_libraray.pyd again!')
warnings.warn('Compile the c_library.pyd again!')
c_library = None
from wafo.transform import TrData
@ -751,8 +750,7 @@ class SpecData1D(WafoData):
tn = paramt[1]
Ntime = paramt[2]
t = linspace(0, tn / A, Ntime) #normalized times
Nstart = max(round(t0 / tn * (Ntime - 1)), 1) #% index to starting point to
#% evaluate
Nstart = max(round(t0 / tn * (Ntime - 1)), 1) #% index to starting point to evaluate
dt = t[1] - t[0]
nr = 2
@ -1055,16 +1053,22 @@ class SpecData1D(WafoData):
if spec.tr is not None:
print(' Transforming data.')
g = spec.tr
G = fliplr(g) #% the invers of g
if derivative:
for i in range(cases):
tmp = tranproc(hstack((x[:, i + 1], xder[:, i + 1])), G)
x[:, i + 1] = tmp[:, 0]
xder[:, i + 1] = tmp[:, 1]
x[:, i + 1], xder[:, i + 1] = g.gauss2dat(x[:, i + 1], xder[:, i + 1])
else:
for i in range(cases):
x[:, i + 1] = tranproc(x[:, i + 1], G)
x[:, i + 1] = g.gauss2dat(x[:, i + 1])
# G = fliplr(g) #% the invers of g
# if derivative:
# for i in range(cases):
# tmp = tranproc(hstack((x[:, i + 1], xder[:, i + 1])), G)
# x[:, i + 1] = tmp[:, 0]
# xder[:, i + 1] = tmp[:, 1]
#
# else:
# for i in range(cases):
# x[:, i + 1] = tranproc(x[:, i + 1], G)
if derivative:
return x, xder

@ -153,12 +153,19 @@ class TrData(WafoData, TrCommon):
>>> g.sigma
array([ 5.])
>>> g = TrData(y,x,mean=1,sigma=5)
>>> g.mean
1
>>> g.sigma
5
Check that the departure from a Gaussian model is zero
>>> g.dist2gauss() < 1e-16
True
"""
def __init__(self, *args, **kwds):
WafoData.__init__(self, *args, **kwds)
super(TrData, self).__init__(*args, **kwds)
self.labels.title = 'Transform'
self.labels.ylab = 'g(x)'
self.labels.xlab = 'x'

Loading…
Cancel
Save