You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
2335 lines
82 KiB
Python
2335 lines
82 KiB
Python
|
|
|
|
# Name: module1
|
|
# Purpose:
|
|
#
|
|
# Author: pab
|
|
#
|
|
# Created: 16.09.2008
|
|
# Copyright: (c) pab 2008
|
|
# Licence: <your licence>
|
|
|
|
# !/usr/bin/env python
|
|
|
|
|
|
from __future__ import absolute_import, division
|
|
from wafo.transform.core import TrData
|
|
from wafo.transform.estimation import TransformEstimator
|
|
from wafo.stats import distributions
|
|
from wafo.misc import (nextpow2, findtp, findrfc, findtc, findcross,
|
|
ecross, JITImport, DotDict, gravity, findrfc_astm,
|
|
detrendma)
|
|
from wafo.interpolate import stineman_interp
|
|
from wafo.containers import PlotData
|
|
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
|
|
from scipy import special
|
|
from scipy.interpolate.interpolate import interp1d
|
|
from scipy.special import ndtr as cdfnorm
|
|
import warnings
|
|
import numpy as np
|
|
from numpy import (inf, pi, zeros, ones, sqrt, where, log, exp, cos, sin,
|
|
arcsin, mod,
|
|
linspace, arange, sort, all, abs, vstack, hstack,
|
|
atleast_1d, finfo, polyfit, r_, nonzero,
|
|
cumsum, ravel, isnan, ceil, diff, array)
|
|
from numpy.fft import fft # @UnusedImport
|
|
from numpy.random import randn
|
|
from matplotlib.mlab import psd, detrend_mean
|
|
from scipy.signal.windows import parzen
|
|
|
|
|
|
floatinfo = finfo(float)
|
|
|
|
_wafocov = JITImport('wafo.covariance')
|
|
_wafocov_estimation = JITImport('wafo.covariance.estimation')
|
|
_wafospec = JITImport('wafo.spectrum')
|
|
|
|
__all__ = ['TimeSeries', 'LevelCrossings', 'CyclePairs', 'TurningPoints',
|
|
'sensortypeid', 'sensortype']
|
|
|
|
|
|
def _invchi2(q, df):
|
|
return special.chdtri(df, q)
|
|
|
|
|
|
class LevelCrossings(PlotData):
|
|
|
|
'''
|
|
Container class for Level crossing data objects in WAFO
|
|
|
|
Member variables
|
|
----------------
|
|
data : array-like
|
|
number of upcrossings or upcrossingintensity
|
|
args : array-like
|
|
crossing levels
|
|
|
|
Examples
|
|
--------
|
|
>>> import wafo.data as wd
|
|
>>> import wafo.objects as wo
|
|
>>> x = wd.sea()
|
|
>>> ts = wo.mat2timeseries(x)
|
|
|
|
>>> tp = ts.turning_points()
|
|
>>> mm = tp.cycle_pairs()
|
|
|
|
>>> 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()
|
|
'''
|
|
|
|
def __init__(self, *args, **kwds):
|
|
options = dict(title='Level crossing spectrum',
|
|
xlab='Levels', ylab='Count',
|
|
plotmethod='semilogy',
|
|
plot_args=['b'],
|
|
plot_args_children=['r--'])
|
|
options.update(**kwds)
|
|
super(LevelCrossings, self).__init__(*args, **options)
|
|
self.intensity = kwds.get('intensity', False)
|
|
self.sigma = kwds.get('sigma')
|
|
self.mean = kwds.get('mean')
|
|
# self.setplotter(plotmethod='step')
|
|
|
|
if self.data is not None:
|
|
i_cmax = self.data.argmax()
|
|
if self.sigma is None or self.mean is None:
|
|
mean, sigma = self.estimate_mean_and_stdev(i_cmax)
|
|
if self.sigma is None:
|
|
# estimated standard deviation of x
|
|
self.sigma = sigma
|
|
if self.mean is None:
|
|
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):
|
|
'''
|
|
Returns an extrapolated level crossing spectrum
|
|
|
|
Parameters
|
|
-----------
|
|
u_min, u_max : real scalars
|
|
extrapolate below u_min and above u_max.
|
|
method : string
|
|
describing the method of estimation. Options are:
|
|
'ml' : Maximum Likelihood method (default)
|
|
'mps': Maximum Product Spacing method
|
|
dist : string
|
|
defining distribution function. Options are:
|
|
genpareto : Generalized Pareto distribution (GPD)
|
|
expon : Exponential distribution (GPD with k=0)
|
|
rayleigh : truncated Rayleigh distribution
|
|
plotflag : scalar integer
|
|
1: Diagnostic plots.
|
|
0: Don't plot diagnostic plots. (default)
|
|
|
|
Returns
|
|
-------
|
|
lc : LevelCrossing object
|
|
with the estimated level crossing spectrum
|
|
Est = Estimated parameters. [struct array]
|
|
|
|
Extrapolates the level crossing spectrum (LC) for high and for low
|
|
levels.
|
|
The tails of the LC is fitted to a survival function of a GPD.
|
|
H(x) = (1-k*x/s)^(1/k) (GPD)
|
|
The use of GPD is motivated by POT methods in extreme value theory.
|
|
For k=0 the GPD is the exponential distribution
|
|
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)
|
|
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'.
|
|
The method 'exp' uses the Exp.
|
|
The method 'ray' uses Ray, and should be used if the load is a
|
|
Gaussian process.
|
|
|
|
Example
|
|
-------
|
|
>>> import wafo.data as wd
|
|
>>> import wafo.objects as wo
|
|
>>> x = wd.sea()
|
|
>>> ts = wo.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)
|
|
>>> lc_exp = lc.extrapolate(-2*s, 2*s, dist='expon')
|
|
>>> lc_ray = lc.extrapolate(-2*s, 2*s, dist='rayleigh')
|
|
|
|
>>> 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
|
|
--------
|
|
cmat2extralc, rfmextrapolate, lc2rfmextreme, extralc, fitgenpar
|
|
|
|
References
|
|
----------
|
|
Johannesson, P., and Thomas, J-.J. (2000):
|
|
Extrapolation of Rainflow Matrices.
|
|
Preprint 2000:82, Mathematical statistics, Chalmers, pp. 18.
|
|
'''
|
|
|
|
i_max = self.data.argmax()
|
|
c_max = self.data[i_max]
|
|
lc_max = self.args[i_max]
|
|
|
|
if u_min is None or u_max is None:
|
|
fraction = sqrt(c_max)
|
|
i = np.flatnonzero(self.data > fraction)
|
|
if u_min is None:
|
|
u_min = self.args[i.min()]
|
|
if u_max is None:
|
|
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)
|
|
# Extrapolate LC for low levels
|
|
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:
|
|
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]))
|
|
x = np.hstack((lc_Low[:, 0], lcx[i_mask], lc_High[:, 0]))
|
|
lc_out = LevelCrossings(f, x, sigma=self.sigma, mean=self.mean)
|
|
lc_out.phat_high = phat_high
|
|
lc_out.phat_low = phat_low
|
|
return lc_out
|
|
|
|
def _extrapolate(self, lcx, lcf, u, offset, method, dist):
|
|
# Extrapolate the level crossing spectra for high levels
|
|
|
|
method = method.lower()
|
|
dist = dist.lower()
|
|
|
|
# Excedences over level u
|
|
Iu = lcx > u
|
|
lcx1, lcf1 = lcx[Iu], lcf[Iu]
|
|
lcf2, lcx2 = self._make_increasing(lcf1[::-1], lcx1[::-1])
|
|
|
|
nim1 = 0
|
|
x = []
|
|
for xk, ni in zip(lcx2.tolist(), lcf2.tolist()):
|
|
ni = int(ni)
|
|
x.append(ones(ni - nim1) * xk)
|
|
nim1 = ni
|
|
|
|
x = np.hstack(x) - u
|
|
|
|
df = 0.01
|
|
xF = np.arange(0.0, 4 + df / 2, df)
|
|
lcu = np.interp(u, lcx, lcf) + 1
|
|
# Estimate tail
|
|
if dist.startswith('gen'):
|
|
genpareto = distributions.genpareto
|
|
phat = genpareto.fit2(x, floc=0, method=method)
|
|
SF = phat.sf(xF)
|
|
|
|
covar = phat.par_cov[::2, ::2]
|
|
# Calculate 90 # confidence region, an ellipse, for (k,s)
|
|
D, B = np.linalg.eig(covar)
|
|
b = phat.par[::2]
|
|
if b[0] > 0:
|
|
phat.upperlimit = u + b[1] / b[0]
|
|
|
|
r = sqrt(-2 * log(1 - 90 / 100)) # 90 # confidence sphere
|
|
Nc = 16 + 1
|
|
ang = linspace(0, 2 * pi, Nc)
|
|
# 90% Circle
|
|
c0 = np.vstack(
|
|
(r * sqrt(D[0]) * sin(ang), r * sqrt(D[1]) * cos(ang)))
|
|
# plot(c0(1,:),c0(2,:))
|
|
|
|
# * ones((1, len(c0))) # Transform to ellipse for (k,s)
|
|
c1 = np.dot(B, c0) + b[:, None]
|
|
# plot(c1(1,:),c1(2,:)), hold on
|
|
|
|
# Calculate conf.int for lcu
|
|
# Assumtion: lcu is Poisson distributed
|
|
# Poissin distr. approximated by normal when calculating conf. int.
|
|
dXX = 1.64 * sqrt(lcu) # 90 # quantile for lcu
|
|
|
|
lcEstCu = zeros((len(xF), Nc))
|
|
lcEstCl = zeros((len(xF), Nc))
|
|
for i in range(Nc):
|
|
k = c1[0, i]
|
|
s = c1[1, i]
|
|
SF2 = genpareto.sf(xF, k, scale=s)
|
|
lcEstCu[:, i] = (lcu + dXX) * (SF2)
|
|
lcEstCl[:, i] = (lcu - dXX) * (SF2)
|
|
# end
|
|
|
|
lcEst = np.vstack((xF + u, lcu * (SF),
|
|
lcEstCl.min(axis=1), lcEstCu.max(axis=1))).T
|
|
elif dist.startswith('exp'):
|
|
expon = distributions.expon
|
|
phat = expon.fit2(x, floc=0, method=method)
|
|
SF = phat.sf(xF)
|
|
lcEst = np.vstack((xF + u, lcu * (SF))).T
|
|
|
|
elif dist.startswith('ray') or dist.startswith('trun'):
|
|
phat = distributions.truncrayleigh.fit2(x, floc=0, method=method)
|
|
SF = phat.sf(xF)
|
|
# if False:
|
|
# n = len(x)
|
|
# Sx = sum((x + offset) ** 2 - offset ** 2)
|
|
# s = sqrt(Sx / n); # Shape parameter
|
|
# F = -np.expm1(-((xF + offset) ** 2 - offset ** 2) / s ** 2)
|
|
lcEst = np.vstack((xF + u, lcu * (SF))).T
|
|
else:
|
|
raise NotImplementedError('Unknown distribution {}'.format(dist))
|
|
|
|
return lcEst, phat
|
|
# End extrapolate
|
|
|
|
def _make_increasing(self, f, t=None):
|
|
# Makes the signal f strictly increasing.
|
|
|
|
n = len(f)
|
|
if t is None:
|
|
t = np.arange(n)
|
|
ff = [f[0], ]
|
|
tt = [t[0], ]
|
|
|
|
for i in range(1, n):
|
|
if f[i] > ff[-1]:
|
|
ff.append(f[i])
|
|
tt.append(t[i])
|
|
|
|
return np.asarray(ff), np.asarray(tt)
|
|
|
|
def sim(self, ns, alpha):
|
|
"""
|
|
Simulates process with given irregularity factor and crossing spectrum
|
|
|
|
Parameters
|
|
----------
|
|
ns : scalar, integer
|
|
number of sample points.
|
|
alpha : real scalar
|
|
irregularity factor, 0<alpha<1, small alpha gives
|
|
irregular process.
|
|
|
|
Returns
|
|
--------
|
|
ts : timeseries object
|
|
with times and values of the simulated process.
|
|
|
|
Example
|
|
-------
|
|
>>> import wafo.spectrum.models as sm
|
|
>>> from wafo.objects import mat2timeseries
|
|
>>> Sj = sm.Jonswap(Hm0=7)
|
|
>>> S = Sj.tospecdata() #Make spectrum object from numerical values
|
|
>>> alpha = S.characteristic('alpha')[0]
|
|
>>> n = 10000
|
|
>>> xs = S.sim(ns=n)
|
|
>>> ts = mat2timeseries(xs)
|
|
>>> tp = ts.turning_points()
|
|
>>> mm = tp.cycle_pairs()
|
|
>>> lc = mm.level_crossings()
|
|
|
|
>>> xs2 = lc.sim(n,alpha)
|
|
>>> ts2 = mat2timeseries(xs2)
|
|
>>> Se = ts2.tospecdata(L=324)
|
|
|
|
>>> alpha2 = Se.characteristic('alpha')[0]
|
|
>>> np.round(alpha2*10)
|
|
array([ 7.])
|
|
>>> np.abs(alpha-alpha2)<0.03
|
|
array([ True], dtype=bool)
|
|
|
|
>>> lc2 = ts2.turning_points().cycle_pairs().level_crossings()
|
|
|
|
>>> 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()
|
|
|
|
"""
|
|
|
|
# TODO: add a good example
|
|
f = linspace(0, 0.49999, 1000)
|
|
rho_st = 2. * sin(f * pi) ** 2 - 1.
|
|
tmp = alpha * arcsin(sqrt((1. + rho_st) / 2))
|
|
tmp = sin(tmp) ** 2
|
|
a2 = (tmp - rho_st) / (1 - tmp)
|
|
y = vstack((a2 + rho_st, 1 - a2)).min(axis=0)
|
|
maxidx = y.argmax()
|
|
# [maximum,maxidx]=max(y)
|
|
|
|
rho_st = rho_st[maxidx]
|
|
a2 = a2[maxidx]
|
|
a1 = 2. * rho_st + a2 - 1.
|
|
r0 = 1.
|
|
r1 = -a1 / (1. + a2)
|
|
r2 = (a1 ** 2 - a2 - a2 ** 2) / (1 + a2)
|
|
sigma2 = r0 + a1 * r1 + a2 * r2
|
|
# randn = np.random.randn
|
|
e = randn(ns) * sqrt(sigma2)
|
|
e[:2] = 0.0
|
|
L0 = randn(1)
|
|
L0 = hstack((L0, r1 * L0 + sqrt(1 - r2 ** 2) * randn(1)))
|
|
# Simulate the process, starting in L0
|
|
z0 = lfilter([1, a1, a2], ones(1), L0)
|
|
L, unused_zf = lfilter(ones(1), [1, a1, a2], e, axis=0, zi=z0)
|
|
|
|
epsilon = 1.01
|
|
min_L = min(L)
|
|
max_L = max(L)
|
|
maxi = max(abs(r_[min_L, max_L])) * epsilon
|
|
mini = -maxi
|
|
nu = 101
|
|
u = linspace(mini, maxi, nu)
|
|
G = cdfnorm(u) # (1 + erf(u / sqrt(2))) / 2
|
|
G = G * (1 - G)
|
|
|
|
x = linspace(0, r1, 100)
|
|
factor1 = 1. / sqrt(1 - x ** 2)
|
|
factor2 = 1. / (1 + x)
|
|
integral = zeros(u.shape, dtype=float)
|
|
for i in range(nu):
|
|
y = factor1 * exp(-u[i] * u[i] * factor2)
|
|
integral[i] = trapz(y, x)
|
|
# end
|
|
G = G - integral / (2 * pi)
|
|
G = G / max(G)
|
|
|
|
Z = ((u >= 0) * 2 - 1) * sqrt(-2 * log(G))
|
|
|
|
sumcr = trapz(self.data, self.args)
|
|
lc = self.data / sumcr
|
|
lc1 = self.args
|
|
mcr = trapz(lc1 * lc, lc1) if self.mean is None else self.mean
|
|
if self.sigma is None:
|
|
scr = sqrt(trapz(lc1 ** 2 * lc, lc1) - mcr ** 2)
|
|
else:
|
|
scr = self.sigma
|
|
lc2 = LevelCrossings(lc, lc1, mean=mcr, sigma=scr, intensity=True)
|
|
|
|
g = lc2.trdata()[0]
|
|
|
|
f = g.gauss2dat(Z)
|
|
G = TrData(f, u)
|
|
|
|
process = G.dat2gauss(L)
|
|
return np.vstack((arange(len(process)), process)).T
|
|
|
|
#
|
|
#
|
|
# Check the result without reference to getrfc:
|
|
# LCe = dat2lc(process)
|
|
# max(lc(:,2))
|
|
# max(LCe(:,2))
|
|
#
|
|
# clf
|
|
# plot(lc(:,1),lc(:,2)/max(lc(:,2)))
|
|
# hold on
|
|
# plot(LCe(:,1),LCe(:,2)/max(LCe(:,2)),'-.')
|
|
# title('Relative crossing intensity')
|
|
#
|
|
# %% Plot made by the function funplot_4, JE 970707
|
|
# %param = [min(process(:,2)) max(process(:,2)) 100]
|
|
# %plot(lc(:,1),lc(:,2)/max(lc(:,2)))
|
|
# %hold on
|
|
# %plot(levels(param),mu/max(mu),'--')
|
|
# %hold off
|
|
# %title('Crossing intensity')
|
|
# %watstamp
|
|
#
|
|
# % Temporarily
|
|
# %funplot_4(lc,param,mu)
|
|
def trdata(self, mean=None, sigma=None, **options):
|
|
'''
|
|
Estimate transformation, g, from observed crossing intensity, version2.
|
|
|
|
Assumption: a Gaussian process, Y, is related to the
|
|
non-Gaussian process, X, by Y = g(X).
|
|
|
|
Parameters
|
|
----------
|
|
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.
|
|
(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 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 (eg. 1000) will
|
|
significantly speed up the estimation for long time series without
|
|
loosing any accuracy. NTR should be chosen greater than PARAM(3).
|
|
(default inf)
|
|
|
|
Returns
|
|
-------
|
|
gs, ge : TrData objects
|
|
smoothed and empirical estimate of the transformation g.
|
|
|
|
|
|
Notes
|
|
-----
|
|
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 GSM or set a larger variance for GVAR.
|
|
If X(t) is likely to cross levels higher than 5 standard deviations
|
|
then the vector param has to be modified. For example if X(t) 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, iseed=10)
|
|
>>> ts = mat2timeseries(xs)
|
|
>>> tp = ts.turning_points()
|
|
>>> mm = tp.cycle_pairs()
|
|
>>> lc = mm.level_crossings()
|
|
>>> g0, g0emp = lc.trdata(plotflag=0)
|
|
>>> g1, g1emp = lc.trdata(gvar=0.5 ) # Equal weight on all points
|
|
>>> g2, g2emp = lc.trdata(gvar=[3.5, 0.5, 3.5]) # Less weight on ends
|
|
>>> int(S.tr.dist2gauss()*100)
|
|
141
|
|
>>> int(g0emp.dist2gauss()*100)
|
|
380995
|
|
>>> int(g0.dist2gauss()*100)
|
|
143
|
|
>>> int(g1.dist2gauss()*100)
|
|
162
|
|
>>> int(g2.dist2gauss()*100)
|
|
120
|
|
|
|
g0.plot() # Check the fit.
|
|
|
|
See also
|
|
troptset, dat2tr, trplot, findcross, smooth
|
|
|
|
NB! the transformated data will be N(0,1)
|
|
|
|
Reference
|
|
---------
|
|
Rychlik , I., Johannesson, P., and Leadbetter, M.R. (1997)
|
|
"Modelling and statistical analysis of ocean wavedata
|
|
using a transformed Gaussian process",
|
|
Marine structures, Design, Construction and Safety,
|
|
Vol 10, pp 13--47
|
|
'''
|
|
estimate = TransformEstimator(**options)
|
|
return estimate._trdata_lc(self, mean, sigma)
|
|
|
|
|
|
class CyclePairs(PlotData):
|
|
'''
|
|
Container class for Cycle Pairs data objects in WAFO
|
|
|
|
Member variables
|
|
----------------
|
|
data : array_like
|
|
args : vector for 1D
|
|
|
|
Examples
|
|
--------
|
|
>>> import wafo.data
|
|
>>> import wafo.objects as wo
|
|
>>> x = wafo.data.sea()
|
|
>>> ts = wo.mat2timeseries(x)
|
|
|
|
>>> tp = ts.turning_points()
|
|
>>> 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')
|
|
'''
|
|
|
|
def __init__(self, *args, **kwds):
|
|
self.kind = kwds.pop('kind', 'min2max')
|
|
self.sigma = kwds.pop('sigma', None)
|
|
self.mean = kwds.pop('mean', None)
|
|
self.time = kwds.pop('time', 1)
|
|
|
|
options = dict(title=self.kind + ' cycle pairs',
|
|
xlab='min', ylab='max',
|
|
plot_args=['b.'])
|
|
options.update(**kwds)
|
|
super(CyclePairs, self).__init__(*args, **options)
|
|
|
|
def amplitudes(self):
|
|
return (self.data - self.args) / 2.
|
|
|
|
def damage(self, beta, K=1):
|
|
"""
|
|
Calculates the total Palmgren-Miner damage of cycle pairs.
|
|
|
|
Parameters
|
|
----------
|
|
beta : array-like, size m
|
|
Beta-values, material parameter.
|
|
K : scalar, optional
|
|
K-value, material parameter.
|
|
|
|
Returns
|
|
-------
|
|
D : ndarray, size m
|
|
Damage.
|
|
|
|
Notes
|
|
-----
|
|
The damage is calculated according to
|
|
D[i] = sum ( K * a**beta[i] ), with a = (max-min)/2
|
|
|
|
Examples
|
|
--------
|
|
>>> import wafo
|
|
>>> from matplotlib import pyplot as plt
|
|
>>> ts = wafo.objects.mat2timeseries(wafo.data.sea())
|
|
>>> tp = ts.turning_points()
|
|
>>> mm = tp.cycle_pairs()
|
|
|
|
>>> bv = range(3,9)
|
|
>>> D = mm.damage(beta=bv)
|
|
>>> 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
|
|
--------
|
|
SurvivalCycleCount
|
|
"""
|
|
amp = abs(self.amplitudes())
|
|
return atleast_1d([K * np.sum(amp ** betai) for betai in beta])
|
|
|
|
def get_minima_and_maxima(self):
|
|
index, = nonzero(self.args <= self.data)
|
|
if index.size == 0:
|
|
index, = nonzero(self.args >= self.data)
|
|
M = self.args[index]
|
|
m = self.data[index]
|
|
else:
|
|
m = self.args[index]
|
|
M = self.data[index]
|
|
return m, M
|
|
|
|
def level_crossings(self, kind='uM', intensity=False):
|
|
""" Return level crossing spectrum from a cycle count.
|
|
|
|
Parameters
|
|
----------
|
|
kind : int or string
|
|
defining crossing type, options are
|
|
0,'u' : only upcrossings.
|
|
1,'uM' : upcrossings and maxima (default).
|
|
2,'umM': upcrossings, minima, and maxima.
|
|
3,'um' : upcrossings and minima.
|
|
intensity : bool
|
|
True if level crossing intensity spectrum
|
|
False if level crossing count spectrum
|
|
Return
|
|
------
|
|
lc : level crossing object
|
|
with levels and number of upcrossings.
|
|
|
|
Calculates the number of upcrossings from a cycle pairs, e.g.
|
|
min2Max cycles or rainflow cycles.
|
|
|
|
Example:
|
|
--------
|
|
>>> import wafo
|
|
>>> ts = wafo.objects.mat2timeseries(wafo.data.sea())
|
|
>>> tp = ts.turning_points()
|
|
>>> mm = tp.cycle_pairs()
|
|
>>> lc = mm.level_crossings()
|
|
|
|
h = mm.plot(marker='.')
|
|
h2 = lc.plot()
|
|
|
|
See also
|
|
--------
|
|
TurningPoints
|
|
LevelCrossings
|
|
"""
|
|
|
|
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, "u", "uM",'
|
|
' "umM", "um"). Got kind = {}'.format(kind))
|
|
|
|
m, M = self.get_minima_and_maxima()
|
|
|
|
ncc = len(m)
|
|
minima = vstack((m, ones(ncc), zeros(ncc), ones(ncc)))
|
|
maxima = vstack((M, -ones(ncc), ones(ncc), zeros(ncc)))
|
|
|
|
extremes = hstack((maxima, minima))
|
|
index = extremes[0].argsort()
|
|
extremes = extremes[:, index]
|
|
|
|
ii = 0
|
|
n = extremes.shape[1]
|
|
extr = zeros((4, n))
|
|
extr[:, 0] = extremes[:, 0]
|
|
for i in range(1, n):
|
|
if extremes[0, i] == extr[0, ii]:
|
|
extr[1:4, ii] = extr[1:4, ii] + extremes[1:4, i]
|
|
else:
|
|
ii += 1
|
|
extr[:, ii] = extremes[:, i]
|
|
|
|
nx = extr[0].argmax() + 1
|
|
levels = extr[0, 0:nx]
|
|
if defnr == 2: # This are upcrossings + maxima
|
|
dcount = cumsum(extr[1, 0:nx]) + extr[2, 0:nx] - extr[3, 0:nx]
|
|
elif defnr == 4: # This are upcrossings + minima
|
|
dcount = cumsum(extr[1, 0:nx])
|
|
dcount[nx - 1] = dcount[nx - 2]
|
|
elif defnr == 1: # This are only upcrossings
|
|
dcount = cumsum(extr[1, 0:nx]) - extr[3, 0:nx]
|
|
elif defnr == 3: # This are upcrossings + minima + maxima
|
|
dcount = cumsum(extr[1, 0:nx]) + extr[2, 0:nx]
|
|
ylab = 'Count'
|
|
if intensity:
|
|
dcount = dcount / self.time
|
|
ylab = 'Intensity [count/sec]'
|
|
return LevelCrossings(dcount, levels, mean=self.mean, sigma=self.sigma,
|
|
ylab=ylab, intensity=intensity)
|
|
|
|
|
|
class TurningPoints(PlotData):
|
|
|
|
'''
|
|
Container class for Turning Points data objects in WAFO
|
|
|
|
Member variables
|
|
----------------
|
|
data : array_like
|
|
args : vector for 1D
|
|
|
|
Examples
|
|
--------
|
|
>>> import wafo.data
|
|
>>> import wafo.objects as wo
|
|
>>> x = wafo.data.sea()
|
|
>>> ts = wo.mat2timeseries(x)
|
|
|
|
>>> tp = ts.turning_points()
|
|
>>> np.allclose(tp.data[:5],
|
|
... [-1.2004945 , 0.83950546, -0.09049454, -0.02049454, -0.09049454])
|
|
True
|
|
|
|
h1 = tp.plot(marker='x')
|
|
'''
|
|
|
|
def __init__(self, *args, **kwds):
|
|
self.name_ = kwds.pop('name', 'WAFO TurningPoints Object')
|
|
self.sigma = kwds.pop('sigma', None)
|
|
self.mean = kwds.pop('mean', None)
|
|
|
|
options = dict(title='Turning points')
|
|
options.update(**kwds)
|
|
super(TurningPoints, self).__init__(*args, **options)
|
|
|
|
if not any(self.args):
|
|
n = len(self.data)
|
|
self.args = range(0, n)
|
|
else:
|
|
self.args = ravel(self.args)
|
|
self.data = ravel(self.data)
|
|
|
|
def rainflow_filter(self, h=0.0, method='clib'):
|
|
'''
|
|
Return rainflow filtered turning points (tp).
|
|
|
|
Parameters
|
|
----------
|
|
h : scalar
|
|
a threshold
|
|
if h<=0, then tp is a sequence of turning points (default)
|
|
if h>0, then all rainflow cycles with height smaller than
|
|
h are removed.
|
|
|
|
Returns
|
|
-------
|
|
tp : TurningPoints object
|
|
with times and turning points.
|
|
|
|
Example:
|
|
>>> import wafo.data
|
|
>>> x = wafo.data.sea()
|
|
>>> x1 = x[:200,:]
|
|
>>> ts1 = mat2timeseries(x1)
|
|
>>> tp = ts1.turning_points(wavetype='Mw')
|
|
>>> tph = tp.rainflow_filter(h=0.3)
|
|
>>> 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.')
|
|
|
|
See also
|
|
---------
|
|
findcross,
|
|
findrfc
|
|
findtp
|
|
'''
|
|
ind = findrfc(self.data, max(h, 0.0), method)
|
|
try:
|
|
t = self.args[ind]
|
|
except:
|
|
t = ind
|
|
mean = self.mean
|
|
sigma = self.sigma
|
|
return TurningPoints(self.data[ind], t, mean=mean, sigma=sigma)
|
|
|
|
def cycle_pairs(self, h=0, kind='min2max', method='clib'):
|
|
""" Return min2Max or Max2min cycle pairs from turning points
|
|
|
|
Parameters
|
|
----------
|
|
kind : string
|
|
type of cycles to return options are 'min2max' or 'max2min'
|
|
method : string
|
|
specify which library to use
|
|
'clib' for wafo's c_library
|
|
'None' for wafo's Python functions
|
|
|
|
Return
|
|
------
|
|
mm : cycles object
|
|
with min2Max or Max2min cycle pairs.
|
|
|
|
Example
|
|
-------
|
|
>>> import wafo
|
|
>>> x = wafo.data.sea()
|
|
>>> ts = wafo.objects.mat2timeseries(x)
|
|
>>> tp = ts.turning_points()
|
|
>>> mM = tp.cycle_pairs()
|
|
>>> np.allclose(mM.data[:5], [ 0.83950546, -0.02049454, -0.04049454,
|
|
... 0.25950546, -0.08049454])
|
|
True
|
|
|
|
>>> h = mM.plot(marker='x')
|
|
|
|
|
|
See also
|
|
--------
|
|
TurningPoints
|
|
SurvivalCycleCount
|
|
"""
|
|
|
|
if h > 0:
|
|
ind = findrfc(self.data, h, method=method)
|
|
data = self.data[ind]
|
|
else:
|
|
data = self.data
|
|
if data[0] > data[1]:
|
|
im = 1
|
|
iM = 0
|
|
else:
|
|
im = 0
|
|
iM = 1
|
|
|
|
# Extract min-max and max-min cycle pairs
|
|
if kind.lower().startswith('min2max'):
|
|
m = data[im:-1:2]
|
|
M = data[im + 1::2]
|
|
else:
|
|
kind = 'max2min'
|
|
M = data[iM:-1:2]
|
|
m = data[iM + 1::2]
|
|
|
|
time = self.args[-1] - self.args[0]
|
|
|
|
return CyclePairs(M, m, kind=kind, mean=self.mean, sigma=self.sigma,
|
|
time=time)
|
|
|
|
def cycle_astm(self):
|
|
"""
|
|
Rainflow counted cycles according to Nieslony's ASTM implementation
|
|
|
|
Parameters
|
|
----------
|
|
|
|
Returns
|
|
-------
|
|
sig_rfc : array-like
|
|
array of shape (n,3) with:
|
|
sig_rfc[:,0] Cycles amplitude
|
|
sig_rfc[:,1] Cycles mean value
|
|
sig_rfc[:,2] Cycle type, half (=0.5) or full (=1.0)
|
|
|
|
References
|
|
----------
|
|
Adam Nieslony, "Determination of fragments of multiaxial service
|
|
loading strongly influencing the fatigue of machine components",
|
|
Mechanical Systems and Signal Processing 23, no. 8 (2009): 2712-2721.
|
|
|
|
and is based on the following standard:
|
|
ASTM E 1049-85 (Reapproved 1997), Standard practices for cycle counting
|
|
in fatigue analysis, in: Annual Book of ASTM Standards,
|
|
vol. 03.01, ASTM, Philadelphia, 1999, pp. 710-718.
|
|
|
|
Copyright (c) 1999-2002 by Adam Nieslony
|
|
Ported to Python by David Verelst
|
|
|
|
Example
|
|
-------
|
|
>>> import wafo
|
|
>>> x = wafo.data.sea()
|
|
>>> sig_ts = wafo.objects.mat2timeseries(x)
|
|
>>> sig_tp = sig_ts.turning_points(h=0, wavetype='astm')
|
|
>>> sig_cp = sig_tp.cycle_astm()
|
|
"""
|
|
|
|
# output of Nieslony's algorithm is organised differently with
|
|
# respect to wafo's approach
|
|
# TODO: integrate ASTM method into the CyclyPairs class?
|
|
return findrfc_astm(self.data)
|
|
|
|
|
|
def mat2timeseries(x):
|
|
"""
|
|
Convert 2D arrays to TimeSeries object
|
|
assuming 1st column is time and the remaining columns contain data.
|
|
"""
|
|
return TimeSeries(x[:, 1::], x[:, 0].ravel())
|
|
|
|
|
|
class TimeSeries(PlotData):
|
|
'''
|
|
Container class for 1D TimeSeries data objects in WAFO
|
|
|
|
Member variables
|
|
----------------
|
|
data : array_like
|
|
args : vector for 1D, list of vectors for 2D, 3D, ...
|
|
|
|
sensortypes : list of integers or strings
|
|
sensor type for time series (default ['n'] : Surface elevation)
|
|
see sensortype for more options
|
|
position : vector of size 3
|
|
instrument position relative to the coordinate system
|
|
|
|
Examples
|
|
--------
|
|
>>> import wafo.data
|
|
>>> import wafo.objects as wo
|
|
>>> x = wafo.data.sea()
|
|
>>> ts = wo.mat2timeseries(x)
|
|
>>> rf = ts.tocovdata(lag=150)
|
|
|
|
>>> S = ts.tospecdata()
|
|
>>> tp = ts.turning_points()
|
|
>>> mm = tp.cycle_pairs()
|
|
>>> lc = mm.level_crossings()
|
|
|
|
h = rf.plot()
|
|
h1 = mm.plot(marker='x')
|
|
h2 = lc.plot()
|
|
'''
|
|
|
|
def __init__(self, *args, **kwds):
|
|
self.name_ = kwds.pop('name', 'WAFO TimeSeries Object')
|
|
self.sensortypes = kwds.pop('sensortypes', ['n', ])
|
|
self.position = kwds.pop('position', [zeros(3), ])
|
|
|
|
super(TimeSeries, self).__init__(*args, **kwds)
|
|
|
|
if not any(self.args):
|
|
n = len(self.data)
|
|
self.args = range(0, n)
|
|
|
|
def sampling_period(self):
|
|
'''
|
|
Returns sampling interval
|
|
|
|
Returns
|
|
-------
|
|
dt : scalar
|
|
sampling interval, unit:
|
|
[s] if lagtype=='t'
|
|
[m] otherwise
|
|
|
|
See also
|
|
'''
|
|
t_vec = self.args
|
|
dt1 = t_vec[1] - t_vec[0]
|
|
n = len(t_vec) - 1
|
|
t = t_vec[-1] - t_vec[0]
|
|
dt = t / n
|
|
if abs(dt - dt1) > 1e-10:
|
|
warnings.warn('Data is not uniformly sampled!')
|
|
return dt
|
|
|
|
def tocovdata(self, lag=None, tr=None, detrend=detrend_mean,
|
|
window='boxcar', flag='biased', norm=False, dt=None):
|
|
'''
|
|
Return auto covariance function from data.
|
|
|
|
Parameters
|
|
----------
|
|
lag : scalar, int
|
|
maximum time-lag for which the ACF is estimated. (Default lag=n-1)
|
|
flag : string, 'biased' or 'unbiased'
|
|
If 'unbiased' scales the raw correlation by 1/(n-abs(k)),
|
|
where k is the index into the result, otherwise scales the raw
|
|
cross-correlation by 1/n. (default)
|
|
norm : bool
|
|
True if normalize output to one
|
|
dt : scalar
|
|
time-step between data points (default see sampling_period).
|
|
|
|
Return
|
|
-------
|
|
R : CovData1D object
|
|
with attributes:
|
|
data : ACF vector length L+1
|
|
args : time lags length L+1
|
|
sigma : estimated large lag standard deviation of the estimate
|
|
assuming x is a Gaussian process:
|
|
if R(k)=0 for all lags k>q then an approximation
|
|
of the variance for large samples due to Bartlett
|
|
var(R(k))=1/N*(R(0)^2+2*R(1)^2+2*R(2)^2+ ..+2*R(q)^2)
|
|
for k>q and where N=length(x). Special case is
|
|
white noise where it equals R(0)^2/N for k>0
|
|
norm : bool
|
|
If false indicating that R is not normalized
|
|
|
|
Example:
|
|
--------
|
|
>>> import wafo.data
|
|
>>> import wafo.objects as wo
|
|
>>> x = wafo.data.sea()
|
|
>>> ts = wo.mat2timeseries(x)
|
|
>>> acf = ts.tocovdata(150)
|
|
>>> np.allclose(acf.data[:3], [ 0.22368637, 0.20838473, 0.17110733])
|
|
True
|
|
|
|
h = acf.plot()
|
|
'''
|
|
estimate_cov = _wafocov_estimation.CovarianceEstimator(
|
|
lag=lag, tr=tr, detrend=detrend, window=window, flag=flag,
|
|
norm=norm, dt=dt)
|
|
return estimate_cov(self)
|
|
|
|
def _get_bandwidth_and_dof(self, wname, n, L, dt, ftype='w'):
|
|
'''Returns bandwidth (rad/sec) and degrees of freedom
|
|
used in chi^2 distribution
|
|
'''
|
|
if isinstance(wname, tuple):
|
|
wname = wname[0]
|
|
dof = int(dict(parzen=3.71, hanning=2.67,
|
|
bartlett=3).get(wname, np.nan) * n/L)
|
|
Be = dict(parzen=1.33, hanning=1,
|
|
bartlett=1.33).get(wname, np.nan) * 2 * pi / (L*dt)
|
|
if ftype == 'f':
|
|
Be = Be / (2 * pi) # bandwidth in Hz
|
|
return Be, dof
|
|
|
|
def tospecdata(self, L=None, tr=None, method='cov', detrend=detrend_mean,
|
|
window='parzen', noverlap=0, ftype='w', alpha=None):
|
|
'''
|
|
Estimate one-sided spectral density from data.
|
|
|
|
Parameters
|
|
----------
|
|
L : scalar integer
|
|
maximum lag size of the window function. As L decreases the
|
|
estimate becomes smoother and Bw increases. If we want to resolve
|
|
peaks in S which is Bf (Hz or rad/sec) apart then Bw < Bf. If no
|
|
value is given the lag size is set to be the lag where the auto
|
|
correlation is less than 2 standard deviations. (maximum 300)
|
|
tr : transformation object
|
|
the transformation assuming that x is a sample of a transformed
|
|
Gaussian process. If g is None then x is a sample of a Gaussian
|
|
process (Default)
|
|
method : string
|
|
defining estimation method. Options are
|
|
'cov' : Frequency smoothing using the window function
|
|
on the estimated autocovariance function. (default)
|
|
'psd' : Welch's averaged periodogram method with no overlapping
|
|
batches
|
|
detrend : function
|
|
defining detrending performed on the signal before estimation.
|
|
(default detrend_mean)
|
|
window : vector of length NFFT or function
|
|
To create window vectors see numpy.blackman, numpy.hamming,
|
|
numpy.bartlett, scipy.signal, scipy.signal.get_window etc.
|
|
noverlap : scalar int
|
|
gives the length of the overlap between segments.
|
|
ftype : character
|
|
defining frequency type: 'w' or 'f' (default 'w')
|
|
|
|
Returns
|
|
---------
|
|
spec : SpecData1D object
|
|
|
|
Example
|
|
-------
|
|
>>> import wafo.data as wd
|
|
>>> import wafo.objects as wo
|
|
>>> x = wd.sea()
|
|
>>> ts = wo.mat2timeseries(x)
|
|
>>> S0 = ts.tospecdata(method='psd')
|
|
>>> np.allclose(S0.data[21:25],
|
|
... (0.2543896 , 0.26366755, 0.23372824, 0.19459349), rtol=1e-2)
|
|
True
|
|
>>> S = ts.tospecdata()
|
|
|
|
>>> np.allclose(S.data[21:25],
|
|
... (0.00207876, 0.0025113 , 0.00300008, 0.00351852), rtol=1e-2)
|
|
True
|
|
>>> h = S.plot()
|
|
|
|
See also
|
|
--------
|
|
dat2tr, dat2cov
|
|
|
|
References:
|
|
-----------
|
|
Georg Lindgren and Holger Rootzen (1986)
|
|
"Stationara stokastiska processer", pp 173--176.
|
|
|
|
Gareth Janacek and Louise Swift (1993)
|
|
"TIME SERIES forecasting, simulation, applications",
|
|
pp 75--76 and 261--268
|
|
|
|
Emanuel Parzen (1962),
|
|
"Stochastic Processes", HOLDEN-DAY,
|
|
pp 66--103
|
|
'''
|
|
|
|
nugget = 1e-12
|
|
rate = 2 # interpolationrate for frequency
|
|
dt = self.sampling_period()
|
|
|
|
yy = self.data.ravel()
|
|
if tr is not None:
|
|
yy = tr.dat2gauss(yy)
|
|
yy = detrend(yy) if hasattr(detrend, '__call__') else yy
|
|
n = len(yy)
|
|
|
|
estimate_L = L is None
|
|
if method == 'cov' or estimate_L:
|
|
tsy = TimeSeries(yy, self.args)
|
|
R = tsy.tocovdata(lag=L, window=window)
|
|
L = len(R.data) - 1
|
|
if method == 'cov':
|
|
# add a nugget effect to ensure that round off errors
|
|
# do not result in negative spectral estimates
|
|
spec = R.tospecdata(rate=rate, nugget=nugget)
|
|
L = min(L, n - 1)
|
|
if method == 'psd':
|
|
nfft = 2 ** nextpow2(L)
|
|
pad_to = rate * nfft # Interpolate the spectrum with rate
|
|
f, S = welch(yy, fs=1.0 / dt, window=window, nperseg=nfft,
|
|
noverlap=noverlap, nfft=pad_to, detrend=detrend,
|
|
return_onesided=True, scaling='density', axis=-1)
|
|
# S, f = psd(yy, Fs=1. / dt, NFFT=nfft, detrend=detrend,
|
|
# window=win, noverlap=noverlap, pad_to=pad_to,
|
|
# scale_by_freq=True)
|
|
fact = 2.0 * pi
|
|
w = fact * f
|
|
spec = _wafospec.SpecData1D(S / fact, w)
|
|
elif method == 'cov':
|
|
pass
|
|
else:
|
|
raise ValueError('Unknown method (%s)' % method)
|
|
|
|
Be, dof = self._get_bandwidth_and_dof(window, n, L, dt, ftype)
|
|
spec.Bw = Be
|
|
|
|
if alpha is not None:
|
|
# Confidence interval constants
|
|
spec.CI = [dof / _invchi2(1 - alpha / 2, dof),
|
|
dof / _invchi2(alpha / 2, dof)]
|
|
|
|
spec.tr = tr
|
|
spec.L = L
|
|
spec.norm = False
|
|
spec.note = 'method=%s' % method
|
|
return spec
|
|
|
|
def trdata(self, method='nonlinear', **options):
|
|
'''
|
|
Estimate transformation, g, from data.
|
|
|
|
Parameters
|
|
----------
|
|
method : string defining transform based on:
|
|
'nonlinear' : smoothed crossing intensity (default)
|
|
'mnonlinear': smoothed marginal distribution
|
|
'hermite' : cubic Hermite polynomial
|
|
'ochi' : exponential function
|
|
'linear' : identity.
|
|
|
|
options : keyword with the following fields:
|
|
csm, gsm : real scalars
|
|
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 (default see lc2tr)
|
|
which defines the region of variation of the data x.
|
|
plotflag : int
|
|
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: int
|
|
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: real scalar
|
|
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)
|
|
|
|
Returns
|
|
-------
|
|
tr, tr_emp : TrData objects
|
|
with the smoothed and empirical transformation, respectively.
|
|
|
|
|
|
TRDATA 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, iseed=10)
|
|
>>> ts = mat2timeseries(xs)
|
|
>>> g0, g0emp = ts.trdata(plotflag=0)
|
|
>>> g1, g1emp = ts.trdata(method='mnonlinear', gvar=0.5 )
|
|
>>> g2, g2emp = ts.trdata(method='nonlinear', gvar=[3.5, 0.5, 3.5])
|
|
>>> 100 < S.tr.dist2gauss()*100 < 200
|
|
True
|
|
>>> 2000 < g0emp.dist2gauss() < 4000
|
|
True
|
|
>>> 80 < g0.dist2gauss()*100 < 150
|
|
True
|
|
>>> 50 < g1.dist2gauss()*100 < 100
|
|
True
|
|
>>> 70 < g2.dist2gauss()*100 < 140
|
|
True
|
|
|
|
See also
|
|
--------
|
|
LevelCrossings.trdata
|
|
wafo.transform.models
|
|
|
|
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
|
|
'''
|
|
estimate = TransformEstimator(method=method, **options)
|
|
return estimate.trdata(self)
|
|
|
|
def turning_points(self, h=0.0, wavetype=None):
|
|
'''
|
|
Return turning points (tp) from data, optionally rainflowfiltered.
|
|
|
|
Parameters
|
|
----------
|
|
h : scalar
|
|
a threshold
|
|
if h<=0, then tp is a sequence of turning points (default)
|
|
if h>0, then all rainflow cycles with height smaller than
|
|
h are removed.
|
|
wavetype : string
|
|
defines the type of wave. Possible options are
|
|
'astm' 'mw' 'Mw' or 'none'.
|
|
If None all rainflow filtered min and max
|
|
will be returned, otherwise only the rainflow filtered
|
|
min and max, which define a wave according to the
|
|
wave definition, will be returned.
|
|
'astm' forces to have the first data point of the load history as
|
|
the first turning point. To be used in combination with
|
|
TurningPoints.cycle_astm()
|
|
|
|
Returns
|
|
-------
|
|
tp : TurningPoints object
|
|
with times and turning points.
|
|
|
|
Example:
|
|
>>> import wafo.data
|
|
>>> x = wafo.data.sea()
|
|
>>> x1 = x[:200,:]
|
|
>>> ts1 = mat2timeseries(x1)
|
|
>>> tp = ts1.turning_points(wavetype='Mw')
|
|
>>> tph = ts1.turning_points(h=0.3,wavetype='Mw')
|
|
>>> np.allclose(tph.data[:3], [ 0.83950546, -0.16049454, 0.25950546])
|
|
True
|
|
|
|
hs = ts1.plot()
|
|
hp = tp.plot('ro')
|
|
hph = tph.plot('k.')
|
|
|
|
See also
|
|
---------
|
|
findcross,
|
|
findrfc
|
|
findtp
|
|
'''
|
|
ind = findtp(self.data, max(h, 0.0), wavetype)
|
|
try:
|
|
t = self.args[ind]
|
|
except:
|
|
t = ind
|
|
mean = self.data.mean()
|
|
sigma = self.data.std()
|
|
return TurningPoints(self.data[ind], t, mean=mean, sigma=sigma)
|
|
|
|
def trough_crest(self, v=None, wavetype=None):
|
|
"""
|
|
Return trough and crest turning points
|
|
|
|
Parameters
|
|
-----------
|
|
v : scalar
|
|
reference level (default v = mean of x).
|
|
|
|
wavetype : string
|
|
defines the type of wave. Possible options are
|
|
'dw', 'uw', 'tw', 'cw' or None.
|
|
If None indices to all troughs and crests will be returned,
|
|
otherwise only the paired ones will be returned
|
|
according to the wavedefinition.
|
|
|
|
Returns
|
|
--------
|
|
tc : TurningPoints object
|
|
with trough and crest turningpoints
|
|
"""
|
|
ind = findtc(self.data, v, wavetype)[0]
|
|
try:
|
|
t = self.args[ind]
|
|
except:
|
|
t = ind
|
|
mean = self.data.mean()
|
|
sigma = self.data.std()
|
|
return TurningPoints(self.data[ind], t, mean=mean, sigma=sigma)
|
|
|
|
def wave_parameters(self, rate=1):
|
|
'''
|
|
Returns several wave parameters from data.
|
|
|
|
Parameters
|
|
----------
|
|
rate : scalar integer
|
|
interpolation rate. Interpolates with spline if greater than one.
|
|
|
|
Returns
|
|
-------
|
|
parameters : dict
|
|
wave parameters such as
|
|
Ac, At : Crest and trough amplitude, respectively
|
|
Tcf, Tcb : Crest front and crest (rear) back period, respectively
|
|
Hu, Hd : zero-up- and down-crossing wave height, respectively.
|
|
Tu, Td : zero-up- and down-crossing wave period, respectively.
|
|
|
|
The definition of g, Ac,At, Tcf, etc. are given in gravity and
|
|
wafo.definitions.
|
|
|
|
Example
|
|
-------
|
|
>>> import wafo.data as wd
|
|
>>> import wafo.objects as wo
|
|
>>> x = wd.sea()
|
|
>>> ts = wo.mat2timeseries(x)
|
|
>>> wp = ts.wave_parameters()
|
|
>>> true_wp = {'Ac':[ 0.25950546, 0.34950546],
|
|
... 'At': [ 0.16049454, 0.43049454],
|
|
... 'Hu': [ 0.69, 0.86],
|
|
... 'Hd': [ 0.42, 0.78],
|
|
... 'Tu': [ 6.10295202, 3.36978685],
|
|
... 'Td': [ 3.84377468, 6.35707656],
|
|
... 'Tcf': [ 0.42656819, 0.57361617],
|
|
... 'Tcb': [ 0.93355982, 1.04063638]}
|
|
>>> for name in ['Ac', 'At', 'Hu', 'Hd', 'Tu', 'Td', 'Tcf', 'Tcb']:
|
|
... np.allclose(wp[name][:2], true_wp[name])
|
|
True
|
|
True
|
|
True
|
|
True
|
|
True
|
|
True
|
|
True
|
|
True
|
|
|
|
import pylab as plt
|
|
h = plt.plot(wp['Td'],wp['Hd'],'.')
|
|
h = plt.xlabel('Td [s]')
|
|
h = plt.ylabel('Hd [m]')
|
|
|
|
|
|
See also
|
|
--------
|
|
wafo.definitions
|
|
'''
|
|
dT = self.sampling_period()/np.maximum(rate, 1)
|
|
xi, ti = self._interpolate(rate)
|
|
|
|
tc_ind, z_ind = findtc(xi, v=0, kind='tw')
|
|
tc_a = xi[tc_ind]
|
|
tc_t = ti[tc_ind]
|
|
Ac = tc_a[1::2] # crest amplitude
|
|
At = -tc_a[0::2] # trough amplitude
|
|
Hu = Ac + At[1:]
|
|
Hd = Ac + At[:-1]
|
|
tu = ecross(ti, xi, z_ind[1::2], v=0)
|
|
Tu = diff(tu) # Period zero-upcrossing waves
|
|
td = ecross(ti, xi, z_ind[::2], v=0)
|
|
Td = diff(td) # Period zero-downcrossing waves
|
|
Tcf = tc_t[1::2] - tu[:-1]
|
|
Tcf[(Tcf == 0)] = dT # avoiding division by zero
|
|
Tcb = td[1:] - tc_t[1::2]
|
|
Tcb[(Tcb == 0)] = dT # avoiding division by zero
|
|
return dict(Ac=Ac, At=At, Hu=Hu, Hd=Hd, Tu=Tu, Td=Td, Tcf=Tcf, Tcb=Tcb)
|
|
|
|
def wave_height_steepness(self, kind='Vcf', rate=1, g=None):
|
|
'''
|
|
Returns waveheights and steepnesses from data.
|
|
|
|
Parameters
|
|
----------
|
|
rate : scalar integer
|
|
interpolation rate. Interpolates with spline if greater than one.
|
|
|
|
kind : scalar integer (default 1)
|
|
0 max(Vcf, Vcb) and corresponding wave height Hd or Hu in H
|
|
1 crest front (rise) speed (Vcf) in S and wave height Hd in H.
|
|
-1 crest back (fall) speed (Vcb) in S and waveheight Hu in H.
|
|
2 crest front steepness in S and the wave height Hd in H.
|
|
-2 crest back steepness in S and the wave height Hu in H.
|
|
3 total wave steepness in S and the wave height Hd in H
|
|
for zero-downcrossing waves.
|
|
-3 total wave steepness in S and the wave height Hu in H.
|
|
for zero-upcrossing waves.
|
|
Returns
|
|
-------
|
|
S, H = Steepness and the corresponding wave height according to kind
|
|
|
|
|
|
The parameters are calculated as follows:
|
|
Crest front speed (velocity) = Vcf = Ac/Tcf
|
|
Crest back speed (velocity) = Vcb = Ac/Tcb
|
|
Crest front steepness = 2*pi*Ac./Td/Tcf/g
|
|
Crest back steepness = 2*pi*Ac./Tu/Tcb/g
|
|
Total wave steepness (zero-downcrossing wave) = 2*pi*Hd./Td.^2/g
|
|
Total wave steepness (zero-upcrossing wave) = 2*pi*Hu./Tu.^2/g
|
|
|
|
The definition of g, Ac,At, Tcf, etc. are given in gravity and
|
|
wafo.definitions.
|
|
|
|
Example
|
|
-------
|
|
>>> import wafo.data as wd
|
|
>>> import wafo.objects as wo
|
|
>>> x = wd.sea()
|
|
>>> ts = wo.mat2timeseries(x)
|
|
>>> true_SH = [
|
|
... [[ 0.01186982, 0.04852534], [ 0.69, 0.86]],
|
|
... [[ 0.02918363, 0.06385979], [ 0.69, 0.86]],
|
|
... [[ 0.27797411, 0.33585743], [ 0.69, 0.86]],
|
|
... [[ 0.60835634, 0.60930197], [ 0.42, 0.78]],
|
|
... [[ 0.60835634, 0.60930197], [ 0.42, 0.78]],
|
|
... [[ 0.10140867, 0.06141156], [ 0.42, 0.78]],
|
|
... [[ 0.01821413, 0.01236672], [ 0.42, 0.78]]]
|
|
>>> for i in range(-3,4):
|
|
... S, H = ts.wave_height_steepness(kind=i)
|
|
... np.allclose((S[:2],H[:2]), true_SH[i+3])
|
|
True
|
|
True
|
|
True
|
|
True
|
|
True
|
|
True
|
|
True
|
|
|
|
import pylab as plt
|
|
h = plt.plot(S,H,'.')
|
|
h = plt.xlabel('S')
|
|
h = plt.ylabel('Hd [m]')
|
|
|
|
See also
|
|
--------
|
|
wafo.definitions
|
|
'''
|
|
|
|
dT = self.sampling_period() / np.maximum(rate, 1)
|
|
if g is None:
|
|
g = gravity() # acceleration of gravity
|
|
|
|
xi, ti = self._interpolate(rate)
|
|
|
|
tc_ind, z_ind = findtc(xi, v=0, kind='tw')
|
|
tc_a = xi[tc_ind]
|
|
tc_t = ti[tc_ind]
|
|
Ac = tc_a[1::2] # crest amplitude
|
|
At = -tc_a[0::2] # trough amplitude
|
|
defnr = dict(maxVcfVcb=0, Vcf=1, Vcb=-1, Scf=2, Scb=-2, StHd=3,
|
|
StHu=-3).get(kind, kind)
|
|
if 0 <= defnr <= 2:
|
|
# time between zero-upcrossing and crest [s]
|
|
tu = ecross(ti, xi, z_ind[1:-1:2], v=0)
|
|
Tcf = tc_t[1::2] - tu
|
|
Tcf[(Tcf == 0)] = dT # avoiding division by zero
|
|
if -2 <= defnr <= 0:
|
|
# time between crest and zero-downcrossing [s]
|
|
td = ecross(ti, xi, z_ind[2::2], v=0)
|
|
Tcb = td - tc_t[1::2]
|
|
Tcb[(Tcb == 0)] = dT
|
|
|
|
if defnr == 0:
|
|
# max(Vcf, Vcr) and the corresponding wave height Hd or Hu in H
|
|
Hu = Ac + At[1:]
|
|
Hd = Ac + At[:-1]
|
|
T = np.where(Tcf < Tcb, Tcf, Tcb)
|
|
S = Ac / T
|
|
H = np.where(Tcf < Tcb, Hd, Hu)
|
|
elif defnr == 1: # extracting crest front velocity [m/s] and
|
|
# Zero-downcrossing wave height [m]
|
|
H = Ac + At[:-1] # Hd
|
|
S = Ac / Tcf
|
|
elif defnr == -1: # extracting crest rear velocity [m/s] and
|
|
# Zero-upcrossing wave height [m]
|
|
H = Ac + At[1:] # Hu
|
|
S = Ac / Tcb
|
|
# crest front steepness in S and the wave height Hd in H.
|
|
elif defnr == 2:
|
|
H = Ac + At[:-1] # Hd
|
|
Td = diff(ecross(ti, xi, z_ind[::2], v=0))
|
|
S = 2 * pi * Ac / Td / Tcf / g
|
|
# crest back steepness in S and the wave height Hu in H.
|
|
elif defnr == -2:
|
|
H = Ac + At[1:]
|
|
Tu = diff(ecross(ti, xi, z_ind[1::2], v=0))
|
|
S = 2 * pi * Ac / Tu / Tcb / g
|
|
elif defnr == 3: # total steepness in S and the wave height Hd in H
|
|
# for zero-downcrossing waves.
|
|
H = Ac + At[:-1]
|
|
# Period zero-downcrossing waves
|
|
Td = diff(ecross(ti, xi, z_ind[::2], v=0))
|
|
S = 2 * pi * H / Td ** 2 / g
|
|
# total steepness in S and the wave height Hu in H for
|
|
elif defnr == -3:
|
|
# zero-upcrossing waves.
|
|
H = Ac + At[1:]
|
|
# Period zero-upcrossing waves
|
|
Tu = diff(ecross(ti, xi, z_ind[1::2], v=0))
|
|
S = 2 * pi * H / Tu ** 2 / g
|
|
|
|
return S, H
|
|
|
|
@staticmethod
|
|
def _default_index(x, vh, wdef, pdef):
|
|
if pdef in ('m2m', 'm2M', 'M2m', 'M2M'):
|
|
index = findtp(x, vh, wdef)
|
|
elif pdef in ('u2u', 'u2d', 'd2u', 'd2d'):
|
|
index = findcross(x, vh, wdef)
|
|
elif pdef in ('t2t', 't2c', 'c2t', 'c2c'):
|
|
index = findtc(x, vh, wdef)[0]
|
|
elif pdef in ('d2t', 't2u', 'u2c', 'c2d', 'all'):
|
|
index, v_ind = findtc(x, vh, wdef)
|
|
# sorting crossings and tp in sequence
|
|
index = sort(r_[index, v_ind])
|
|
else:
|
|
raise ValueError('Unknown pdef option! {}'.format(str(pdef)))
|
|
return index
|
|
|
|
def _get_start_index(self, pdef, down_crossing_or_max):
|
|
if down_crossing_or_max:
|
|
if pdef in ('d2t', 'M2m', 'c2t', 'd2u', 'M2M', 'c2c', 'd2d',
|
|
'all'):
|
|
start = 1
|
|
elif pdef in ('t2u', 'm2M', 't2c', 'u2d', 'm2m', 't2t', 'u2u'):
|
|
start = 2
|
|
elif pdef in ('u2c'):
|
|
start = 3
|
|
elif pdef in ('c2d'):
|
|
start = 4
|
|
else:
|
|
raise ValueError('Unknown pdef option!')
|
|
# else first is up-crossing or min
|
|
elif pdef in ('all', 'u2c', 'm2M', 't2c', 'u2d', 'm2m', 't2t', 'u2u'):
|
|
start = 0
|
|
elif pdef in ('c2d', 'M2m', 'c2t', 'd2u', 'M2M', 'c2c', 'd2d'):
|
|
start = 1
|
|
elif pdef in ('d2t'):
|
|
start = 2
|
|
elif pdef in ('t2u'):
|
|
start = 3
|
|
else:
|
|
raise ValueError('Unknown pdef option!')
|
|
return start
|
|
|
|
def _get_step(self, pdef):
|
|
# determine the steps between wanted periods
|
|
if pdef in ('d2t', 't2u', 'u2c', 'c2d'):
|
|
step = 4
|
|
elif pdef in ('all'):
|
|
step = 1 # secret option!
|
|
else:
|
|
step = 2
|
|
return step
|
|
|
|
def _interpolate(self, rate):
|
|
if rate > 1: # interpolate with spline
|
|
n = ceil(self.data.size * rate)
|
|
ti = linspace(self.args[0], self.args[-1], n)
|
|
x = stineman_interp(ti, self.args, self.data.ravel())
|
|
# xi = interp1d(self.args, self.data.ravel(), kind='cubic')(ti)
|
|
else:
|
|
x = self.data.ravel()
|
|
ti = self.args
|
|
return x, ti
|
|
|
|
def wave_periods(self, vh=None, pdef='d2d', wdef=None, index=None, rate=1):
|
|
"""
|
|
Return sequence of wave periods/lengths from data.
|
|
|
|
Parameters
|
|
----------
|
|
vh : scalar
|
|
reference level ( default v=mean(x(:,2)) ) or
|
|
rainflow filtering height (default h=0)
|
|
pdef : string
|
|
defining type of waveperiod (wavelength) returned:
|
|
Level v separated 't2c', 'c2t', 't2t' or 'c2c' -waveperiod.
|
|
Level v 'd2d', 'u2u', 'd2u' or 'u2d' -waveperiod.
|
|
Rain flow filtered (with height greater than h)
|
|
'm2M', 'M2m', 'm2m' or 'M2M' -waveperiod.
|
|
Explanation to the abbreviations:
|
|
M=Max, m=min, d=down-crossing, u=up-crossing ,
|
|
t=trough and c=crest.
|
|
Thus 'd2d' means period between a down-crossing to the
|
|
next down-crossing and 'u2c' means period between a
|
|
u-crossing to the following crest.
|
|
wdef : string
|
|
defining type of wave. Possible options are
|
|
'mw','Mw','dw', 'uw', 'tw', 'cw' or None.
|
|
If wdef is None all troughs and crests will be used,
|
|
otherwise only the troughs and crests which define a
|
|
wave according to the wavedefinition are used.
|
|
|
|
index : vector
|
|
index sequence of one of the following :
|
|
-level v-crossings (indices to "du" are required to
|
|
calculate 'd2d', 'd2u', 'u2d' or 'u2u' waveperiods)
|
|
-level v separated trough and crest turningpoints
|
|
(indices to 'tc' are required to calculate
|
|
't2t', 't2c', 'c2t' or 'c2c' waveperiods)
|
|
-level v crossings and level v separated trough and
|
|
crest turningpoints (indices to "dutc" are
|
|
required to calculate t2u, u2c, c2d or d2t
|
|
waveperiods)
|
|
-rainflow filtered turningpoints with minimum rfc height h
|
|
(indices to "mMtc" are required to calculate
|
|
'm2m', 'm2M', 'M2m' or 'M2M' waveperiods)
|
|
|
|
rate : scalar
|
|
interpolation rate. If rate larger than one, then x is
|
|
interpolated before extrating T
|
|
|
|
Returns
|
|
--------
|
|
T : vector
|
|
sequence of waveperiods (or wavelengths).
|
|
index : vector
|
|
of indices
|
|
|
|
|
|
Example:
|
|
--------
|
|
Histogram of crest2crest waveperiods
|
|
>>> import wafo.data as wd
|
|
>>> import wafo.objects as wo
|
|
>>> import pylab as plb
|
|
>>> x = wd.sea()
|
|
>>> ts = wo.mat2timeseries(x[0:400,:])
|
|
>>> T, ix = ts.wave_periods(vh=0.0, pdef='c2c')
|
|
>>> np.allclose(T[:3], [-0.27, -0.08, 0.32])
|
|
True
|
|
|
|
h = plb.hist(T)
|
|
|
|
See also:
|
|
--------
|
|
findtp,
|
|
findtc,
|
|
findcross, perioddef
|
|
"""
|
|
|
|
x, ti = self._interpolate(rate)
|
|
|
|
if vh is None:
|
|
if pdef[0] in ('m', 'M'):
|
|
vh = 0
|
|
print(' The minimum rfc height, h, is set to: %g' % vh)
|
|
else:
|
|
vh = x.mean()
|
|
print(' The level l is set to: %g' % vh)
|
|
|
|
if index is None:
|
|
index = self._default_index(x, vh, wdef, pdef)
|
|
|
|
down_crossing_or_max = (x[index[0]] > x[index[1]])
|
|
start = self._get_start_index(pdef, down_crossing_or_max)
|
|
step = self._get_step(pdef)
|
|
|
|
# determine the distance between min2min, t2t etc..
|
|
if pdef in ('m2m', 't2t', 'u2u', 'M2M', 'c2c', 'd2d'):
|
|
dist = 2
|
|
else:
|
|
dist = 1
|
|
|
|
nn = len(index)
|
|
|
|
if pdef[0] in ('u', 'd'):
|
|
t0 = ecross(ti, x, index[start:(nn - dist):step], vh)
|
|
else: # min, Max, trough, crest or all crossings wanted
|
|
t0 = x[index[start:(nn - dist):step]]
|
|
|
|
if pdef[2] in ('u', 'd'):
|
|
t1 = ecross(ti, x, index[(start + dist):nn:step], vh)
|
|
else: # min, Max, trough, crest or all crossings wanted
|
|
t1 = x[index[(start + dist):nn:step]]
|
|
|
|
T = t1 - t0
|
|
return T, index
|
|
|
|
def reconstruct(self, inds=None, Nsim=20, L=None, def_='nonlinear',
|
|
**options):
|
|
'''
|
|
function [y,g,g2,test,tobs,mu1o, mu1oStd] = reconstruct(x,)
|
|
RECONSTRUCT reconstruct the spurious/missing points of timeseries
|
|
|
|
CALL: [y,g,g2,test,tobs,mu1o,mu1oStd]=
|
|
reconstruct(x,inds,Nsim,L,def,options)
|
|
|
|
Returns
|
|
-------
|
|
y = reconstructed signal
|
|
g,g2 = smoothed and empirical transformation, respectively
|
|
test, tobs = test observator int(g(u)-u)^2 du and
|
|
int(g_new(u)-g_old(u))^2 du,
|
|
respectively, where int limits is given by param in lc2tr.
|
|
Test is a measure of departure from the Gaussian model for
|
|
the data. Tobs is a measure of the convergence of the
|
|
estimation of g.
|
|
mu1o = expected surface elevation of the Gaussian model process.
|
|
mu1o_std = standarddeviation of mu1o.
|
|
|
|
Parameters
|
|
----------
|
|
x : 2 column timeseries
|
|
first column sampling times [sec]
|
|
second column surface elevation [m]
|
|
inds : integer array
|
|
indices to spurious points of x
|
|
Nsim = the maximum # of iterations before we stop
|
|
|
|
L = lag size of the Parzen window function.
|
|
If no value is given the lag size is set to
|
|
be the lag where the auto correlation is less than
|
|
2 standard deviations. (maximum 200)
|
|
def :
|
|
'nonlinear' : transform from smoothed crossing intensity (default)
|
|
'mnonlinear': transform from smoothed marginal distribution
|
|
'linear' : identity.
|
|
options = options structure defining how the estimation of g is
|
|
done, see troptset.
|
|
|
|
In order to reconstruct the data a transformed Gaussian random process
|
|
is used for modelling and simulation of the missing/removed data
|
|
conditioned on the other known observations.
|
|
|
|
Estimates of standarddeviations of y is obtained by a call to tranproc
|
|
Std = tranproc(mu1o+/-mu1oStd,fliplr(g));
|
|
|
|
See also
|
|
--------
|
|
troptset, findoutliers, cov2csdat, dat2cov, dat2tr, detrendma
|
|
|
|
Reference
|
|
---------
|
|
Brodtkorb, P, Myrhaug, D, and Rue, H (2001)
|
|
"Joint distribution of wave height and wave crest velocity from
|
|
reconstructed data with application to ringing"
|
|
Int. Journal of Offshore and Polar Engineering, Vol 11, No. 1,
|
|
pp 23--32
|
|
|
|
Brodtkorb, P, Myrhaug, D, and Rue, H (1999)
|
|
"Joint distribution of wave height and wave crest velocity from
|
|
reconstructed data
|
|
in Proceedings of 9th ISOPE Conference, Vol III, pp 66-73
|
|
'''
|
|
|
|
opt = DotDict(chkder=True, plotflag=False, csm=0.9, gsm=.05,
|
|
param=(-5, 5, 513), delay=2, linextrap=True, ntr=10000,
|
|
ne=7, gvar=1)
|
|
opt.update(options)
|
|
|
|
_xn = self.data.copy().ravel()
|
|
# n = len(xn)
|
|
#
|
|
# if n < 2:
|
|
# raise ValueError('The vector must have more than 2 elements!')
|
|
#
|
|
# param = opt.param
|
|
# plotflags = dict(none=0, off=0, final=1, iter=2)
|
|
# plotflag = plotflags.get(opt.plotflag, opt.plotflag)
|
|
#
|
|
# olddef = def_
|
|
# method = 'approx'
|
|
# ptime = opt.delay # pause for ptime sec if plotflag=2
|
|
#
|
|
# expect1 = 1 # first reconstruction 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)
|
|
#
|
|
# cmvmax = 100
|
|
# # if number of consecutive missing values (cmv) are longer they
|
|
# # are not used in estimation of g, due to the fact that the
|
|
# # conditional expectation approaches zero as the length to
|
|
# # the closest known points increases, see below in the for loop
|
|
# dT = self.sampling_period()
|
|
#
|
|
# Lm = np.minimum([n, 200, int(200/dT)]) # Lagmax 200 seconds
|
|
# if L is not None:
|
|
# Lm = max(L, Lm)
|
|
# # Lma: size of the moving average window used for detrending the
|
|
# # reconstructed signal
|
|
# Lma = 1500
|
|
# if inds is not None:
|
|
# xn[inds] = np.nan
|
|
#
|
|
# inds = isnan(xn)
|
|
# if not inds.any():
|
|
# raise ValueError('No spurious data given')
|
|
#
|
|
# endpos = np.diff(inds)
|
|
# strtpos = np.flatnonzero(endpos > 0)
|
|
# endpos = np.flatnonzero(endpos < 0)
|
|
#
|
|
# indg = np.flatnonzero(1-inds) # indices to good points
|
|
# inds = np.flatnonzero(inds) # indices to spurious points
|
|
#
|
|
# indNaN = [] # indices to points omitted in the covariance estimation
|
|
# indr = np.arange(n) # indices to point used in the estimation of g
|
|
#
|
|
# # Finding more than cmvmax consecutive spurios points.
|
|
# # They will not be used in the estimation of g and are thus removed
|
|
# # from indr.
|
|
#
|
|
# if strtpos.size > 0 and (endpos.size == 0 or
|
|
# endpos[-1] < strtpos[-1]):
|
|
# if (n - strtpos[-1]) > cmvmax:
|
|
# indNaN = indr[strtpos[-1]+1:n]
|
|
# indr = indr[:strtpos[-1]+1]
|
|
# strtpos = strtpos[:-1]
|
|
#
|
|
# if endpos.size > 0 and (strtpos.size == 0 or endpos[0] < strtpos[0]):
|
|
# if endpos[0] > cmvmax:
|
|
# indNaN = np.hstack((indNaN, indr[:endpos[0]]))
|
|
# indr = indr[endpos[0]:]
|
|
#
|
|
# strtpos = strtpos-endpos[0]
|
|
# endpos = endpos-endpos[0]
|
|
# endpos = endpos[1:]
|
|
#
|
|
# for ix in range(len(strtpos)-1, -1, -1):
|
|
# if (endpos[ix]-strtpos[ix] > cmvmax):
|
|
# indNaN = np.hstack((indNaN, indr[strtpos[ix]+1:endpos[ix]]))
|
|
# # remove this when estimating the transform
|
|
# del indr[strtpos[ix]+1:endpos[ix]]
|
|
#
|
|
# if len(indr) < 0.1*n:
|
|
# raise ValueError('Not possible to reconstruct signal')
|
|
#
|
|
# if indNaN.any():
|
|
# indNaN = np.sort(indNaN)
|
|
#
|
|
# # initial reconstruction attempt
|
|
# xn[indg, 1] = detrendma(xn[indg, 1], 1500)
|
|
# g, test, cmax, irr, g2 = dat2tr(xn[indg, :], def_, opt)
|
|
# xnt = xn.copy()
|
|
# xnt[indg,:] = dat2gaus(xn[indg,:], g)
|
|
# xnt[inds, 1] = np.nan
|
|
# rwin = findrwin(xnt, Lm, L)
|
|
# print('First reconstruction attempt, e(g-u) = {}'.format(test))
|
|
# # old simcgauss
|
|
# [samp ,mu1o, mu1oStd] = cov2csdat(xnt(:,2),rwin,1,method,inds);
|
|
# if expect1,# reconstruction by expectation
|
|
# xnt(inds,2) =mu1o;
|
|
# else
|
|
# xnt(inds,2) =samp;
|
|
# end
|
|
# xn=gaus2dat(xnt,g);
|
|
# xn(:,2)=detrendma(xn(:,2),Lma); # detrends the signal with a moving
|
|
# # average of size Lma
|
|
# g_old=g;
|
|
#
|
|
# bias = mean(xn(:,2));
|
|
# xn(:,2)=xn(:,2)-bias; # bias correction
|
|
#
|
|
# if plotflag==2
|
|
# clf
|
|
# mind=1:min(1500,n);
|
|
# waveplot(xn(mind,:),x(inds(mind),:), 6,1)
|
|
# subplot(111)
|
|
# pause(ptime)
|
|
# end
|
|
#
|
|
# test0=0;
|
|
# for ix=1:Nsim,
|
|
# # if 0,#ix==2,
|
|
# # rwin=findrwin(xn,Lm,L);
|
|
# # xs=cov2sdat(rwin,[n 100 dT]);
|
|
# # [g0 test0 cmax irr g2] = dat2tr(xs,def,opt);
|
|
# # [test0 ind0]=sort(test0);
|
|
# # end
|
|
|
|
# if 1, #test>test0(end-5),
|
|
# # 95# sure the data comes from a non-Gaussian process
|
|
# def = olddef; #Non Gaussian process
|
|
# else
|
|
# def = 'linear'; # Gaussian process
|
|
# end
|
|
# # used for isope article
|
|
# # indr =[1:27000 30000:39000];
|
|
# # Too many consecutive missing values will influence the
|
|
# # estimation of g. By default do not use consecutive missing
|
|
# # values if there are more than cmvmax.
|
|
#
|
|
# [g test cmax irr g2] = dat2tr(xn(indr,:),def,opt);
|
|
# if plotflag==2,
|
|
# pause(ptime)
|
|
# end
|
|
#
|
|
# #tobs=sqrt((param(2)-param(1))/(param(3)-1)*
|
|
# sum((g_old(:,2)-g(:,2)).^2))
|
|
# # new call
|
|
# tobs=sqrt((param(2)-param(1))/(param(3)-1)
|
|
# *sum((g(:,2)-interp1(g_old(:,1)-bias, g_old(:,2),g(:,1),
|
|
# 'spline')).^2));
|
|
#
|
|
# if ix>1
|
|
# if tol>tobs2 && tol>tobs,
|
|
# break, #estimation of g converged break out of for loop
|
|
# end
|
|
# end
|
|
#
|
|
# tobs2=tobs;
|
|
#
|
|
# xnt=dat2gaus(xn,g);
|
|
# if ~isempty(indNaN), xnt(indNaN,2)=NaN; end
|
|
# rwin=findrwin(xnt,Lm,L);
|
|
# disp(['Simulation nr: ', int2str(ix), ' of ' num2str(Nsim),
|
|
# ' e(g-g_old)=', num2str(tobs), ', e(g-u)=', num2str(test)])
|
|
# [samp ,mu1o, mu1oStd] = cov2csdat(xnt(:,2),rwin,1,method,inds);
|
|
#
|
|
# if expect,
|
|
# xnt(inds,2) =mu1o;
|
|
# else
|
|
# xnt(inds,2) =samp;
|
|
# end
|
|
#
|
|
# xn=gaus2dat(xnt,g);
|
|
# if ix<Nsim
|
|
# bias=mean(xn(:,2));
|
|
# xn(:,2) = (xn(:,2)-bias); # bias correction
|
|
# end
|
|
# g_old=g;# saving the last transform
|
|
# if plotflag==2
|
|
# waveplot(xn(mind,:),x(inds(mind),:),6,1,[])
|
|
# subplot(111)
|
|
# pause(ptime)
|
|
# end
|
|
# end # for loop
|
|
#
|
|
# if 1, #test>test0(end-5)
|
|
# xnt=dat2gaus(xn,g);
|
|
# [samp ,mu1o, mu1oStd] = cov2csdat(xnt(:,2),rwin,1,method,inds);
|
|
# xnt(inds,2) =samp;
|
|
# xn=gaus2dat(xnt,g);
|
|
# bias=mean(xn(:,2));
|
|
# xn(:,2) = (xn(:,2)-bias); # bias correction
|
|
# g(:,1)=g(:,1)-bias;
|
|
# g2(:,1)=g2(:,1)-bias;
|
|
# gn=trangood(g);
|
|
#
|
|
# #mu1o=mu1o-tranproc(bias,gn);
|
|
# muUStd=tranproc(mu1o+2*mu1oStd,fliplr(gn));#
|
|
# muLStd=tranproc(mu1o-2*mu1oStd,fliplr(gn));#
|
|
# else
|
|
# muLStd=mu1o-2*mu1oStd;
|
|
# muUStd=mu1o+2*mu1oStd;
|
|
# end
|
|
#
|
|
# if plotflag==2 && length(xn)<10000,
|
|
# waveplot(xn,[xn(inds,1) muLStd ;xn(inds,1) muUStd ],
|
|
# 6,round(n/3000),[])
|
|
# legend('reconstructed','2 stdev')
|
|
# #axis([770 850 -1 1])
|
|
# #axis([1300 1325 -1 1])
|
|
# end
|
|
# y=xn;
|
|
# toc
|
|
#
|
|
# return
|
|
#
|
|
# def findrwin(xnt, Lm, L=None):
|
|
# r = dat2cov(xnt, Lm) # computes ACF
|
|
# # finding where ACF is less than 2 st. deviations .
|
|
# # in order to find a better L value
|
|
# if L is None:
|
|
# L = np.flatnonzero(np.abs(r.R) > 2 * r.stdev)
|
|
# if len(L) == 0:
|
|
# L = Lm;
|
|
# else:
|
|
# L = min(np.floor(4/3*(L[-1] + 1), Lm)
|
|
# win = parzen(2 * L - 1)
|
|
# r.R[:L] = win[L:2*L-1] * r.R[:L]
|
|
# r.R[L:] = 0
|
|
# return r
|
|
|
|
def plot_wave(self, sym1='k.', ts=None, sym2='k+', nfig=None, nsub=None,
|
|
sigma=None, vfact=3):
|
|
'''
|
|
Plots the surface elevation of timeseries.
|
|
|
|
Parameters
|
|
----------
|
|
sym1, sym2 : string
|
|
plot symbol and color for data and ts, respectively
|
|
(see PLOT) (default 'k.' and 'k+')
|
|
ts : TimeSeries or TurningPoints object
|
|
to overplot data. default zero-separated troughs and crests.
|
|
nsub : scalar integer
|
|
Number of subplots in each figure. By default nsub is such that
|
|
there are about 20 mean down crossing waves in each subplot.
|
|
If nfig is not given and nsub is larger than 6 then nsub is
|
|
changed to nsub=min(6,ceil(nsub/nfig))
|
|
nfig : scalar integer
|
|
Number of figures. By default nfig=ceil(Nsub/6).
|
|
sigma : real scalar
|
|
standard deviation of data.
|
|
vfact : real scalar
|
|
how large in stdev the vertical scale should be (default 3)
|
|
|
|
|
|
Example
|
|
-------
|
|
Plot x1 with red lines and mark troughs and crests with blue circles.
|
|
>>> import wafo
|
|
>>> x = wafo.data.sea()
|
|
>>> ts150 = wafo.objects.mat2timeseries(x[:150,:])
|
|
|
|
>>> h = ts150.plot_wave('r-', sym2='bo')
|
|
|
|
See also
|
|
--------
|
|
findtc, plot
|
|
'''
|
|
|
|
nw = 20
|
|
tn = self.args
|
|
xn = self.data.ravel()
|
|
indmiss = isnan(xn) # indices to missing points
|
|
indg = where(1 - indmiss)[0]
|
|
if ts is None:
|
|
tc_ix = findtc(xn[indg], 0, 'tw')[0]
|
|
xn2 = xn[tc_ix]
|
|
tn2 = tn[tc_ix]
|
|
else:
|
|
xn2 = ts.data
|
|
tn2 = ts.args
|
|
|
|
if sigma is None:
|
|
sigma = xn[indg].std()
|
|
|
|
if nsub is None:
|
|
# about Nw mdc waves in each plot
|
|
nsub = int(len(xn2) / (2 * nw)) + 1
|
|
if nfig is None:
|
|
nfig = int(ceil(nsub / 6))
|
|
nsub = min(6, int(ceil(nsub / nfig)))
|
|
|
|
n = len(xn)
|
|
Ns = int(n / (nfig * nsub))
|
|
ind = r_[0:Ns]
|
|
if all(xn >= 0):
|
|
vscale = [0, 2 * sigma * vfact] # @UnusedVariable
|
|
else:
|
|
vscale = array([-1, 1]) * vfact * sigma # @UnusedVariable
|
|
|
|
XlblTxt = 'Time [sec]'
|
|
dT = 1
|
|
timespan = tn[ind[-1]] - tn[ind[0]]
|
|
if abs(timespan) > 18000: # more than 5 hours
|
|
dT = 1 / (60 * 60)
|
|
XlblTxt = 'Time (hours)'
|
|
elif abs(timespan) > 300: # more than 5 minutes
|
|
dT = 1 / 60
|
|
XlblTxt = 'Time (minutes)'
|
|
|
|
if np.max(abs(xn[indg])) > 5 * sigma:
|
|
XlblTxt = XlblTxt + ' (Spurious data since max > 5 std.)'
|
|
|
|
plot = plt.plot
|
|
subplot = plt.subplot
|
|
figs = []
|
|
for unused_iz in range(nfig):
|
|
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)
|
|
h_scale = array([tn[ind[0]], tn[ind[-1]]])
|
|
ind2 = where((h_scale[0] <= tn2) & (tn2 <= h_scale[1]))[0]
|
|
plot(tn[ind] * dT, xn[ind], sym1)
|
|
if len(ind2) > 0:
|
|
plot(tn2[ind2] * dT, xn2[ind2], sym2)
|
|
plot(h_scale * dT, [0, 0], 'k-')
|
|
# plt.axis([h_scale*dT, v_scale])
|
|
for iy in [-2, 2]:
|
|
plot(h_scale * dT, iy * sigma * ones(2), ':')
|
|
ind = ind + Ns
|
|
plt.xlabel(XlblTxt)
|
|
|
|
return figs
|
|
|
|
def plot_sp_wave(self, wave_idx_, *args, **kwds):
|
|
"""
|
|
Plot specified wave(s) from timeseries
|
|
|
|
Parameters
|
|
----------
|
|
wave_idx : integer vector
|
|
of indices to waves we want to plot, i.e., wave numbers.
|
|
tz_idx : integer vector
|
|
of indices to the beginning, middle and end of
|
|
defining wave, i.e. for zero-downcrossing waves, indices to
|
|
zerocrossings (default trough2trough wave)
|
|
|
|
Examples
|
|
--------
|
|
Plot waves nr. 6,7,8 and waves nr. 12,13,...,17
|
|
>>> import wafo
|
|
>>> x = wafo.data.sea()
|
|
>>> ts = wafo.objects.mat2timeseries(x[0:500,...])
|
|
|
|
>>> h = ts.plot_sp_wave(np.r_[6:9,12:18])
|
|
|
|
See also
|
|
--------
|
|
plot_wave, findtc
|
|
"""
|
|
wave_idx = atleast_1d(wave_idx_).flatten()
|
|
tz_idx = kwds.pop('tz_idx', None)
|
|
if tz_idx is None:
|
|
# finding trough to trough waves
|
|
unused_tc_ind, tz_idx = findtc(self.data, 0, 'tw')
|
|
|
|
dw = nonzero(abs(diff(wave_idx)) > 1)[0]
|
|
Nsub = dw.size + 1
|
|
Nwp = zeros(Nsub, dtype=int)
|
|
if Nsub > 1:
|
|
dw = dw + 1
|
|
Nwp[Nsub - 1] = wave_idx[-1] - wave_idx[dw[-1]] + 1
|
|
wave_idx[dw[-1] + 1:] = -2
|
|
for ix in range(Nsub - 2, 1, -2):
|
|
# of waves pr subplot
|
|
Nwp[ix] = wave_idx[dw[ix] - 1] - wave_idx[dw[ix - 1]] + 1
|
|
wave_idx[dw[ix - 1] + 1:dw[ix]] = -2
|
|
|
|
Nwp[0] = wave_idx[dw[0] - 1] - wave_idx[0] + 1
|
|
wave_idx[1:dw[0]] = -2
|
|
wave_idx = wave_idx[wave_idx > -1]
|
|
else:
|
|
Nwp[0] = wave_idx[-1] - wave_idx[0] + 1
|
|
|
|
Nsub = min(6, Nsub)
|
|
Nfig = int(ceil(Nsub / 6))
|
|
Nsub = min(6, int(ceil(Nsub / Nfig)))
|
|
figs = []
|
|
for unused_iy in range(Nfig):
|
|
figs.append(plt.figure())
|
|
for ix in range(Nsub):
|
|
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
|
|
plt.plot(self.args[ind], self.data[ind], *args, **kwds)
|
|
plt.hold('on')
|
|
xi = [self.args[ind[0]], self.args[ind[-1]]]
|
|
plt.plot(xi, [0, 0])
|
|
|
|
if Nwp[ix] == 1:
|
|
plt.ylabel('Wave %d' % wave_idx[ix])
|
|
else:
|
|
plt.ylabel(
|
|
'Wave %d - %d' % (wave_idx[ix],
|
|
wave_idx[ix] + Nwp[ix] - 1))
|
|
plt.xlabel('Time [sec]')
|
|
# wafostamp
|
|
return figs
|
|
|
|
|
|
def test_docstrings():
|
|
import doctest
|
|
print('Testing docstrings in %s' % __file__)
|
|
doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
|
|
|
|
|
|
if __name__ == '__main__':
|
|
test_docstrings()
|