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

updated doctests

Browse Source
master
pbrod 9 years ago
parent b7ce91334e
commit 7b354864c4
4 changed files with 66 additions and 56 deletions
Show Stats Download Patch File Download Diff File

10
wafo/covariance/core.py
Unescape Escape View File

@ -68,7 +68,7 @@ def rndnormnd(mean, cov, cases=1):
Example Example
------- -------
>>> mu = [0, 5] >>> mu = [0, 5]
>>> S = [[1 0.45], [0.45, 0.25]] >>> S = [[1, 0.45], [0.45, 0.25]]
>>> r = rndnormnd(mu, S, 1) >>> r = rndnormnd(mu, S, 1)
plot(r(:,1),r(:,2),'.') plot(r(:,1),r(:,2),'.')
@ -211,11 +211,9 @@ class CovData1D(PlotData):
>>> abs(S1.data-S.data).max() < 1e-4 >>> abs(S1.data-S.data).max() < 1e-4
True True
>>> S1.plot('r-') S1.plot('r-')
>>> S.plot('b:') S.plot('b:')
>>> pylab.show() pylab.show()
>>> all(abs(S1.data-S.data)<1e-4)
See also See also
-------- --------

6
wafo/gaussian.py
Unescape Escape View File

@ -109,9 +109,11 @@ class Rind(object):
>>> Bup, Blo = np.atleast_2d(Bup,Blo) >>> Bup, Blo = np.atleast_2d(Bup,Blo)
>>> Bup[0,ind] = np.minimum(Bup[0,ind] , infinity*dev[indI[ind+1]]) >>> Bup[0,ind] = np.minimum(Bup[0,ind] , infinity*dev[indI[ind+1]])
>>> Blo[0,ind] = np.maximum(Blo[0,ind] ,-infinity*dev[indI[ind+1]]) >>> Blo[0,ind] = np.maximum(Blo[0,ind] ,-infinity*dev[indI[ind+1]])
>>> np.allclose(rind(Sc,m,Blo,Bup,indI, xc, nt=0), >>> val, err, terr = rind(Sc,m,Blo,Bup,indI, xc, nt=0)
... ([0.05494076], [ 0.00083066], [ 1.00000000e-10]), rtol=1e-3) >>> np.allclose(val, 0.05494076, rtol=1e-3)
True True
>>> err < 1e-3, terr< 1e-7
True, True
Compute expectation E( X1^{+}*X2^{+} ) with random Compute expectation E( X1^{+}*X2^{+} ) with random
correlation coefficient,Cov(X1,X2) = rho2. correlation coefficient,Cov(X1,X2) = rho2.

11
wafo/objects.py
Unescape Escape View File

@ -17,7 +17,7 @@ from wafo.transform.core import TrData
from wafo.transform.estimation import TransformEstimator from wafo.transform.estimation import TransformEstimator
from wafo.stats import distributions from wafo.stats import distributions
from wafo.misc import (nextpow2, findtp, findrfc, findtc, findcross, from wafo.misc import (nextpow2, findtp, findrfc, findtc, findcross,
ecross, JITImport, DotDict, gravity, findrfc_astm) ecross, JITImport, DotDict, gravity, findrfc_astm)
from wafo.interpolate import stineman_interp from wafo.interpolate import stineman_interp
from wafo.containers import PlotData from wafo.containers import PlotData
from wafo.plotbackend import plotbackend from wafo.plotbackend import plotbackend
@ -1965,10 +1965,11 @@ class TimeSeries(PlotData):
expect = 1 # reconstruct by expectation? 1=yes 0=no expect = 1 # reconstruct by expectation? 1=yes 0=no
tol = 0.001 # absolute tolerance of e(g_new-g_old) tol = 0.001 # absolute tolerance of e(g_new-g_old)
cmvmax = 100 # if number of consecutive missing values (cmv) are longer they cmvmax = 100
# are not used in estimation of g, due to the fact that the # if number of consecutive missing values (cmv) are longer they
# conditional expectation approaches zero as the length to # are not used in estimation of g, due to the fact that the
# the closest known points increases, see below in the for loop # conditional expectation approaches zero as the length to
# the closest known points increases, see below in the for loop
dT = self.sampling_period() dT = self.sampling_period()
Lm = np.minimum([n, 200, int(200/dT)]) # Lagmax 200 seconds Lm = np.minimum([n, 200, int(200/dT)]) # Lagmax 200 seconds

