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

Fixed a bug and added more tests

Browse Source
master
pbrod 9 years ago
parent d6a5d0e807
commit e17d08fa61
3 changed files with 230 additions and 169 deletions
Show Stats Download Patch File Download Diff File

6
wafo/misc.py
Unescape Escape View File

@ -962,7 +962,7 @@ def findcross(x, v=0.0, kind=None, method='clib'):
if xor(is_odd, kind in ('dw', 'uw')):
ind = ind[:-1]
else:
raise ValueError('Unknown wave/crossing definition!')
raise ValueError('Unknown wave/crossing definition! {}'.format(kind))
return ind
@ -1973,7 +1973,7 @@ def findtc(x_in, v=None, kind=None):
if (2 * n_tc + 1 < n_c) and (kind in (None, 'tw')):
# trough
ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]].argmin()
ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]+1].argmin()
else: # the first is a up-crossing
for i in range(n_tc):
@ -1985,7 +1985,7 @@ def findtc(x_in, v=None, kind=None):
if (2 * n_tc + 1 < n_c) and (kind in (None, 'cw')):
# crest
ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]].argmax()
ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]+1].argmax()
return v_ind[:n_c - 1] + ind + 1, v_ind

262
wafo/objects.py
Unescape Escape View File

@ -21,7 +21,7 @@ from wafo.misc import (nextpow2, findtp, findrfc, findtc, findcross,
detrendma)
from wafo.interpolate import stineman_interp
from wafo.containers import PlotData
from wafo.plotbackend import plotbackend
from wafo.plotbackend import plotbackend as plt
from scipy.integrate import trapz
from scipy.signal import welch, lfilter
from scipy.signal.windows import get_window # @UnusedImport
@ -69,9 +69,9 @@ class LevelCrossings(PlotData):
Examples
--------
>>> import wafo.data
>>> import wafo.data as wd
>>> import wafo.objects as wo
>>> x = wafo.data.sea()
>>> x = wd.sea()
>>> ts = wo.mat2timeseries(x)
>>> tp = ts.turning_points()
@ -80,8 +80,14 @@ class LevelCrossings(PlotData):
>>> lc = mm.level_crossings()
>>> np.allclose(lc.data[:5], [ 0., 1., 2., 2., 3.])
True
>>> m, s = lc.estimate_mean_and_stdev()
>>> np.allclose([m, s], (0.033974280952584639, 0.48177752818956326))
True
>>> np.allclose((lc.mean, lc.sigma),
... (1.5440875692709283e-09, 0.47295493383306714))
True
h2 = lc.plot()
>>> h2 = lc.plot()
'''
def __init__(self, *args, **kwds):
@ -97,27 +103,37 @@ class LevelCrossings(PlotData):
self.mean = kwds.get('mean')
# self.setplotter(plotmethod='step')
icmax = self.data.argmax()
if self.data is not None:
i_cmax = self.data.argmax()
if self.sigma is None or self.mean is None:
logcros = where(self.data == 0.0, inf, -log(self.data))
logcmin = logcros[icmax]
logcros = sqrt(2 * abs(logcros - logcmin))
logcros[0:icmax + 1] = 2 * logcros[
icmax] - logcros[0:icmax + 1]
ncr = 10
# least square fit
p = polyfit(self.args[ncr:-ncr], logcros[ncr:-ncr], 1)
mean, sigma = self.estimate_mean_and_stdev(i_cmax)
if self.sigma is None:
# estimated standard deviation of x
self.sigma = 1.0 / p[0]
self.sigma = sigma
if self.mean is None:
self.mean = -p[1] / p[0] # self.args[icmax]
cmax = self.data[icmax]
self.mean = mean
cmax = self.data[i_cmax]
x = (self.args - self.mean) / self.sigma
y = cmax * exp(-x ** 2 / 2.0)
self.children = [PlotData(y, self.args)]
def estimate_mean_and_stdev(self, i_cmax=None):
"""
Return mean and standard deviation of process x estimated from crossing
"""
if i_cmax is None:
i_cmax = self.data.argmax()
logcros = where(self.data == 0.0, inf, -log(self.data))
logcmin = logcros[i_cmax]
logcros = sqrt(2 * abs(logcros - logcmin))
logcros[0:i_cmax + 1] = 2 * logcros[i_cmax] - logcros[0:i_cmax + 1]
ncr = 10
# least square fit
p = polyfit(self.args[ncr:-ncr], logcros[ncr:-ncr], 1)
sigma = 1.0 / p[0]
mean = -p[1] / p[0] # self.args[i_cmax]
return mean, sigma
def extrapolate(self, u_min=None, u_max=None, method='ml', dist='genpar',
plotflag=0):
'''
@ -137,8 +153,8 @@ class LevelCrossings(PlotData):
expon : Exponential distribution (GPD with k=0)
rayleigh : truncated Rayleigh distribution
plotflag : scalar integer
1: Diagnostic plots. (default)
0: Don't plot diagnostic plots.
1: Diagnostic plots.
0: Don't plot diagnostic plots. (default)
Returns
-------
@ -155,7 +171,7 @@ class LevelCrossings(PlotData):
H(x) = exp(-x/s), k=0 (expon)
The tails with the survival function of a truncated Rayleigh
distribution.
H(x) = exp(-((x+x0).^2-x0^2)/s^2) (rayleigh)
H(x) = exp(-((x+x0)**2-x0^2)/s**2) (rayleigh)
where x0 is the distance from the truncation level to where the LC has
its maximum.
The method 'gpd' uses the GPD. We recommend the use of 'gpd,ml'.
@ -165,24 +181,36 @@ class LevelCrossings(PlotData):
Example
-------
>>> import wafo.data
>>> import wafo.data as wd
>>> import wafo.objects as wo
>>> x = wafo.data.sea()
>>> x = wd.sea()
>>> ts = wo.mat2timeseries(x)
>>> tp = ts.turning_points()
>>> mm = tp.cycle_pairs()
>>> lc = mm.level_crossings()
>>> s = x[:,1].std()
>>> s = x[:, 1].std()
>>> lc_gpd = lc.extrapolate(-2*s, 2*s)
>>> lc_exp = lc.extrapolate(-2*s, 2*s, dist='expon')
>>> lc_ray = lc.extrapolate(-2*s, 2*s, dist='rayleigh')
lc.plot()
lc_gpd.plot()
lc_exp.plot()
lc_ray.plot()
>>> n = 3
>>> np.allclose([lc_gpd.data[:n], lc_gpd.data[-n:]],
... [[ 0., 0., 0.], [ 0., 0., 0.]])
True
>>> np.allclose([lc_exp.data[:n], lc_exp.data[-n:]],
... [[ 6.51864195e-12, 7.02339889e-12, 7.56724060e-12],
... [ 1.01040335e-05, 9.70417448e-06, 9.32013956e-06]])
True
>>> np.allclose([lc_ray.data[:n], lc_ray.data[-n:]],
... [[ 1.78925398e-37, 2.61098785e-37, 3.80712964e-37],
... [ 1.28140956e-13, 1.11668143e-13, 9.72878135e-14]])
True
>>> h0 = lc.plot()
>>> h1 = lc_gpd.plot()
>>> h2 = lc_exp.plot()
>>> h3 = lc_ray.plot()
See also
@ -209,17 +237,16 @@ class LevelCrossings(PlotData):
u_max = self.args[i.max()]
lcf, lcx = self.data, self.args
# Extrapolate LC for high levels
[lc_High, phat_high] = self._extrapolate(lcx, lcf, u_max,
u_max - lc_max, method, dist)
lc_High, phat_high = self._extrapolate(lcx, lcf, u_max, u_max - lc_max,
method, dist)
# Extrapolate LC for low levels
[lcEst1, phat_low] = self._extrapolate(-lcx[::-1], lcf[::-1], -u_min,
lcEst1, phat_low = self._extrapolate(-lcx[::-1], lcf[::-1], -u_min,
lc_max - u_min, method, dist)
lc_Low = lcEst1[::-1, :] # [-lcEst1[::-1, 0], lcEst1[::-1, 1::]]
lc_Low[:, 0] *= -1
if plotflag:
plotbackend.semilogx(lcf, lcx,
lc_High[:, 1], lc_High[:, 0],
plt.semilogx(lcf, lcx, lc_High[:, 1], lc_High[:, 0],
lc_Low[:, 1], lc_Low[:, 0])
i_mask = (u_min < lcx) & (lcx < u_max)
f = np.hstack((lc_Low[:, 1], lcf[i_mask], lc_High[:, 1]))
@ -267,7 +294,7 @@ class LevelCrossings(PlotData):
r = sqrt(-2 * log(1 - 90 / 100)) # 90 # confidence sphere
Nc = 16 + 1
ang = linspace(0, 2 * pi, Nc)
# 90# Circle
# 90% Circle
c0 = np.vstack(
(r * sqrt(D[0]) * sin(ang), r * sqrt(D[1]) * cos(ang)))
# plot(c0(1,:),c0(2,:))
@ -309,7 +336,7 @@ class LevelCrossings(PlotData):
# F = -np.expm1(-((xF + offset) ** 2 - offset ** 2) / s ** 2)
lcEst = np.vstack((xF + u, lcu * (SF))).T
else:
raise ValueError()
raise NotImplementedError('Unknown distribution {}'.format(dist))
return lcEst, phat
# End extrapolate
@ -373,14 +400,14 @@ class LevelCrossings(PlotData):
>>> lc2 = ts2.turning_points().cycle_pairs().level_crossings()
import pylab as plt
h0 = S.plot('b')
h1 = Se.plot('r')
>>> import pylab as plt
>>> h0 = S.plot('b')
>>> h1 = Se.plot('r')
h = plt.subplot(211)
h2 = lc2.plot()
h = plt.subplot(212)
h0 = lc.plot()
>>> h = plt.subplot(211)
>>> h2 = lc2.plot()
>>> h = plt.subplot(212)
>>> h0 = lc.plot()
"""
@ -438,8 +465,7 @@ class LevelCrossings(PlotData):
lc1 = self.args
mcr = trapz(lc1 * lc, lc1) if self.mean is None else self.mean
if self.sigma is None:
scr = trapz(lc1 ** 2 * lc, lc1)
scr = sqrt(scr - mcr ** 2)
scr = sqrt(trapz(lc1 ** 2 * lc, lc1) - mcr ** 2)
else:
scr = self.sigma
lc2 = LevelCrossings(lc, lc1, mean=mcr, sigma=scr, intensity=True)
@ -582,19 +608,6 @@ class LevelCrossings(PlotData):
return estimate._trdata_lc(self, mean, sigma)
def test_levelcrossings_extrapolate():
import wafo.data
x = wafo.data.sea()
ts = mat2timeseries(x)
tp = ts.turning_points()
mm = tp.cycle_pairs()
lc = mm.level_crossings()
s = x[:, 1].std()
lc_gpd = lc.extrapolate(-2 * s, 2 * s, dist='rayleigh') # @UnusedVariable
class CyclePairs(PlotData):
'''
Container class for Cycle Pairs data objects in WAFO
@ -612,12 +625,22 @@ class CyclePairs(PlotData):
>>> ts = wo.mat2timeseries(x)
>>> tp = ts.turning_points()
>>> mm = tp.cycle_pairs()
>>> np.allclose(mm.data[:5],
>>> mM = tp.cycle_pairs(kind='min2max')
>>> np.allclose(mM.data[:5],
... [ 0.83950546, -0.02049454, -0.04049454, 0.25950546, -0.08049454])
True
>>> np.allclose(mM.args[:5],
... [-1.2004945 , -0.09049454, -0.09049454, -0.16049454, -0.43049454])
True
>>> Mm = tp.cycle_pairs(kind='max2min')
>>> np.allclose(Mm.data[:5],
... [ 0.83950546, -0.02049454, -0.04049454, 0.25950546, -0.08049454])
True
>>> np.allclose(Mm.args[:5],
... [-0.09049454, -0.09049454, -0.16049454, -0.43049454, -0.21049454])
True
h1 = mm.plot(marker='x')
>>> h1 = mM.plot(marker='x')
'''
def __init__(self, *args, **kwds):
@ -666,12 +689,11 @@ class CyclePairs(PlotData):
>>> bv = range(3,9)
>>> D = mm.damage(beta=bv)
>>> D
array([ 138.5238799 , 117.56050788, 108.99265423, 107.86681126,
112.3791076 , 122.08375071])
h = mm.plot(marker='.')
h = plt.plot(bv,D,'x-')
>>> np.allclose(D, [ 138.5238799 , 117.56050788, 108.99265423,
... 107.86681126, 112.3791076 , 122.08375071])
True
>>> h = mm.plot(marker='.')
>>> h = plt.plot(bv, D, 'x-')
See also
--------
@ -730,14 +752,10 @@ class CyclePairs(PlotData):
LevelCrossings
"""
if isinstance(kind, str):
t = dict(u=0, uM=1, umM=2, um=3)
defnr = t.get(kind, 1)
else:
defnr = kind
defnr = dict(u=0, uM=1, umM=2, um=3).get(kind, kind)
if defnr not in [1, 2, 3, 4]:
raise ValueError('kind must be one of (1,2,3,4).')
raise ValueError('kind must be one of (1, 2, 3, 4, "u", "uM",'
' "umM", "um"). Got kind = {}'.format(kind))
m, M = self.get_minima_and_maxima()
@ -847,11 +865,13 @@ class TurningPoints(PlotData):
>>> np.allclose(tph.data[:5],
... [-0.16049454, 0.25950546, -0.43049454, -0.08049454, -0.42049454])
True
>>> np.allclose(tph.args[:5],
... [ 7.05, 7.8 , 9.8 , 11.8 , 12.8 ])
True
hs = ts1.plot()
hp = tp.plot('ro')
hph = tph.plot('k.')
>>> hs = ts1.plot()
>>> hp = tp.plot('ro')
>>> hph = tph.plot('k.')
See also
---------
@ -896,7 +916,7 @@ class TurningPoints(PlotData):
... 0.25950546, -0.08049454])
True
h = mM.plot(marker='x')
>>> h = mM.plot(marker='x')
See also
@ -1205,9 +1225,15 @@ class TimeSeries(PlotData):
Example
-------
x = load('sea.dat');
S = dat2spec(x);
specplot(S)
>>> import wafo.data as wd
>>> import wafo.objects as wo
>>> x = wd.sea()
>>> ts = wo.mat2timeseries(x)
>>> S = ts.tospecdata()
>>> np.allclose(S.data[21:25],
... (0.00207876, 0.0025113 , 0.00300008, 0.00351852))
True
>>> h = S.plot()
See also
--------
@ -1717,7 +1743,7 @@ class TimeSeries(PlotData):
# sorting crossings and tp in sequence
index = sort(r_[index, v_ind])
else:
raise ValueError('Unknown pdef option!')
raise ValueError('Unknown pdef option! {}'.format(str(pdef)))
return index
def _get_start_index(self, pdef, down_crossing_or_max):
@ -2221,7 +2247,7 @@ class TimeSeries(PlotData):
>>> x = wafo.data.sea()
>>> ts150 = wafo.objects.mat2timeseries(x[:150,:])
h = ts150.plot_wave('r-', sym2='bo')
>>> h = ts150.plot_wave('r-', sym2='bo')
See also
--------
@ -2272,12 +2298,12 @@ class TimeSeries(PlotData):
if np.max(abs(xn[indg])) > 5 * sigma:
XlblTxt = XlblTxt + ' (Spurious data since max > 5 std.)'
plot = plotbackend.plot
subplot = plotbackend.subplot
plot = plt.plot
subplot = plt.subplot
figs = []
for unused_iz in range(nfig):
figs.append(plotbackend.figure())
plotbackend.title('Surface elevation from mean water level (MWL).')
figs.append(plt.figure())
plt.title('Surface elevation from mean water level (MWL).')
for ix in range(nsub):
if nsub > 1:
subplot(nsub, 1, ix)
@ -2287,11 +2313,11 @@ class TimeSeries(PlotData):
if len(ind2) > 0:
plot(tn2[ind2] * dT, xn2[ind2], sym2)
plot(h_scale * dT, [0, 0], 'k-')
# plotbackend.axis([h_scale*dT, v_scale])
# plt.axis([h_scale*dT, v_scale])
for iy in [-2, 2]:
plot(h_scale * dT, iy * sigma * ones(2), ':')
ind = ind + Ns
plotbackend.xlabel(XlblTxt)
plt.xlabel(XlblTxt)
return figs
@ -2315,7 +2341,7 @@ class TimeSeries(PlotData):
>>> x = wafo.data.sea()
>>> ts = wafo.objects.mat2timeseries(x[0:500,...])
h = ts.plot_sp_wave(np.r_[6:9,12:18])
>>> h = ts.plot_sp_wave(np.r_[6:9,12:18])
See also
--------
@ -2350,75 +2376,33 @@ class TimeSeries(PlotData):
Nsub = min(6, int(ceil(Nsub / Nfig)))
figs = []
for unused_iy in range(Nfig):
figs.append(plotbackend.figure())
figs.append(plt.figure())
for ix in range(Nsub):
plotbackend.subplot(Nsub, 1, mod(ix, Nsub) + 1)
plt.subplot(Nsub, 1, mod(ix, Nsub) + 1)
ind = r_[tz_idx[2 * wave_idx[ix] - 1]:tz_idx[
2 * wave_idx[ix] + 2 * Nwp[ix] - 1]]
# indices to wave
plotbackend.plot(self.args[ind], self.data[ind], *args, **kwds)
plotbackend.hold('on')
plt.plot(self.args[ind], self.data[ind], *args, **kwds)
plt.hold('on')
xi = [self.args[ind[0]], self.args[ind[-1]]]
plotbackend.plot(xi, [0, 0])
plt.plot(xi, [0, 0])
if Nwp[ix] == 1:
plotbackend.ylabel('Wave %d' % wave_idx[ix])
plt.ylabel('Wave %d' % wave_idx[ix])
else:
plotbackend.ylabel(
plt.ylabel(
'Wave %d - %d' % (wave_idx[ix],
wave_idx[ix] + Nwp[ix] - 1))
plotbackend.xlabel('Time [sec]')
plt.xlabel('Time [sec]')
# wafostamp
return figs
def main():
import wafo
ts = wafo.objects.mat2timeseries(wafo.data.sea())
_S = ts.tospecdata(method='psd')
tp = ts.turning_points()
mm = tp.cycle_pairs()
lc = mm.level_crossings()
lc.plot()
T = ts.wave_periods(vh=0.0, pdef='c2c') # @UnusedVariable
# main()
import wafo.spectrum.models as sm
Sj = sm.Jonswap()
S = Sj.tospecdata()
R = S.tocovdata()
x = R.sim(ns=1000, dt=0.2) # @UnusedVariable
S.characteristic(['hm0', 'tm02'])
ns = 1000
dt = .2
x1 = S.sim(ns, dt=dt)
ts = TimeSeries(x1[:, 1], x1[:, 0])
tp = ts.turning_points(0.0)
x = np.arange(-2, 2, 0.2)
# Plot 2 objects in one call
d2 = PlotData(np.sin(x), x, xlab='x', ylab='sin', title='sinus')
d0 = d2.copy()
d0.data = d0.data * 0.9
d1 = d2.copy()
d1.data = d1.data * 1.2
d1.children = [d0]
d2.children = [d1]
d2.plot()
print('Done')
def test_docstrings():
import doctest
print('Testing docstrings in %s' % __file__)
doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
if __name__ == '__main__':
test_docstrings()
# main()
# test_levelcrossings_extrapolate()

127
wafo/tests/test_objects.py
Unescape Escape View File

@ -14,11 +14,13 @@ import wafo.spectrum.models as sm
import wafo.transform.models as tm
class TestObjects(TestCase):
def test_timeseries(self):
class TestTimeSeries(TestCase):
def setUp(self):
x = wafo.data.sea()
ts = wo.mat2timeseries(x)
self.ts = wo.mat2timeseries(x)
def test_timeseries(self):
ts = self.ts
assert_array_almost_equal(ts.sampling_period(), 0.25)
S = ts.tospecdata()
@ -54,12 +56,57 @@ class TestObjects(TestCase):
self.assert_(0.54 < g1.dist2gauss() < 0.95)
self.assert_(1.5 < g2.dist2gauss() < 1.9)
def test_cycles_and_levelcrossings(self):
def test_timeseries_wave_periods(self):
true_t = ([-0.69, -0.86, -1.05],
[0.42, 0.78, 1.37],
[0.09, 0.51, -0.85],
[-0.27, -0.08, 0.32],
[3.84377468, 3.60707656, 3.90490909],
[6.25273295, 6.35295202, 6.11978685],
[2.48364668, 2.49282402, 2.50553431],
[3.76908628, 3.860128, 3.61425254],
[-5.05027968, -5.41405436, -5.35113091],
[7.53392635, 7.90687837, 7.85666522],
[-0.2811934, -0.86392635, -0.87687837],
[4.05027968, 4.72405436, 4.49113092],
[2.03999996, 0.07, 0.05],
[-0.93, -0.07, -0.12],
[1.10999996, 0., -0.07],
[-0.86, -0.02, 0.3],
[0.93, -0.8, -0.2])
pdefs = ['t2c', 'c2t', 't2t', 'c2c',
'd2d', 'u2u', 'd2u', 'u2d',
'd2t', 't2u', 'u2c', 'c2d',
'm2M', 'M2m', 'm2m', 'M2M', 'all',
]
ts = wo.TimeSeries(self.ts.data[0:400, :], self.ts.args[:400])
for pdef, truth in zip(pdefs, true_t):
T, _ix = ts.wave_periods(vh=0.0, pdef=pdef)
# print(T[:3, 0])
assert_array_almost_equal(T[:3, 0], truth)
true_t2 = ([1.10999996, 0., - 0.07],
[-0.02, 0.3, - 0.34],
[6.10295202, 5.86978685, 6.08501107],
[6.25273295, 6.35295202, 6.11978685],
[-0.27, -0.08, 0.32],
[-0.27, -0.08, 0.32])
wdefs = ['mw', 'Mw', 'dw', 'uw', 'tw', 'cw', ]
for wdef, truth in zip(wdefs, true_t2):
pdef = '{0}2{0}'.format(wdef[0].lower())
T, _ix = ts.wave_periods(vh=0.0, pdef=pdef, wdef=wdef)
# print(T[:3, 0])
assert_array_almost_equal(T[:3, 0], truth)
class TestObjects(TestCase):
def setUp(self):
x = wafo.data.sea()
ts = wo.mat2timeseries(x)
self.ts = wo.mat2timeseries(x)
tp = ts.turning_points()
def test_cycles_and_levelcrossings(self):
tp = self.ts.turning_points()
assert_array_almost_equal(tp.data[:10],
[-1.200495, 0.839505, -0.090495, -0.020495,
-0.090495, -0.040495, -0.160495, 0.259505,
@ -71,29 +118,59 @@ class TestObjects(TestCase):
[0.839505, -0.020495, -0.040495, 0.259505,
-0.080495, -0.080495, 0.349505, 0.859505,
0.009505, 0.319505])
lc = mm.level_crossings()
assert_array_almost_equal(lc.data[:10],
[0., 1., 2., 2., 3., 4.,
5., 6., 7., 9.])
true_lcs = (([0., 1., 2., 2., 3., 4., 5., 6., 7., 9.],
[-1.7504945, -1.4404945, -1.4204945, -1.4004945,
-1.3704945, -1.3204945, -1.2704945, -1.2604945,
-1.2504945, -1.2004945]),
([0., 1., 2., 3., 3., 4., 5., 6., 7., 9.],
[-1.7504945, -1.4404945, -1.4204945, -1.4004945,
-1.3704945, -1.3204945, -1.2704945, -1.2604945,
-1.2504945, -1.2004945]),
([1., 2., 3., 4., 4., 5., 6., 7., 9., 11.],
[-1.7504945, -1.4404945, -1.4204945, -1.4004945,
-1.3704945, -1.3204945, -1.2704945, -1.2604945,
-1.2504945, -1.2004945]),
([1., 2., 3., 3., 4., 5., 6., 7., 9., 11.],
[-1.7504945, -1.4404945, -1.4204945, -1.4004945,
-1.3704945, -1.3204945, -1.2704945, -1.2604945,
-1.2504945, -1.2004945]))
for i, true_lc in enumerate(true_lcs):
true_count, true_levels = true_lc
lc = mm.level_crossings(kind=i+1)
assert_array_almost_equal(lc.data[:10], true_count)
assert_array_almost_equal(lc.args[:10], true_levels)
def test_levelcrossings_extrapolate(self):
x = wafo.data.sea()
ts = wo.mat2timeseries(x)
tp = ts.turning_points()
tp = self.ts.turning_points()
mm = tp.cycle_pairs()
lc = mm.level_crossings()
s = x[:, 1].std()
lc_gpd = lc.extrapolate(-2 * s, 2 * s, dist='rayleigh')
assert_array_almost_equal(lc_gpd.data[:10],
[1.789254e-37, 2.610988e-37,
3.807130e-37, 5.546901e-37,
8.075384e-37, 1.174724e-36,
1.707531e-36, 2.480054e-36,
3.599263e-36, 5.219466e-36])
s = lc.sigma # x[:, 1].std()
ix = slice(0, 1000, 100)
lc_ray = lc.extrapolate(-2 * s, 2 * s, dist='rayleigh')
assert_array_almost_equal(lc_ray.data[ix],
[1.78925398e-37, 9.61028192e-23,
2.05282628e-11, 1.74389448e-03,
5.89169345e+01, 5.240000e+02,
6.72609651e+01, 4.46086175e-01,
2.23463577e-04, 8.45526153e-09])
lc_exp = lc.extrapolate(-2 * s, 2 * s, dist='expon')
lc_gpd = lc.extrapolate(-2 * s, 2 * s, dist='genpareto')
assert_array_almost_equal(lc_exp.data[ix],
[6.51864195e-12, 1.13025876e-08,
1.95974080e-05, 3.39796881e-02,
5.89169345e+01, 5.24000000e+02,
6.43476951e+01, 1.13478843e+00,
2.00122906e-02, 3.52921977e-04])
assert_array_almost_equal(lc_gpd.data[ix],
[0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00,
5.89169345e+01, 5.24000000e+02,
6.80484770e+01, 1.41019390e-01,
0.00000000e+00, 0.00000000e+00])
if __name__ == "__main__":

Write Preview
Loading…
Cancel
Save