95
wafo/spectrum/core.py
Unescape Escape View File

@ -651,7 +651,6 @@ class SpecData1D(PlotData):
objects objects
''' '''
dt, rate = self._get_default_dt_and_rate(dt) dt, rate = self._get_default_dt_and_rate(dt)
self._check_dt(dt) self._check_dt(dt)
@ -970,8 +969,8 @@ class SpecData1D(PlotData):
Hw1 = exp(interp1d(log(abs(S1.data / S2.data)), S2.args)(freq)) Hw1 = exp(interp1d(log(abs(S1.data / S2.data)), S2.args)(freq))
else: else:
# Geometric mean # Geometric mean
Hw1 = exp((interp1d(log(abs(S1.data / S2.data)), S2.args)(freq) fun = interp1d(log(abs(S1.data / S2.data)), S2.args)
+ log(Hw2)) / 2) Hw1 = exp((fun(freq) + log(Hw2)) / 2)
# end # end
# Hw1 = (interp1q( S2.w,abs(S1.S./S2.S),freq)+Hw2)/2 # Hw1 = (interp1q( S2.w,abs(S1.S./S2.S),freq)+Hw2)/2
# plot(freq, abs(Hw11-Hw1),'g') # plot(freq, abs(Hw11-Hw1),'g')
@ -1549,7 +1548,7 @@ class SpecData1D(PlotData):
# semi-definitt, since the circulant spectrum are the eigenvalues of # semi-definitt, since the circulant spectrum are the eigenvalues of
# the circulant covariance matrix. # the circulant covariance matrix.
callFortran = 0 # callFortran = 0
# %options.method<0 # %options.method<0
# if callFortran, % call fortran # if callFortran, % call fortran
# ftmp = cov2mmtpdfexe(R,dt,u,defnr,Nstart,hg,options) # ftmp = cov2mmtpdfexe(R,dt,u,defnr,Nstart,hg,options)
@ -1732,15 +1731,17 @@ class SpecData1D(PlotData):
between Maxima, Minima and Max to level v crossing given the Max and between Maxima, Minima and Max to level v crossing given the Max and
the min is returned the min is returned
Y= Y = [Xt, Xd, Xc]
X'(t2)..X'(ts)..X'(tn-1)||X''(t1) X''(tn) X'(ts)|| X'(t1) X'(tn) X(t1) X(tn) X(ts) where
= [ Xt Xd Xc ] Xt = X'(t2)..X'(ts)..X'(tn-1)
Xd = ||X''(t1) X''(tn) X'(ts)||
Xc = X'(t1) X'(tn) X(t1) X(tn) X(ts)
Nt = tn-2, Nd = 3, Nc = 5 Nt = tn-2, Nd = 3, Nc = 5
Xt= contains Nt time points in the indicator function Xt = contains Nt time points in the indicator function
Xd= " Nd derivatives Xd = " Nd derivatives
Xc= " Nc variables to condition on Xc = " Nc variables to condition on
There are 4 (NI=5) regions with constant barriers: There are 4 (NI=5) regions with constant barriers:
(indI(1)=0); for i\in (indI(1),indI(2)] Y(i)<0. (indI(1)=0); for i\in (indI(1),indI(2)] Y(i)<0.
@ -1748,7 +1749,7 @@ class SpecData1D(PlotData):
(indI(3)=Nt+1); for i\in (indI(3)+1,indI(4)], Y(i)>0 (deriv. X''(tn)) (indI(3)=Nt+1); for i\in (indI(3)+1,indI(4)], Y(i)>0 (deriv. X''(tn))
(indI(4)=Nt+2); for i\in (indI(4)+1,indI(5)], Y(i)<0 (deriv. X'(ts)) (indI(4)=Nt+2); for i\in (indI(4)+1,indI(5)], Y(i)<0 (deriv. X'(ts))
''' '''
R0, R1, R2, R3, R4 = R[:, :5].T R0, _R1, R2, _R3, R4 = R[:, :5].T
covinput = self._covinput_mmt_pdf covinput = self._covinput_mmt_pdf
Ntime = len(R0) Ntime = len(R0)
Nx0 = max(1, len(hg)) Nx0 = max(1, len(hg))
@ -2110,12 +2111,14 @@ class SpecData1D(PlotData):
The order of the variables in the covariance matrix are organized as The order of the variables in the covariance matrix are organized as
follows: follows:
for ts <= 1: for ts <= 1:
X'(t2)..X'(ts),...,X'(tn-1) X''(t1),X''(tn) X'(t1),X'(tn),X(t1),X(tn) Xt = X'(t2)..X'(ts),...,X'(tn-1)
= [ Xt | Xd | Xc ] Xd = X''(t1), X''(tn), X'(t1), X'(tn)
Xc = X(t1),X(tn)
for ts > =2: for ts > =2:
X'(t2)..X'(ts),...,X'(tn-1) X''(t1),X''(tn) X'(ts) X'(t1),X'(tn),X(t1),X(tn) X(ts) Xt = X'(t2)..X'(ts),...,X'(tn-1)
= [ Xt | Xd | Xc ] Xd = X''(t1), X''(tn), X'(ts), X'(t1), X'(tn),
Xc = X(t1),X(tn) X(ts)
where where
@ -2123,9 +2126,11 @@ class SpecData1D(PlotData):
Xd = derivatives Xd = derivatives
Xc = variables to condition on Xc = variables to condition on
Computations of all covariances follows simple rules: Cov(X(t),X(s)) = r(t,s), Computations of all covariances follows simple rules:
Cov(X(t),X(s)) = r(t,s),
then Cov(X'(t),X(s))=dr(t,s)/dt. Now for stationary X(t) we have then Cov(X'(t),X(s))=dr(t,s)/dt. Now for stationary X(t) we have
a function r(tau) such that Cov(X(t),X(s))=r(s-t) (or r(t-s) will give the same result). a function r(tau) such that Cov(X(t),X(s))=r(s-t) (or r(t-s) will give
the same result).
Consequently Cov(X'(t),X(s)) = -r'(s-t) = -sign(s-t)*r'(|s-t|) Consequently Cov(X'(t),X(s)) = -r'(s-t) = -sign(s-t)*r'(|s-t|)
Cov(X'(t),X'(s)) = -r''(s-t) = -r''(|s-t|) Cov(X'(t),X'(s)) = -r''(s-t) = -r''(|s-t|)
@ -2141,11 +2146,11 @@ class SpecData1D(PlotData):
# for # for
i = np.arange(tn - 2) # 1:tn-2 i = np.arange(tn - 2) # 1:tn-2
# j = abs(i+1-ts) # j = abs(i+1-ts)
# BIG(i,N) = -sign(R1(j+1),R1(j+1)*dble(ts-i-1)) %cov(X'(ti+1),X(ts)) # BIG(i,N) = -sign(R1(j+1),R1(j+1)*dble(ts-i-1))
j = i + 1 - ts j = i + 1 - ts
tau = abs(j) tau = abs(j)
# BIG(i,N) = abs(R1(tau)).*sign(R1(tau).*j.') # BIG(i,N) = abs(R1(tau)).*sign(R1(tau).*j.')
BIG[i, N] = R1[tau] * sign(j) BIG[i, N] = R1[tau] * sign(j) # cov(X'(ti+1),X(ts))
# end do # end do
# Cov(Xc) # Cov(Xc)
BIG[N, N] = R0[0] # cov(X(ts),X(ts)) BIG[N, N] = R0[0] # cov(X(ts),X(ts))
@ -2154,8 +2159,8 @@ class SpecData1D(PlotData):
BIG[tn + shft + 3, N] = R0[ts] # cov(X(t1),X(ts)) BIG[tn + shft + 3, N] = R0[ts] # cov(X(t1),X(ts))
BIG[tn + shft + 4, N] = R0[tn - ts] # cov(X(tn),X(ts)) BIG[tn + shft + 4, N] = R0[tn - ts] # cov(X(tn),X(ts))
# Cov(Xd,Xc) # Cov(Xd,Xc)
BIG[tn - 1, N] = R2[ts] # %cov(X''(t1),X(ts)) BIG[tn - 1, N] = R2[ts] # cov(X''(t1),X(ts))
BIG[tn, N] = R2[tn - ts] # %cov(X''(tn),X(ts)) BIG[tn, N] = R2[tn - ts] # cov(X''(tn),X(ts))
# ADD a level u crossing at ts # ADD a level u crossing at ts
@ -2163,12 +2168,13 @@ class SpecData1D(PlotData):
# for # for
i = np.arange(tn - 2) # 1:tn-2 i = np.arange(tn - 2) # 1:tn-2
j = abs(i + 1 - ts) j = abs(i + 1 - ts)
BIG[i, tn + shft] = -R2[j] # %cov(X'(ti+1),X'(ts))
BIG[i, tn + shft] = -R2[j] # cov(X'(ti+1),X'(ts))
# end do # end do
# Cov(Xd) # Cov(Xd)
BIG[tn + shft, tn + shft] = -R2[0] # %cov(X'(ts),X'(ts)) BIG[tn + shft, tn + shft] = -R2[0] # cov(X'(ts),X'(ts))
BIG[tn - 1, tn + shft] = R3[ts] # %cov(X''(t1),X'(ts)) BIG[tn - 1, tn + shft] = R3[ts] # cov(X''(t1),X'(ts))
BIG[tn, tn + shft] = -R3[tn - ts] # %cov(X''(tn),X'(ts)) BIG[tn, tn + shft] = -R3[tn - ts] # cov(X''(tn),X'(ts))
# Cov(Xd,Xc) # Cov(Xd,Xc)
BIG[tn + shft, N] = 0.0 # %cov(X'(ts),X(ts)) BIG[tn + shft, N] = 0.0 # %cov(X'(ts),X(ts))
@ -2260,7 +2266,9 @@ class SpecData1D(PlotData):
# BIG(i,j) = BIG(j,i) # BIG(i,j) = BIG(j,i)
# end #do # end #do
# end #do # end #do
lp = np.flatnonzero(np.tril(ones(BIG.shape))) # indices to lower triangular part
# indices to lower triangular part:
lp = np.flatnonzero(np.tril(ones(BIG.shape)))
BIGT = BIG.T BIGT = BIG.T
BIG[lp] = BIGT[lp] BIG[lp] = BIGT[lp]
return BIG return BIG
@ -2304,7 +2312,7 @@ class SpecData1D(PlotData):
# compiled cov2mmtpdf.f with rind70.f # compiled cov2mmtpdf.f with rind70.f
# dos([ wafoexepath 'cov2mmtpdf.exe']) # dos([ wafoexepath 'cov2mmtpdf.exe'])
dens = 1 # load('dens.out') dens = 1 # load('dens.out')
self._cleanup(*filenames) self._cleanup(*filenames)
self._cleanup(*filenames2) self._cleanup(*filenames2)
@ -2586,7 +2594,7 @@ class SpecData1D(PlotData):
>>> import numpy as np >>> import numpy as np
>>> import scipy.stats as st >>> import scipy.stats as st
>>> x2, x1 = S.sim_nl(ns=20000,cases=20) >>> x2, x1 = S.sim_nl(ns=20000,cases=20, output='data')
>>> truth1 = [0,np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk') >>> truth1 = [0,np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk')
>>> truth1[-1] = truth1[-1]-3 >>> truth1[-1] = truth1[-1]-3
>>> np.round(truth1, 3) >>> np.round(truth1, 3)
@ -2597,7 +2605,7 @@ class SpecData1D(PlotData):
... res = fun(x2[:,1::], axis=0) ... res = fun(x2[:,1::], axis=0)
... m = res.mean() ... m = res.mean()
... sa = res.std() ... sa = res.std()
... #trueval, m, sa ... # trueval, m, sa
... np.abs(m-trueval)<sa ... np.abs(m-trueval)<sa
True True
True True
@ -2606,7 +2614,7 @@ class SpecData1D(PlotData):
>>> x = [] >>> x = []
>>> for i in range(20): >>> for i in range(20):
... x2, x1 = S.sim_nl(ns=20000,cases=1) ... x2, x1 = S.sim_nl(ns=20000,cases=1, output='data')
... x.append(x2[:,1::]) ... x.append(x2[:,1::])
>>> x2 = np.hstack(x) >>> x2 = np.hstack(x)
>>> truth1 = [0,np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk') >>> truth1 = [0,np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk')
@ -2619,7 +2627,7 @@ class SpecData1D(PlotData):
... res = fun(x2, axis=0) ... res = fun(x2, axis=0)
... m = res.mean() ... m = res.mean()
... sa = res.std() ... sa = res.std()
... #trueval, m, sa ... # trueval, m, sa
... np.abs(m-trueval)<sa ... np.abs(m-trueval)<sa
True True
True True
@ -2687,31 +2695,32 @@ class SpecData1D(PlotData):
# interpolate for freq. [1:(N/2)-1]*df and create 2-sided, uncentered # interpolate for freq. [1:(N/2)-1]*df and create 2-sided, uncentered
# spectra # spectra
f = arange(1, ns / 2.) * df ns2 = ns // 2
f_u = hstack((0., f_i, df * ns / 2.)) f = arange(1, ns2) * df
w = 2. * pi * hstack((0., f, df * ns / 2.)) f_u = hstack((0., f_i, df * ns2))
w = 2. * pi * hstack((0., f, df * ns2))
kw = w2k(w, 0., water_depth, g)[0] kw = w2k(w, 0., water_depth, g)[0]
s_u = hstack((0., abs(s_i) / 2., 0.)) s_u = hstack((0., abs(s_i) / 2., 0.))
s_i = interp(f, f_u, s_u) s_i = interp(f, f_u, s_u)
nmin = (s_i > s_max * reltol).argmax() nmin = (s_i > s_max * reltol).argmax()
nmax = flatnonzero(s_i > 0).max() nmax = flatnonzero(s_i > 0).max()
s_u = hstack((0., s_i, 0, s_i[(ns / 2) - 2::-1])) s_u = hstack((0., s_i, 0, s_i[ns2 - 2::-1]))
del(s_i, f_u) del(s_i, f_u)
# Generate standard normal random numbers for the simulations # Generate standard normal random numbers for the simulations
randn = random.randn randn = random.randn
z_r = randn((ns / 2) + 1, cases) z_r = randn(ns2 + 1, cases)
z_i = vstack((zeros((1, cases)), z_i = vstack((zeros((1, cases)),
randn((ns / 2) - 1, cases), randn(ns2 - 1, cases),
zeros((1, cases)))) zeros((1, cases))))
amp = zeros((ns, cases), dtype=complex) amp = zeros((ns, cases), dtype=complex)
amp[0:(ns / 2 + 1), :] = z_r - 1j * z_i amp[0:(ns2 + 1), :] = z_r - 1j * z_i
del(z_r, z_i) del(z_r, z_i)
amp[(ns / 2 + 1):ns, :] = amp[ns / 2 - 1:0:-1, :].conj() amp[(ns2 + 1):ns, :] = amp[ns2 - 1:0:-1, :].conj()
amp[0, :] = amp[0, :] * sqrt(2.) amp[0, :] = amp[0, :] * sqrt(2.)
amp[(ns / 2), :] = amp[(ns / 2), :] * sqrt(2.) amp[(ns2), :] = amp[(ns2), :] * sqrt(2.)
# Make simulated time series # Make simulated time series
T = (ns - 1) * d_t T = (ns - 1) * d_t
@ -2719,8 +2728,8 @@ class SpecData1D(PlotData):
if method.startswith('apd'): # apdeterministic if method.startswith('apd'): # apdeterministic
# Deterministic amplitude and phase # Deterministic amplitude and phase
amp[1:(ns / 2), :] = amp[1, 0] amp[1:(ns2), :] = amp[1, 0]
amp[(ns / 2 + 1):ns, :] = amp[1, 0].conj() amp[(ns2 + 1):ns, :] = amp[1, 0].conj()
amp = sqrt(2) * Ssqr[:, newaxis] * \ amp = sqrt(2) * Ssqr[:, newaxis] * \
exp(1J * arctan2(amp.imag, amp.real)) exp(1J * arctan2(amp.imag, amp.real))
elif method.startswith('ade'): # adeterministic elif method.startswith('ade'): # adeterministic
@ -3038,7 +3047,7 @@ class SpecData1D(PlotData):
>>> ys = wo.mat2timeseries(S.sim(ns=2**13)) >>> ys = wo.mat2timeseries(S.sim(ns=2**13))
>>> g0, gemp = ys.trdata() >>> g0, gemp = ys.trdata()
>>> t0 = g0.dist2gauss() >>> t0 = g0.dist2gauss()
>>> t1 = S0.testgaussian(ns=2**13, t0=t0, cases=50) >>> t1 = S0.testgaussian(ns=2**13, cases=50)
>>> sum(t1 > t0) < 5 >>> sum(t1 > t0) < 5
True True

Write Preview
Loading…
Cancel
Save