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

Cleanup of code

Browse Source
master
Per A Brodtkorb 9 years ago
parent e6ab39cbe7
commit 9522c6c24f
15 changed files with 831 additions and 1163 deletions
Show Stats Download Patch File Download Diff File

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

@ -14,7 +14,7 @@ note : Memorandum string.
date : Date and time of creation or change. date : Date and time of creation or change.
''' '''
from __future__ import division from __future__ import division, absolute_import
import warnings import warnings
import numpy as np import numpy as np
from numpy import (zeros, ones, sqrt, inf, where, nan, from numpy import (zeros, ones, sqrt, inf, where, nan,
@ -27,9 +27,9 @@ from scipy.linalg import toeplitz, lstsq
from scipy import sparse from scipy import sparse
from pylab import stineman_interp from pylab import stineman_interp
from wafo.containers import PlotData from ..containers import PlotData
from wafo.misc import sub_dict_select, nextpow2 # , JITImport from ..misc import sub_dict_select, nextpow2 # , JITImport
import wafo.spectrum as _wafospec from .. import spectrum as _wafospec
from scipy.sparse.linalg.dsolve.linsolve import spsolve from scipy.sparse.linalg.dsolve.linsolve import spsolve
from scipy.sparse.base import issparse from scipy.sparse.base import issparse
from scipy.signal.windows import parzen from scipy.signal.windows import parzen

1
wafo/covariance/tests/__init__.py
Unescape Escape View File

@ -1 +0,0 @@

12
wafo/gaussian.py
Unescape Escape View File

@ -355,8 +355,8 @@ class Rind(object):
dev = sqrt(diag(BIG)) # std dev = sqrt(diag(BIG)) # std
ind = nonzero(indI[1:] > -1)[0] ind = nonzero(indI[1:] > -1)[0]
infin = repeat(2, len(indI) - 1) infin = repeat(2, len(indI) - 1)
infin[ind] = (2 - (Bup[0, ind] > infinity * dev[indI[ind + 1]]) infin[ind] = (2 - (Bup[0, ind] > infinity * dev[indI[ind + 1]]) -
- 2 * (Blo[0, ind] < -infinity * dev[indI[ind + 1]])) 2 * (Blo[0, ind] < -infinity * dev[indI[ind + 1]]))
Bup[0, ind] = minimum(Bup[0, ind], infinity * dev[indI[ind + 1]]) Bup[0, ind] = minimum(Bup[0, ind], infinity * dev[indI[ind + 1]])
Blo[0, ind] = maximum(Blo[0, ind], -infinity * dev[indI[ind + 1]]) Blo[0, ind] = maximum(Blo[0, ind], -infinity * dev[indI[ind + 1]])
@ -992,10 +992,10 @@ def prbnorm2d(a, b, r):
infin = np.repeat(2, 2) - (upper > infinity) - 2 * (lower < -infinity) infin = np.repeat(2, 2) - (upper > infinity) - 2 * (lower < -infinity)
if np.all(infin == 2): if np.all(infin == 2):
return (bvd(lower[0], lower[1], correl) return (bvd(lower[0], lower[1], correl) -
- bvd(upper[0], lower[1], correl) bvd(upper[0], lower[1], correl) -
- bvd(lower[0], upper[1], correl) bvd(lower[0], upper[1], correl) +
+ bvd(upper[0], upper[1], correl)) bvd(upper[0], upper[1], correl))
elif (infin[0] == 2 and infin[1] == 1): elif (infin[0] == 2 and infin[1] == 1):
return (bvd(lower[0], lower[1], correl) - return (bvd(lower[0], lower[1], correl) -
bvd(upper[0], lower[1], correl)) bvd(upper[0], lower[1], correl))

226
wafo/integrate.py
Unescape Escape View File

@ -266,8 +266,8 @@ def romberg(fun, a, b, releps=1e-3, abseps=1e-3):
fp[i] = 4 * fp[i - 1] fp[i] = 4 * fp[i - 1]
# Richardson extrapolation # Richardson extrapolation
for k in range(i): for k in range(i):
rom[two, k + 1] = rom[two, k] + \ rom[two, k + 1] = (rom[two, k] +
(rom[two, k] - rom[one, k]) / (fp[k] - 1) (rom[two, k] - rom[one, k]) / (fp[k] - 1))
Ih1 = Ih2 Ih1 = Ih2
Ih2 = Ih4 Ih2 = Ih4
@ -1119,6 +1119,8 @@ def quadgr(fun, a, b, abseps=1e-5, max_iter=17):
''' '''
# Author: jonas.lundgren@saabgroup.com, 2009. license BSD # Author: jonas.lundgren@saabgroup.com, 2009. license BSD
# Order limits (required if infinite limits) # Order limits (required if infinite limits)
a = np.asarray(a)
b = np.asarray(b)
if a == b: if a == b:
Q = b - a Q = b - a
err = b - a err = b - a
@ -1138,17 +1140,17 @@ def quadgr(fun, a, b, abseps=1e-5, max_iter=17):
# Change of variable # Change of variable
if np.isfinite(a) & np.isinf(b): if np.isfinite(a) & np.isinf(b):
# a to inf # a to inf
fun1 = lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2 [Q, err] = quadgr(lambda t: fun(a + t / (1 - t)) / (1 - t) ** 2,
[Q, err] = quadgr(fun1, 0, 1, abseps) 0, 1, abseps)
elif np.isinf(a) & np.isfinite(b): elif np.isinf(a) & np.isfinite(b):
# -inf to b # -inf to b
fun2 = lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2 [Q, err] = quadgr(lambda t: fun(b + t / (1 + t)) / (1 + t) ** 2,
[Q, err] = quadgr(fun2, -1, 0, abseps) -1, 0, abseps)
else: # -inf to inf else: # -inf to inf
fun1 = lambda t: fun(t / (1 - t)) / (1 - t) ** 2 [Q1, err1] = quadgr(lambda t: fun(t / (1 - t)) / (1 - t) ** 2,
fun2 = lambda t: fun(t / (1 + t)) / (1 + t) ** 2 0, 1, abseps / 2)
[Q1, err1] = quadgr(fun1, 0, 1, abseps / 2) [Q2, err2] = quadgr(lambda t: fun(t / (1 + t)) / (1 + t) ** 2,
[Q2, err2] = quadgr(fun2, -1, 0, abseps / 2) -1, 0, abseps / 2)
Q = Q1 + Q2 Q = Q1 + Q2
err = err1 + err2 err = err1 + err2
@ -1187,8 +1189,8 @@ def quadgr(fun, a, b, abseps=1e-5, max_iter=17):
hh = hh / 2 hh = hh / 2
x = np.hstack([x + a, x + b]) / 2 x = np.hstack([x + a, x + b]) / 2
# Quadrature # Quadrature
Q0[k] = hh * \ Q0[k] = hh * np.sum(wq * np.sum(np.reshape(fun(x), (-1, nq)), axis=0),
np.sum(wq * np.sum(np.reshape(fun(x), (-1, nq)), axis=0), axis=0) axis=0)
# Richardson extrapolation # Richardson extrapolation
if k >= 5: if k >= 5:
@ -1424,207 +1426,7 @@ def test_docstrings():
doctest.testmod() doctest.testmod()
# def levin_integrate():
# ''' An oscillatory integral
# Sheehan Olver, December 2010
#
#
# (Chebfun example quad/LevinIntegrate.m)
#
# This example computes the highly oscillatory integral of
#
# f * exp( 1i * w * g ),
#
# over (0,1) using the Levin method [1]. This method computes the integral
# by rewriting it as an ODE
#
# u' + 1i * w * g' u = f,
#
# so that the indefinite integral of f * exp( 1i * w * g ) is
#
# u * exp( 1i * w * g ).
#
#
#
# We use as an example
#
# f = 1 / ( x + 2 );
# g = cos( x - 2 );
# w = 100000;
#
# #
# References:
#
# [1] Levin, D., Procedures for computing one and two-dimensional integrals
# of functions with rapid irregular oscillations, Maths Comp., 38 (1982) 531--538
# '''
# exp = np.exp
# domain=[0, 1]
# x = Chebfun.identity(domain=domain)
# f = 1./(x+2)
# g = np.cos(x-2)
# D = np.diff(domain)
#
#
# # Here is are plots of this integrand, with w = 100, in complex space
# w = 100;
# line_opts = dict(line_width=1.6)
# font_opts = dict(font_size= 14)
# #
#
# intg = f*exp(1j*w*g)
# xs, ys, xi, yi, d = intg.plot_data(1000)
# #intg.plot(with_interpolation_points=True)
# #xi = np.linspace(0, 1, 1024)
# # plt.plot(xs, ys) # , **line_opts)
# # plt.plot(xi, yi, 'r.')
# # #axis equal
# # plt.title('Complex plot of integrand') #,**font_opts)
# # plt.show('hold')
# ##
# # and of just the real part
# # intgr = np.real(intg)
# # xs, ys, xi, yi, d = intgr.plot_data(1000)
# #intgr.plot()
# # plt.plot(xs, np.real(ys)) # , **line_opts)
# # plt.plot(xi, np.real(yi), 'r.')
# #axis equal
# # plt.title('Real part of integrand') #,**font_opts)
# # plt.show('hold')
#
# ##
# # The Levin method will be accurate for large and small w, and the time
# # taken is independent of w. Here we take a reasonably large value of w.
# w = 1000;
# intg = f*exp(1j*w*g)
# val0 = np.sum(intg)
# # val1 = sum(intg)
# print(val0)
# ##
# # Start timing
# #tic
#
# ##
# # Construct the operator L
# L = D + 1j*w*np.diag(g.differentiate())
#
# ##
# # From asymptotic analysis, we know that there exists a solution to the
# # equation which is non-oscillatory, though we do not know what initial
# # condition it satisfies. Thus we find a particular solution to this
# # equation with no boundary conditions.
#
# u = L / f
#
# ##
# # Because L is a differential operator with derivative order 1, \ expects
# # it to be given a boundary condition, which is why the warning message is
# # displayed. However, this doesn't cause any problems: though there are,
# # in fact, a family of solutions to the ODE without boundary conditions
# # due to the kernel
# #
# # exp(- 1i * w * g),
# #
# # it does not actually matter which particular solution is computed.
# # Non-uniqueness is also not an issue: \ in matlab is least squares, hence
# # does not require uniqueness. The existence of a non-oscillatory solution
# # ensures that \ converges to a u with length independent of w.
# #
# # One could prevent the warning by applying a boundary condition consistent
# # with the rest of the system, that is
# # L.lbc = {L(1,:),f(0)};
#
# ##
# # Now we evaluate the antiderivative at the endpoints to obtain the
# # integral.
#
# u(1)*exp(1j*w*g(1)) - u(0)*exp(1j*w*g(0))
#
# #toc
#
#
# ##
# # Here is a way to compute the integral using Clenshaw--Curtis quadrature.
# # As w becomes large, this takes an increasingly long time as the
# # oscillations must be resolved.
#
# #tic
# sum( f*exp(1j*w*g) )
# #toc
# aLevinTQ[omega_,a_,b_,f_,g_,nu_,wprec_,prm_,test_,basis_,np_]:=
#
# Module[{r,betam,A,AA,BB,S,F,w=N[omega, wprec]},
# M=Length[nu]-1;
# PB[k_,t_]:=If[basis==1,t^k,ChebyshevT[k,t]];
#
# ff[t_]:=((b-a)/2)*f[(b-a)*t/2+(a+b)/2];
#
# gg[t_]:=g[(b-a)*t/2+(a+b)/2];
# dgg[t_]:=Derivative[1][gg][t];
#
# If[test==0, betam=Min[Abs[dgg[-1]*w], Abs[dgg[1]*w]];
# While[prm*M/betam >=1, betam=2*betam]];
# If[test>0,x[k_]:=N[Cos[k*Pi/M], wprec],x[k_]:=
# Which[k<prm*M, N[-1+k/betam, wprec], k==Ceiling[prm*M],0,
# k>prm*M, N[1-(M-k)/betam, wprec]]];
#
# Psi[k_,t_]:=Derivative[0,1][PB][k,t]+I*w*dgg[t]*PB[k,t];
#
# ne[j_]:=nu[[j+1]]; S[-1]=0; S[j_]:=Sum[ne[i],{i,0,j}];
# nn=S[M]-1;
# A=ConstantArray[0,{nn+1,nn+1}];
# F=ConstantArray[0,nn+1]; r=0;
# While[r<M+1, Do[Do[ AA[j,k]=
# Limit[Derivative[0,Mod[j-S[r-1],ne[r]]][Psi][k,t],t->x[r]],
# {k,0,S[M]-1}],{j,S[r-1],S[r]-1}];
#
# Do[BB[j]=Limit[Derivative[Mod[j-S[r-1],ne[r]]][ff][t],
# t->x[r]],{j,S[r-1],S[r]-1}];
# Do[F[[j]]=BB[j-1],{j,S[r-1]+1,S[r]}];
# Do[Do[A[[j,k]]=AA[j-1,k-1],{k,1,S[M]}],{j,S[r-1]+1,S[r]}];
# r=r+1;]; (*sv=SingularValueList[N[A,wprec]];
# con=sv[[1]]/sv[[-1]]; Print["cond2(A)= ",N[con,3]];*)
# LS=Block[{MaxExtraPrecision=0},LeastSquares[N[A, wprec],F]];
# vvv[t_]:=Sum[LS[[k+1]]*PB[k,t], {k,0,nn}];
# NR=vvv[1]*Exp[I*w*gg[1]]-vvv[-1]*Exp[I*w*gg[-1]];
# Print["n=", np+ii+2s-2, ", Result= ", N[NR, wprec/2+5]];
# If[ii==0,PR=NR];];
# (* End of subroutine aLevinTQ /A.I. Hascelik, July 2013/ *)
#
# def main_levin():
# a=1; b=2;
# omega=100;
# prm=1/2;
# f[t_]=Exp[4t]
# g[t_]=t+Exp[4t]*Gamma[t]
#
# dg[t_]:=Derivative[1][g][t];
#
# prec=16
# wprec=2*prec
# delta = min(abs(omega*dg(a)), abs(omega*dg(b)))
# alpha = min(abs(omega*g(a)), abs(omega*g(b)))
# s=1; #(*if s>1, the integral is computed by Q_s^L*)
# test= 1 if delta<10 or alpha <=10 or s>1 else 0
#
# m = 1 if s>1 else np.floor(prec/max(np.log10(beta+1),1)+2)
# nc = 2*m+1 #(*or np=2m, number of collocation points*)
# basis=1; # (*take basis=0 for the Chebysev basis*)
# for ii in range(0, 2, 2):
# nu = np.ones((nc+ii,)) # ConstantArray[1,nc+ii];
# nu[0] = s
# nu[-1] = s
# #nu[[1]]=s;
# #nu[[-1]]=s;
# aLevinTQ[omega,a,b,f,g,nu,wprec,prm,test,basis,nc],
# #{ii,0,2,2}];
# Print["Error= ",Abs[NR-PR]];
if __name__ == '__main__': if __name__ == '__main__':
# levin_integrate()
test_docstrings() test_docstrings()
# qdemo(np.exp, 0, 3, plot_error=True) # qdemo(np.exp, 0, 3, plot_error=True)
# plt.show('hold') # plt.show('hold')

2
wafo/spectrum/__init__.py
Unescape Escape View File

@ -7,4 +7,4 @@ Spectrum package in WAFO Toolbox.
from __future__ import absolute_import from __future__ import absolute_import
from .core import * from .core import *
from . import models from . import models
from wafo.wave_theory import dispersion_relation from ..wave_theory import dispersion_relation

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

@ -1,6 +1,4 @@
from __future__ import absolute_import, division from __future__ import absolute_import, division
from wafo.misc import meshgrid, gravity, cart2polar, polar2cart
from wafo.objects import TimeSeries, mat2timeseries
import warnings import warnings
import os import os
import numpy as np import numpy as np
@ -15,15 +13,18 @@ from scipy.integrate import simps, trapz
from scipy.special import erf from scipy.special import erf
from scipy.linalg import toeplitz from scipy.linalg import toeplitz
import scipy.interpolate as interpolate import scipy.interpolate as interpolate
from scipy.interpolate.interpolate import interp1d, interp2d
from ..misc import meshgrid, gravity, cart2polar, polar2cart
from ..objects import TimeSeries, mat2timeseries
from ..interpolate import stineman_interp from ..interpolate import stineman_interp
from ..wave_theory.dispersion_relation import w2k # , k2w from ..wave_theory.dispersion_relation import w2k # , k2w
from ..containers import PlotData, now from ..containers import PlotData, now
# , tranproc
from ..misc import sub_dict_select, nextpow2, discretize, JITImport from ..misc import sub_dict_select, nextpow2, discretize, JITImport
# from wafo.graphutil import cltext
from ..kdetools import qlevels from ..kdetools import qlevels
from scipy.interpolate.interpolate import interp1d
# from wafo.transform import TrData
from ..transform.models import TrLinear
from ..plotbackend import plotbackend
try: try:
from ..gaussian import Rind from ..gaussian import Rind
@ -40,20 +41,16 @@ except ImportError:
warnings.warn('Compile the cov2mod.pyd again!') warnings.warn('Compile the cov2mod.pyd again!')
cov2mod = None cov2mod = None
# from wafo.transform import TrData
from ..transform.models import TrLinear
from ..plotbackend import plotbackend
# Trick to avoid error due to circular import # Trick to avoid error due to circular import
_WAFOCOV = JITImport('wafo.covariance') _WAFOCOV = JITImport('wafo.covariance')
__all__ = ['SpecData1D', 'SpecData2D', 'plotspec'] __all__ = ['SpecData1D', 'SpecData2D', 'plotspec']
_EPS = np.finfo(float).eps
def _set_seed(iseed): def _set_seed(iseed):
'''Set seed of random generator''' '''Set seed of random generator'''
if iseed is not None: if iseed is not None:
@ -157,64 +154,63 @@ def qtf(w, h=inf, g=9.81):
def plotspec(specdata, linetype='b-', flag=1): def plotspec(specdata, linetype='b-', flag=1):
'''
PLOTSPEC Plot a spectral density
Parameters
----------
S : SpecData1D or SpecData2D object
defining spectral density.
linetype : string
defining color and linetype, see plot for possibilities
flag : scalar integer
defining the type of plot
1D:
1 plots the density, S, (default)
2 plot 10log10(S)
3 plots both the above plots
2D:
Directional spectra: S(w,theta), S(f,theta)
1 polar plot S (default)
2 plots spectral density and the directional
spreading, int S(w,theta) dw or int S(f,theta) df
3 plots spectral density and the directional
spreading, int S(w,theta)/S(w) dw or int S(f,theta)/S(f) df
4 mesh of S
5 mesh of S in polar coordinates
6 contour plot of S
7 filled contour plot of S
Wavenumber spectra: S(k1,k2)
1 contour plot of S (default)
2 filled contour plot of S
Example
-------
>>> import numpy as np
>>> import wafo.spectrum as ws
>>> Sj = ws.models.Jonswap(Hm0=3, Tp=7)
>>> S = Sj.tospecdata()
>>> ws.plotspec(S,flag=1)
S = demospec('dir'); S2 = mkdspec(jonswap,spreading);
plotspec(S,2), hold on
# Same as previous fig. due to frequency independent spreading
plotspec(S,3,'g')
# Not the same as previous figs. due to frequency dependent spreading
plotspec(S2,2,'r')
plotspec(S2,3,'m')
# transform from angular frequency and radians to frequency and degrees
Sf = ttspec(S,'f','d'); clf
plotspec(Sf,2),
See also
dat2spec, createspec, simpson
'''
pass pass
# ''' # # label the contour levels
# PLOTSPEC Plot a spectral density
#
# Parameters
# ----------
# S : SpecData1D or SpecData2D object
# defining spectral density.
# linetype : string
# defining color and linetype, see plot for possibilities
# flag : scalar integer
# defining the type of plot
# 1D:
# 1 plots the density, S, (default)
# 2 plot 10log10(S)
# 3 plots both the above plots
# 2D:
# Directional spectra: S(w,theta), S(f,theta)
# 1 polar plot S (default)
# 2 plots spectral density and the directional
# spreading, int S(w,theta) dw or int S(f,theta) df
# 3 plots spectral density and the directional
# spreading, int S(w,theta)/S(w) dw or int S(f,theta)/S(f) df
# 4 mesh of S
# 5 mesh of S in polar coordinates
# 6 contour plot of S
# 7 filled contour plot of S
# Wavenumber spectra: S(k1,k2)
# 1 contour plot of S (default)
# 2 filled contour plot of S
#
# Example
# -------
# >>> import numpy as np
# >>> import wafo.spectrum as ws
# >>> Sj = ws.models.Jonswap(Hm0=3, Tp=7)
# >>> S = Sj.tospecdata()
# >>> ws.plotspec(S,flag=1)
#
# S = demospec('dir'); S2 = mkdspec(jonswap,spreading);
# plotspec(S,2), hold on
# # Same as previous fig. due to frequency independent spreading
# plotspec(S,3,'g')
# # Not the same as previous figs. due to frequency dependent spreading
# plotspec(S2,2,'r')
# plotspec(S2,3,'m')
# % transform from angular frequency and radians to frequency and degrees
# Sf = ttspec(S,'f','d'); clf
# plotspec(Sf,2),
#
# See also dat2spec, createspec, simpson
# '''
#
# label the contour levels
# txtFlag = 0 # txtFlag = 0
# LegendOn = 1 # LegendOn = 1
# #
#
# ftype = specdata.freqtype # options are 'f' and 'w' and 'k' # ftype = specdata.freqtype # options are 'f' and 'w' and 'k'
# data = specdata.data # data = specdata.data
# if data.ndim == 2: # if data.ndim == 2:
@ -270,10 +266,10 @@ def plotspec(specdata, linetype='b-', flag=1):
# Fn = freq[-1] # Nyquist frequency # Fn = freq[-1] # Nyquist frequency
# indm = findpeaks(data, n=4) # indm = findpeaks(data, n=4)
# maxS = data.max() # maxS = data.max()
# if isfield(S,'CI') && ~isempty(S.CI), # if isfield(S,'CI') && ~isempty(S.CI):
# maxS = maxS*S.CI(2) # maxS = maxS*S.CI(2)
# txtCI = [num2str(100*S.p), '% CI'] # txtCI = [num2str(100*S.p), '% CI']
# end # #end
# #
# Fp = freq[indm]# %peak frequency/wave number # Fp = freq[indm]# %peak frequency/wave number
# #
@ -293,7 +289,7 @@ def plotspec(specdata, linetype='b-', flag=1):
# ':', label=txt) # ':', label=txt)
# plotbackend.plot(freq, data, linetype) # plotbackend.plot(freq, data, linetype)
# specdata.labels.labelfig() # specdata.labels.labelfig()
# if isfield(S,'CI'), # if isfield(S,'CI'):
# plot(freq,S.S*S.CI(1), 'r:' ) # plot(freq,S.S*S.CI(1), 'r:' )
# plot(freq,S.S*S.CI(2), 'r:' ) # plot(freq,S.S*S.CI(2), 'r:' )
# #
@ -320,10 +316,10 @@ def plotspec(specdata, linetype='b-', flag=1):
# maxS)).repeat(len(Fp)), # maxS)).repeat(len(Fp)),
# 10 * log10(data.take(indm) / maxS))), ':',label=txt) # 10 * log10(data.take(indm) / maxS))), ':',label=txt)
# hold on # hold on
# if isfield(S,'CI'), # if isfield(S,'CI'):
# plot(freq(ind),10*log10(S.S(ind)*S.CI(1)/maxS), 'r:' ) # plot(freq(ind),10*log10(S.S(ind)*S.CI(1)/maxS), 'r:' )
# plot(freq(ind),10*log10(S.S(ind)*S.CI(2)/maxS), 'r:' ) # plot(freq(ind),10*log10(S.S(ind)*S.CI(2)/maxS), 'r:' )
# end #
# plotbackend.plot(freq[ind], 10 * log10(data[ind] / maxS), linetype) # plotbackend.plot(freq[ind], 10 * log10(data[ind] / maxS), linetype)
# #
# a = plotbackend.axis() # a = plotbackend.axis()
@ -371,7 +367,7 @@ def plotspec(specdata, linetype='b-', flag=1):
# xlabel(xlbl_txt) # xlabel(xlbl_txt)
# ylabel(xlbl_txt) # ylabel(xlbl_txt)
# title(ylbl4_txt) # title(ylbl4_txt)
# %return # # return
# km=max([-freq(1) freq(end) S.k2(1) -S.k2(end)]) # km=max([-freq(1) freq(end) S.k2(1) -S.k2(end)])
# axis([-km km -km km]) # axis([-km km -km km])
# hold on # hold on
@ -380,8 +376,8 @@ def plotspec(specdata, linetype='b-', flag=1):
# axis('square') # axis('square')
# #
# #
# %cltext(z_level) # # cltext(z_level)
# %axis('square') # # axis('square')
# if ~ih, hold off,end # if ~ih, hold off,end
# case {'dir'} # case {'dir'}
# thmin = S.theta(1)-phi;thmax=S.theta(end)-phi # thmin = S.theta(1)-phi;thmax=S.theta(end)-phi
@ -514,120 +510,13 @@ def plotspec(specdata, linetype='b-', flag=1):
# #
# if ~ih, hold off, end # if ~ih, hold off, end
# #
# % The following two commands install point-and-click editing of # # The following two commands install point-and-click editing of
# % all the text objects (title, xlabel, ylabel) of the current figure: # # all the text objects (title, xlabel, ylabel) of the current figure:
# #
# %set(findall(gcf,'type','text'),'buttondownfcn','edtext') # #set(findall(gcf,'type','text'),'buttondownfcn','edtext')
# %set(gcf,'windowbuttondownfcn','edtext(''hide'')') # #set(gcf,'windowbuttondownfcn','edtext(''hide'')')
# #
# return # return
#
#
# function fixthetalabels(thmin,thmax,xy,dim)
# %FIXTHETALABELS pretty prints the ticklabels and x or y labels for theta
# %
# % CALL fixthetalabels(thmin,thmax,xy,dim)
# %
# % thmin, thmax = minimum and maximum value for theta (wave direction)
# % xy = 'x' if theta is plotted on the x-axis
# % 'y' if theta is plotted on the y-axis
# % dim = specifies the dimension of the plot (ie number of axes shown 2 or 3)
# % If abs(thmax-thmin)<3*pi it is assumed that theta is given in radians
# % otherwise degrees
#
# ind = [('x' == xy) ('y' == xy) ];
# yx = 'yx';
# yx = yx(ind);
# if nargin<4||isempty(dim),
# dim=2;
# end
# %drawnow
# %pause
#
# if abs(thmax-thmin)<3*pi, %Radians given. Want xticks given as fractions of pi
# %Trick to update the axis
# if xy=='x'
# if dim<3,
# axis([thmin,thmax 0 inf ])
# else
# axis([thmin,thmax 0 inf 0 inf])
# end
# else
# if dim<3,
# axis([0 inf thmin,thmax ])
# else
# axis([0 inf thmin,thmax 0 inf])
# end
# end
#
# set(gca,[xy 'tick'],pi*(thmin/pi:0.25:thmax/pi));
# set(gca,[xy 'ticklabel'],[]);
# x = get(gca,[xy 'tick']);
# y = get(gca,[yx 'tick']);
# y1 = y(1);
# dy = y(2)-y1;
# yN = y(end)+dy;
# ylim = [y1 yN];
# dy1 = diff(ylim)/40;
# %ylim=get(gca,[yx 'lim'])%,ylim=ylim(2);
#
# if xy=='x'
# for j=1:length(x)
# xtxt = num2pistr(x(j));
# figtext(x(j),y1-dy1,xtxt,'data','data','center','top');
# end
# % ax = [thmin thmax 0 inf];
# ax = [thmin thmax ylim];
# if dim<3,
# figtext(mean(x),y1-7*dy1,'Wave directions (rad)','data','data','center','top')
# else
# ax = [ax 0 inf];
# xlabel('Wave directions (rad)')
# end
# else
# %ax = [0 inf thmin thmax];
# ax = [ylim thmin thmax];
#
# if dim<3,
# for j=1:length(x)
# xtxt = num2pistr(x(j));
# figtext(y1-dy1/2,x(j),xtxt,'data','data','right');
# end
# set(gca,'DefaultTextRotation',90)
# %ylabel('Wave directions (rad)')
# figtext(y1-3*dy1,mean(x),'Wave directions (rad)','data','data','center','bottom')
# set(gca,'DefaultTextRotation',0)
# else
# for j=1:length(x)
# xtxt = num2pistr(x(j));
# figtext(y1-3*dy1,x(j),xtxt,'data','data','right');
# end
# ax = [ax 0 inf];
# ylabel('Wave directions (rad)')
# end
# end
# %xtxt = num2pistr(x(j));
# %for j=2:length(x)
# % xtxt = strvcat(xtxt,num2pistr(x(j)));
# %end
# %set(gca,[xy 'ticklabel'],xtxt)
# else % Degrees given
# set(gca,[xy 'tick'],thmin:45:thmax)
# if xy=='x'
# ax=[thmin thmax 0 inf];
# if dim>=3, ax=[ax 0 inf]; end
# xlabel('Wave directions (deg)')
# else
# ax=[0 inf thmin thmax ];
# if dim>=3, ax=[ax 0 inf]; end
# ylabel('Wave directions (deg)')
# end
# end
# axis(ax)
# return
#
class SpecData1D(PlotData): class SpecData1D(PlotData):
@ -684,6 +573,43 @@ class SpecData1D(PlotData):
self.setlabels() self.setlabels()
def _get_default_dt_and_rate(self, dt):
dt_old = self.sampling_period()
if dt is None:
return dt_old, 1
rate = max(round(dt_old * 1. / dt), 1.)
return dt, rate
def _check_dt(self, dt):
freq = self.args
checkdt = 1.2 * min(diff(freq)) / 2. / pi
if self.freqtype in 'f':
checkdt *= 2 * pi
if (checkdt < 2. ** -16 / dt):
print('Step dt = %g in computation of the density is ' +
'too small.' % dt)
print('The computed covariance (by FFT(2^K)) may differ from the')
print('theoretical. Solution:')
raise ValueError('use larger dt or sparser grid for spectrum.')
def _check_cov_matrix(self, acfmat, nt, dt):
eps0 = 0.0001
if nt + 1 >= 5:
cc2 = acfmat[0, 0] - acfmat[4, 0] * (acfmat[4, 0] / acfmat[0, 0])
if (cc2 < eps0):
warnings.warn('Step dt = %g in computation of the density ' +
'is too small.' % dt)
cc1 = acfmat[0, 0] - acfmat[1, 0] * (acfmat[1, 0] / acfmat[0, 0])
if (cc1 < eps0):
warnings.warn('Step dt = %g is small, and may cause numerical ' +
'inaccuracies.' % dt)
@property
def lagtype(self):
if self.freqtype in 'k': # options are 'f' and 'w' and 'k'
return 'x'
return 't'
def tocov_matrix(self, nr=0, nt=None, dt=None): def tocov_matrix(self, nr=0, nt=None, dt=None):
''' '''
Computes covariance function and its derivatives, alternative version Computes covariance function and its derivatives, alternative version
@ -725,59 +651,30 @@ class SpecData1D(PlotData):
objects objects
''' '''
ftype = self.freqtype # %options are 'f' and 'w' and 'k'
dt, rate = self._get_default_dt_and_rate(dt)
self._check_dt(dt)
freq = self.args freq = self.args
n_f = len(freq) n_f = len(freq)
dt_old = self.sampling_period()
if dt is None:
dt = dt_old
rate = 1
else:
rate = max(round(dt_old * 1. / dt), 1.)
if nt is None: if nt is None:
nt = rate * (n_f - 1) nt = rate * (n_f - 1)
else: # %check if Nt is ok else: # %check if Nt is ok
nt = minimum(nt, rate * (n_f - 1)) nt = minimum(nt, rate * (n_f - 1))
checkdt = 1.2 * min(diff(freq)) / 2. / pi
if ftype in 'k':
lagtype = 'x'
else:
lagtype = 't'
if ftype in 'f':
checkdt = checkdt * 2 * pi
msg1 = 'Step dt = %g in computation of the density is too small.' % dt
msg2 = 'Step dt = %g is small, and may cause numerical inaccuracies.' % dt
if (checkdt < 2. ** -16 / dt):
print(msg1)
print('The computed covariance (by FFT(2^K)) may differ from the')
print('theoretical. Solution:')
raise ValueError('use larger dt or sparser grid for spectrum.')
# Calculating covariances
#~~~~~~~~~~~~~~~~~~~~~~~~
spec = self.copy() spec = self.copy()
spec.resample(dt) spec.resample(dt)
acf = spec.tocovdata(nr, nt, rate=1) acf = spec.tocovdata(nr, nt, rate=1)
acfmat = zeros((nt + 1, nr + 1), dtype=float) acfmat = zeros((nt + 1, nr + 1), dtype=float)
acfmat[:, 0] = acf.data[0:nt + 1] acfmat[:, 0] = acf.data[0:nt + 1]
fieldname = 'R' + lagtype * nr fieldname = 'R' + self.lagtype * nr
for i in range(1, nr + 1): for i in range(1, nr + 1):
fname = fieldname[:i + 1] fname = fieldname[:i + 1]
r_i = getattr(acf, fname) r_i = getattr(acf, fname)
acfmat[:, i] = r_i[0:nt + 1] acfmat[:, i] = r_i[0:nt + 1]
eps0 = 0.0001 self._check_cov_matrix(acfmat, nt, dt)
if nt + 1 >= 5:
cc2 = acfmat[0, 0] - acfmat[4, 0] * (acfmat[4, 0] / acfmat[0, 0])
if (cc2 < eps0):
warnings.warn(msg1)
cc1 = acfmat[0, 0] - acfmat[1, 0] * (acfmat[1, 0] / acfmat[0, 0])
if (cc1 < eps0):
warnings.warn(msg2)
return acfmat return acfmat
def tocovdata(self, nr=0, nt=None, rate=None): def tocovdata(self, nr=0, nt=None, rate=None):
@ -1657,8 +1554,8 @@ class SpecData1D(PlotData):
# ftmp = cov2mmtpdfexe(R,dt,u,defnr,Nstart,hg,options) # ftmp = cov2mmtpdfexe(R,dt,u,defnr,Nstart,hg,options)
# err = repmat(nan,size(ftmp)) # err = repmat(nan,size(ftmp))
# else # else
[ftmp, err, terr, options] = self._cov2mmtpdf( ftmp, err, terr, options = self._cov2mmtpdf(R, dt, u, defnr, Nstart,
R, dt, u, defnr, Nstart, hg, options) hg, options)
# end # end
note = '' note = ''
@ -1850,8 +1747,8 @@ 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, R5 = R[:, :5].T R0, R1, R2, R3, R4 = R[:, :5].T
covinput = self._covinput_mmt_pdf
Ntime = len(R0) Ntime = len(R0)
Nx0 = max(1, len(hg)) Nx0 = max(1, len(hg))
Nx1 = Nx0 Nx1 = Nx0
@ -1876,14 +1773,15 @@ class SpecData1D(PlotData):
if def_nr <= 1: # just plain Mm if def_nr <= 1: # just plain Mm
Nx = Nx1 * (Nx1 - 1) / 2 Nx = Nx1 * (Nx1 - 1) / 2
IJ = (Nx1 + isOdd) / 2 IJ = (Nx1 + isOdd) / 2
if (hg[0] + hg[Nx1 - 1] == 0 and (hg[IJ - 1] == 0 or hg[IJ - 1] + hg[IJ] == 0)): if (hg[0] + hg[Nx1 - 1] == 0 and (hg[IJ - 1] == 0 or
hg[IJ - 1] + hg[IJ] == 0)):
symmetry = 0 symmetry = 0
print(' Integration region symmetric') print(' Integration region symmetric')
# May save Nx1-isOdd integrations in each time step # May save Nx1-isOdd integrations in each time step
# This is not implemented yet. # This is not implemented yet.
# Nx = Nx1*(Nx1-1)/2-Nx1+isOdd # Nx = Nx1*(Nx1-1)/2-Nx1+isOdd
# normalizing constant:
# CC = normalizing constant = 1/ expected number of zero-up-crossings of X' # CC = 1/ expected number of zero-up-crossings of X'
# CC = 2*pi*sqrt(-R2[0]/R4[0]) # CC = 2*pi*sqrt(-R2[0]/R4[0])
# XcScale = log(CC) # XcScale = log(CC)
XcScale = log(2 * pi * sqrt(-R2[0] / R4[0])) XcScale = log(2 * pi * sqrt(-R2[0] / R4[0]))
@ -1902,16 +1800,15 @@ class SpecData1D(PlotData):
Nc = 5 Nc = 5
NI = 5 NI = 5
Nd = 3 Nd = 3
# CC= normalizing constant= 1/ expected number of u-up-crossings of X # CC = 1/ expected number of u-up-crossings of X
# CC = 2*pi*sqrt(-R0(1)/R2(1))*exp(0.5D0*u*u/R0(1)) # CC = 2*pi*sqrt(-R0(1)/R2(1))*exp(0.5D0*u*u/R0(1))
XcScale = log(2 * pi * sqrt(-R0[0] / R2[0])) + 0.5 * u * u / R0[0] XcScale = log(2 * pi * sqrt(-R0[0] / R2[0])) + 0.5 * u * u / R0[0]
options['xcscale'] = XcScale options['xcscale'] = XcScale
opt0 = struct2cell(options) # opt0 = [options[n] for n in ('SCIS', 'XcScale', 'ABSEPS', 'RELEPS',
# opt0 = opt0(1:10) # 'COVEPS', 'MAXPTS', 'MINPTS', 'seed',
# seed = [] # 'NIT1')]
# opt0 = {SCIS,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS,seed,NIT1} rind = Rind(**options)
if (Nx > 1): if (Nx > 1):
# (M,m) or (M,m)v distribution wanted # (M,m) or (M,m)v distribution wanted
if ((def_nr == 0 or def_nr == 2)): if ((def_nr == 0 or def_nr == 2)):
@ -1982,10 +1879,10 @@ class SpecData1D(PlotData):
indI[2] = Nt + 1 indI[2] = Nt + 1
indI[3] = Ntd indI[3] = Ntd
# positive wave period # positive wave period
BIG[:Ntdc, :Ntdc] = covinput( # self._covinput_mmt_pdf(BIG, R, tn, ts, tnold)
BIG[:Ntdc, :Ntdc], R0, R1, R2, R3, R4, Ntd, 0) BIG[:Ntdc, :Ntdc] = covinput(BIG[:Ntdc, :Ntdc], R, Ntd, 0)
[fxind, err0, terr0] = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc], [fxind, err0, terr0] = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc],
a_lo, a_up, indI, xc, Nt, *opt0) a_lo, a_up, indI, xc, Nt)
# fxind = CC*rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),xc,Nt,NIT1, # fxind = CC*rind(BIG(1:Ntdc,1:Ntdc),ex(1:Ntdc),xc,Nt,NIT1,
# speed1,indI,a_lo,a_up) # speed1,indI,a_lo,a_up)
if (Nx < 2): if (Nx < 2):
@ -2001,12 +1898,9 @@ class SpecData1D(PlotData):
if def_nr in [-2, -1, 0]: if def_nr in [-2, -1, 0]:
for i in range(1, Nx1): for i in range(1, Nx1):
J = IJ + i J = IJ + i
pdf[:i, i, 0] = pdf[:i, i, 0] + \ pdf[:i, i, 0] += fxind[IJ:J].T * dt # *CC
fxind[IJ:J].T * dt # *CC err[:i, i, 0] += (err0[IJ + 1:J].T * dt) ** 2
err[:i, i, 0] = err[:i, i, 0] + \ terr[:i, i, 0] += (terr0[IJ:J].T * dt)
(err0[IJ + 1:J].T * dt) ** 2
terr[:i, i, 0] = terr[
:i, i, 0] + (terr0[IJ:J].T * dt)
IJ = J IJ = J
elif def_nr == 1: # joint density of (M,m,TMm) elif def_nr == 1: # joint density of (M,m,TMm)
for i in range(1, Nx1): for i in range(1, Nx1):
@ -2020,24 +1914,18 @@ class SpecData1D(PlotData):
elif def_nr == 2: elif def_nr == 2:
for i in range(1, Nx1): for i in range(1, Nx1):
J = IJ + Nx1 J = IJ + Nx1
pdf[1:Nx1, i, 0] = pdf[1:Nx1, i, 0] + \ pdf[1:Nx1, i, 0] += fxind[IJ:J].T * dt # %*CC
fxind[IJ:J].T * dt # %*CC err[1:Nx1, i, 0] += (err0[IJ:J].T * dt) ** 2
err[1:Nx1, i, 0] = err[1:Nx1, i, 0] + \ terr[1:Nx1, i, 0] += (terr0[IJ:J].T * dt)
(err0[IJ:J].T * dt) ** 2
terr[1:Nx1, i, 0] = terr[
1:Nx1, i, 0] + (terr0[IJ:J].T * dt)
IJ = J IJ = J
# end %do # end %do
elif def_nr == 3: elif def_nr == 3:
# joint density of level v separated (M,m,TMm)v # joint density of level v separated (M,m,TMm)v
for i in range(1, Nx1): for i in range(1, Nx1):
J = IJ + Nx1 J = IJ + Nx1
pdf[1:Nx1, i, Ntd] = pdf[ pdf[1:Nx1, i, Ntd] += fxind[IJ:J].T # %*CC
1:Nx1, i, Ntd] + fxind[IJ:J].T # %*CC err[1:Nx1, i, Ntd] += (err0[IJ:J].T) ** 2
err[1:Nx1, i, Ntd] = err[ terr[1:Nx1, i, Ntd] += (terr0[IJ:J].T)
1:Nx1, i, Ntd] + (err0[IJ:J].T) ** 2
terr[1:Nx1, i, Ntd] = terr[
1:Nx1, i, Ntd] + (terr0[IJ:J].T)
IJ = J IJ = J
# end %do # end %do
# end % SELECT # end % SELECT
@ -2067,10 +1955,9 @@ class SpecData1D(PlotData):
for ts in range(1, tn - 1): # = 2:tn-1: for ts in range(1, tn - 1): # = 2:tn-1:
# positive wave period # positive wave period
BIG[:Ntdc, :Ntdc] = covinput(BIG[:Ntdc, :Ntdc], BIG[:Ntdc, :Ntdc] = covinput(BIG[:Ntdc, :Ntdc],
R0, R1, R2, R3, R4, R, tn, ts, tnold)
tn, ts, tnold) fxind, err0, terr0 = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc],
[fxind, err0, terr0] = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc], a_lo, a_up, indI, xc, Nt)
a_lo, a_up, indI, xc, Nt, *opt0)
# tnold = tn # tnold = tn
if def_nr in [3, 4]: if def_nr in [3, 4]:
@ -2087,12 +1974,9 @@ class SpecData1D(PlotData):
IJ = 0 IJ = 0
for i in range(1, Nx1): for i in range(1, Nx1):
J = IJ + Nx1 J = IJ + Nx1
pdf[1:Nx1, i, ts] = pdf[ pdf[1:Nx1, i, ts] += fxind[IJ:J].T * dt
1:Nx1, i, ts] + fxind[IJ:J].T * dt err[1:Nx1, i, ts] += (err0[IJ:J].T * dt) ** 2
err[1:Nx1, i, ts] = err[ terr[1:Nx1, i, ts] += (terr0[IJ:J].T * dt)
1:Nx1, i, ts] + (err0[IJ:J].T * dt) ** 2
terr[1:Nx1, i, ts] = terr[
1:Nx1, i, ts] + (terr0[IJ:J].T * dt)
IJ = J IJ = J
# end %do # end %do
# end # end
@ -2111,12 +1995,9 @@ class SpecData1D(PlotData):
for i in range(1, Nx1): # = 2:Nx1 for i in range(1, Nx1): # = 2:Nx1
J = IJ + Nx1 J = IJ + Nx1
# %*CC # %*CC
pdf[1:Nx1, i, tn - ts] = pdf[1:Nx1, pdf[1:Nx1, i, tn - ts] += fxind[IJ:J].T * dt
i, tn - ts] + fxind[IJ:J].T * dt err[1:Nx1, i, tn - ts] += (err0[IJ:J].T * dt) ** 2
err[1:Nx1, i, tn - ts] = err[1:Nx1, i, terr[1:Nx1, i, tn - ts] += (terr0[IJ:J].T * dt)
tn - ts] + (err0[IJ:J].T * dt) ** 2
terr[
1:Nx1, i, tn - ts] = terr[1:Nx1, i, tn - ts] + (terr0[IJ:J].T * dt)
IJ = J IJ = J
# end %do # end %do
# end # end
@ -2129,10 +2010,9 @@ class SpecData1D(PlotData):
# linear. # linear.
# positive wave period # positive wave period
BIG[:Ntdc, :Ntdc] = covinput(BIG[:Ntdc, :Ntdc], BIG[:Ntdc, :Ntdc] = covinput(BIG[:Ntdc, :Ntdc],
R0, R1, R2, R3, R4, R, tn, ts, tnold)
tn, ts, tnold) fxind, err0, terr0 = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc],
[fxind, err0, terr0] = rind(BIG[:Ntdc, :Ntdc], ex[:Ntdc], a_lo, a_up, indI, xc, Nt)
a_lo, a_up, indI, xc, Nt, *opt0)
#[fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex,a_lo,a_up,indI, xc,Nt,opt0{:}) #[fxind,err0] = rind(BIG(1:Ntdc,1:Ntdc),ex,a_lo,a_up,indI, xc,Nt,opt0{:})
# tnold = tn # tnold = tn
@ -2156,21 +2036,15 @@ class SpecData1D(PlotData):
# Max to the crossing of level u (M,m,TMd). # Max to the crossing of level u (M,m,TMd).
for i in range(1, Nx1): for i in range(1, Nx1):
J = IJ + Nx1 J = IJ + Nx1
pdf[1:Nx1, i, ts] = pdf[ pdf[1:Nx1, i, ts] += fxind[IJ:J] * dt # %*CC
1:Nx1, i, ts] + fxind[IJ:J] * dt # %*CC err[1:Nx1, i, ts] += (err0[IJ:J] * dt) ** 2
err[1:Nx1, i, ts] = err[ terr[1:Nx1, i, ts] += (terr0[IJ:J] * dt)
1:Nx1, i, ts] + (err0[IJ:J] * dt) ** 2
terr[1:Nx1, i, ts] = terr[
1:Nx1, i, ts] + (terr0[IJ:J] * dt)
if (ts < tn - ts): if (ts < tn - ts):
# exploiting the symmetry # exploiting the symmetry
# %*CC # %*CC
pdf[i, 1:Nx1, tn - ts] = pdf[i, pdf[i, 1:Nx1, tn - ts] += fxind[IJ:J] * dt
1:Nx1, tn - ts] + fxind[IJ:J] * dt err[i, 1:Nx1, tn - ts] += (err0[IJ:J] * dt) ** 2
err[i, 1:Nx1, tn - ts] = err[i, 1:Nx1, terr[i, 1:Nx1, tn - ts] += (terr0[IJ:J] * dt)
tn - ts] + (err0[IJ:J] * dt) ** 2
terr[
i, 1:Nx1, tn - ts] = terr[i, 1:Nx1, tn - ts] + (terr0[IJ:J] * dt)
# end # end
IJ = J IJ = J
# end %do # end %do
@ -2179,20 +2053,15 @@ class SpecData1D(PlotData):
# from the crossing of level u to min (M,m,Tdm). # from the crossing of level u to min (M,m,Tdm).
for i in range(1, Nx1): # = 2:Nx1, for i in range(1, Nx1): # = 2:Nx1,
J = IJ + Nx1 J = IJ + Nx1
pdf[1:Nx1, i, tn - ts] = pdf[1:Nx1, pdf[1:Nx1, i, tn - ts] += fxind[IJ:J] * dt
i, tn - ts] + fxind[IJ:J] * dt err[1:Nx1, i, tn - ts] += (err0[IJ:J] * dt) ** 2
err[1:Nx1, i, tn - ts] = err[1:Nx1, i,
tn - ts] + (err0[IJ:J] * dt) ** 2
terr[ terr[
1:Nx1, i, tn - ts] = terr[1:Nx1, i, tn - ts] + (terr0[IJ:J] * dt) 1:Nx1, i, tn - ts] += (terr0[IJ:J] * dt)
if (ts < tn - ts + 1): if (ts < tn - ts + 1):
# exploiting the symmetry # exploiting the symmetry
pdf[i, 1:Nx1, ts] = pdf[ pdf[i, 1:Nx1, ts] += fxind[IJ:J] * dt
i, 1:Nx1, ts] + fxind[IJ:J] * dt err[i, 1:Nx1, ts] += (err0[IJ:J] * dt) ** 2
err[i, 1:Nx1, ts] = err[ terr[i, 1:Nx1, ts] += (terr0[IJ:J] * dt)
i, 1:Nx1, ts] + (err0[IJ:J] * dt) ** 2
terr[i, 1:Nx1, ts] = terr[
i, 1:Nx1, ts] + (terr0[IJ:J] * dt)
# end %ENDIF # end %ENDIF
IJ = J IJ = J
# end %do # end %do
@ -2263,7 +2132,7 @@ class SpecData1D(PlotData):
Cov(X''(t),X(s)) = r''(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|) Cov(X''(t),X''(s)) = r''''(s-t) = r''''(|s-t|)
""" """
R0, R1, R2, R3, R4 = R.T R0, R1, R2, R3, R4 = R[:, :5].T
if (ts > 1): if (ts > 1):
shft = 1 shft = 1
N = tn + 5 + shft N = tn + 5 + shft
@ -4006,17 +3875,17 @@ class SpecData2D(PlotData):
# ftype = self.freqtype # ftype = self.freqtype
freq = self.args[1] freq = self.args[1]
theta = linspace(-pi, pi, ntOld) theta = linspace(-pi, pi, ntOld)
[F, T] = meshgrid(freq, theta) # [F, T] = meshgrid(freq, theta)
dtheta = self.theta[1] - self.theta[0] dtheta = self.theta[1] - self.theta[0]
self.theta[nt] = self.theta[nt - 1] + dtheta self.theta[nt] = self.theta[nt - 1] + dtheta
self.data[nt, :] = self.data[0, :] self.data[nt, :] = self.data[0, :]
self.data = interp2(freq, self.data = interp2d(freq,
np.vstack([self.theta[0] - dtheta, np.vstack([self.theta[0] - dtheta,
self.theta]), self.theta]),
np.vstack([self.data[nt, :], np.vstack([self.data[nt, :],
self.data]), F, T, self.data]),
method) kind=method)(freq, theta)
self.args[0] = theta self.args[0] = theta
elif stype == 'k2d': elif stype == 'k2d':
@ -4274,33 +4143,16 @@ class SpecData2D(PlotData):
self.labels.zlab = labels[2] self.labels.zlab = labels[2]
def _test_specdata():
import wafo.spectrum.models as sm
Sj = sm.Jonswap()
S = Sj.tospecdata()
me, va, sk, ku = S.stats_nl(moments='mvsk')
def main(): def main():
import matplotlib import matplotlib
matplotlib.interactive(True) matplotlib.interactive(True)
from wafo.spectrum import models as sm from wafo.spectrum import models as sm
w = linspace(0, 3, 100)
Sj = sm.Jonswap() Sj = sm.Jonswap()
S = Sj.tospecdata() S = Sj.tospecdata()
f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7)
f.err
f.plot()
f.show()
# pdfplot(f)
# hold on,
# plot(f.x{:}, f.f+f.err,'r',f.x{:}, f.f-f.err) estimated error bounds
# hold off
# S = SpecData1D(Sj(w),w)
R = S.tocovdata(nr=1) R = S.tocovdata(nr=1)
S1 = S.copy()
Si = R.tospecdata() Si = R.tospecdata()
ns = 5000 ns = 5000
dt = .2 dt = .2

81
wafo/spectrum/models.py
Unescape Escape View File

@ -1,3 +1,4 @@
#!/usr/bin/env python
""" """
Models module Models module
------------- -------------
@ -26,17 +27,7 @@ Spreading - Directional spreading function.
""" """
# Name: models from __future__ import absolute_import, division
# Purpose: Interface to various spectrum models
#
# Author: pab
#
# Created: 29.08.2008
# Copyright: (c) pab 2008
# Licence: <your licence>
#!/usr/bin/env python
from __future__ import division
import warnings import warnings
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
@ -50,17 +41,20 @@ from numpy import (inf, atleast_1d, newaxis, any, minimum, maximum, array,
cos, abs, sinh, isfinite, mod, expm1, tanh, cosh, finfo, cos, abs, sinh, isfinite, mod, expm1, tanh, cosh, finfo,
ones, ones_like, isnan, zeros_like, flatnonzero, sinc, ones, ones_like, isnan, zeros_like, flatnonzero, sinc,
hstack, vstack, real, flipud, clip) hstack, vstack, real, flipud, clip)
from wafo.wave_theory.dispersion_relation import w2k, k2w # @UnusedImport from ..wave_theory.dispersion_relation import w2k, k2w # @UnusedImport
from wafo.spectrum import SpecData1D, SpecData2D from .core import SpecData1D, SpecData2D
sech = lambda x: 1.0 / cosh(x)
eps = finfo(float).eps
__all__ = ['Bretschneider', 'Jonswap', 'Torsethaugen', 'Wallop', 'McCormick', __all__ = ['Bretschneider', 'Jonswap', 'Torsethaugen', 'Wallop', 'McCormick',
'OchiHubble', 'Tmaspec', 'jonswap_peakfact', 'jonswap_seastate', 'OchiHubble', 'Tmaspec', 'jonswap_peakfact', 'jonswap_seastate',
'spreading', 'w2k', 'k2w', 'phi1'] 'spreading', 'w2k', 'k2w', 'phi1']
_EPS = finfo(float).eps
def sech(x):
return 1.0 / cosh(x)
def _gengamspec(wn, N=5, M=4): def _gengamspec(wn, N=5, M=4):
''' Return Generalized gamma spectrum in dimensionless form ''' Return Generalized gamma spectrum in dimensionless form
@ -543,7 +537,7 @@ class Jonswap(ModelSpectrum):
self.method = method self.method = method
self.wnc = wnc self.wnc = wnc
if self.gamma == None or not isfinite(self.gamma) or self.gamma < 1: if self.gamma is None or not isfinite(self.gamma) or self.gamma < 1:
self.gamma = jonswap_peakfact(Hm0, Tp) self.gamma = jonswap_peakfact(Hm0, Tp)
self._preCalculateAg() self._preCalculateAg()
@ -580,7 +574,7 @@ class Jonswap(ModelSpectrum):
if self.gamma == 1: if self.gamma == 1:
self.Ag = 1.0 self.Ag = 1.0
self.method = 'parametric' self.method = 'parametric'
elif self.Ag != None: elif self.Ag is not None:
self.method = 'custom' self.method = 'custom'
if self.Ag <= 0: if self.Ag <= 0:
raise ValueError('Ag must be larger than 0!') raise ValueError('Ag must be larger than 0!')
@ -615,7 +609,7 @@ class Jonswap(ModelSpectrum):
# elseif N == 5 && M == 4, # elseif N == 5 && M == 4,
# options.Ag = (1+1.0*log(gammai).**1.16)/gammai # options.Ag = (1+1.0*log(gammai).**1.16)/gammai
# options.Ag = (1-0.287*log(gammai)) # options.Ag = (1-0.287*log(gammai))
### options.normalizeMethod = 'Three' # options.normalizeMethod = 'Three'
# elseif N == 4 && M == 4, # elseif N == 4 && M == 4,
# options.Ag = (1+1.1*log(gammai).**1.19)/gammai # options.Ag = (1+1.1*log(gammai).**1.19)/gammai
else: else:
@ -624,7 +618,7 @@ class Jonswap(ModelSpectrum):
if self.sigmaA != 0.07 or self.sigmaB != 0.09: if self.sigmaA != 0.07 or self.sigmaB != 0.09:
warnings.warn('Use integration to calculate Ag when ' + warnings.warn('Use integration to calculate Ag when ' +
'sigmaA~=0.07 or sigmaB~=0.09') 'sigmaA!=0.07 or sigmaB!=0.09')
def peak_e_factor(self, wn): def peak_e_factor(self, wn):
''' PEAKENHANCEMENTFACTOR ''' PEAKENHANCEMENTFACTOR
@ -790,9 +784,9 @@ class Tmaspec(Jonswap):
self.type = 'TMA' self.type = 'TMA'
def phi(self, w, h=None, g=None): def phi(self, w, h=None, g=None):
if h == None: if h is None:
h = self.h h = self.h
if g == None: if g is None:
g = self.g g = self.g
return phi1(w, h, g) return phi1(w, h, g)
@ -1095,14 +1089,14 @@ class McCormick(Bretschneider):
self.type = 'McCormick' self.type = 'McCormick'
self.Hm0 = Hm0 self.Hm0 = Hm0
self.Tp = Tp self.Tp = Tp
if Tz == None: if Tz is None:
Tz = 0.8143 * Tp Tz = 0.8143 * Tp
self.Tz = Tz self.Tz = Tz
if chk_seastate: if chk_seastate:
self.chk_seastate() self.chk_seastate()
if M == None and self.Hm0 > 0: if M is None and self.Hm0 > 0:
self._TpdTz = Tp / Tz self._TpdTz = Tp / Tz
M = 1.0 / optimize.fminbound(self._localoptfun, 0.01, 5) M = 1.0 / optimize.fminbound(self._localoptfun, 0.01, 5)
self.M = M self.M = M
@ -1410,7 +1404,7 @@ class Spreading(object):
5 to 30 being a function of dimensionless wind speed. 5 to 30 being a function of dimensionless wind speed.
However, Goda and Suzuki (1975) proposed SP = 10 for wind waves, SP = 25 However, Goda and Suzuki (1975) proposed SP = 10 for wind waves, SP = 25
for swell with short decay distance and SP = 75 for long decay distance. for swell with short decay distance and SP = 75 for long decay distance.
Compared to experiments Krogstad et al. (1998) found that m_a = 5 +/- eps Compared to experiments Krogstad et al. (1998) found that m_a = 5 +/- _EPS
and that -1< m_b < -3.5. and that -1< m_b < -3.5.
Values given in the litterature: [s_a s_b m_a m_b wn_lo wn_c wn_up] Values given in the litterature: [s_a s_b m_a m_b wn_lo wn_c wn_up]
(Mitsuyasu: s_a == s_b) (cos-2s) [15 15 5 -2.5 0 1 3 ] (Mitsuyasu: s_a == s_b) (cos-2s) [15 15 5 -2.5 0 1 3 ]
@ -1510,7 +1504,7 @@ class Spreading(object):
if not self.method[0] in methods: if not self.method[0] in methods:
raise ValueError('Unknown method') raise ValueError('Unknown method')
self.method = methods[self.method[0]] self.method = methods[self.method[0]]
elif self.method == None: elif self.method is None:
pass pass
else: else:
if method < 0 or 3 < method: if method < 0 or 3 < method:
@ -1561,7 +1555,7 @@ class Spreading(object):
'The length of theta0 must equal to 1 or the length of w') 'The length of theta0 must equal to 1 or the length of w')
else: else:
TH = mod(theta - th0 + pi, 2 * pi) - pi # make sure -pi<=TH<pi TH = mod(theta - th0 + pi, 2 * pi) - pi # make sure -pi<=TH<pi
if self.method != None: # frequency dependent spreading if self.method is not None: # frequency dependent spreading
TH = TH[:, newaxis] TH = TH[:, newaxis]
# Get spreading parameter # Get spreading parameter
@ -1574,7 +1568,7 @@ class Spreading(object):
r1 = abs(s / (s + 1)) r1 = abs(s / (s + 1))
# Find distribution parameter from first Fourier coefficient. # Find distribution parameter from first Fourier coefficient.
s_par = self.fourier2distpar(r1) s_par = self.fourier2distpar(r1)
if self.method != None: if self.method is not None:
s_par = s_par[newaxis, :] s_par = s_par[newaxis, :]
return s_par, TH, phi0, Nt return s_par, TH, phi0, Nt
@ -1823,7 +1817,7 @@ class Spreading(object):
A[ix] = Ai + 0.5 * (da[ix] - Ai) * (Ai <= 0.0) A[ix] = Ai + 0.5 * (da[ix] - Ai) * (Ai <= 0.0)
ix = flatnonzero( ix = flatnonzero(
(abs(da) > sqrt(eps) * abs(A)) * (abs(da) > sqrt(eps))) (abs(da) > sqrt(_EPS) * abs(A)) * (abs(da) > sqrt(_EPS)))
if ix.size == 0: if ix.size == 0:
if any(A > pi): if any(A > pi):
raise ValueError( raise ValueError(
@ -1837,8 +1831,10 @@ class Spreading(object):
''' '''
Returns the solution of R1 = besseli(1,K)/besseli(0,K), Returns the solution of R1 = besseli(1,K)/besseli(0,K),
''' '''
def fun0(x):
return sp.ive(1, x) / sp.ive(0, x)
K0 = hstack((linspace(0, 10, 513), linspace(10.00001, 100))) K0 = hstack((linspace(0, 10, 513), linspace(10.00001, 100)))
fun0 = lambda x: sp.ive(1, x) / sp.ive(0, x)
funK = interp1d(fun0(K0), K0) funK = interp1d(fun0(K0), K0)
K0 = funK(r1.ravel()) K0 = funK(r1.ravel())
k1 = flatnonzero(isnan(K0)) k1 = flatnonzero(isnan(K0))
@ -1848,15 +1844,14 @@ class Spreading(object):
ix0 = flatnonzero(r1 != 0.0) ix0 = flatnonzero(r1 != 0.0)
K = zeros_like(r1) K = zeros_like(r1)
fun = lambda x: fun0(x) - r1[ix]
for ix in ix0: for ix in ix0:
K[ix] = optimize.fsolve(fun, K0[ix]) K[ix] = optimize.fsolve(lambda x: fun0(x) - r1[ix], K0[ix])
return K return K
def fourier2b(self, r1): def fourier2b(self, r1):
''' Returns the solution of R1 = pi/(2*B*sinh(pi/(2*B)). ''' Returns the solution of R1 = pi/(2*B*sinh(pi/(2*B)).
''' '''
B0 = hstack((linspace(eps, 5, 513), linspace(5.0001, 100))) B0 = hstack((linspace(_EPS, 5, 513), linspace(5.0001, 100)))
funB = interp1d(self._r1ofsech2(B0), B0) funB = interp1d(self._r1ofsech2(B0), B0)
B0 = funB(r1.ravel()) B0 = funB(r1.ravel())
@ -1867,7 +1862,9 @@ class Spreading(object):
ix0 = flatnonzero(r1 != 0.0) ix0 = flatnonzero(r1 != 0.0)
B = zeros_like(r1) B = zeros_like(r1)
fun = lambda x: 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix]
def fun(x):
return 0.5 * pi / (sinh(.5 * pi / x)) - x * r1[ix]
for ix in ix0: for ix in ix0:
B[ix] = abs(optimize.fsolve(fun, B0[ix])) B[ix] = abs(optimize.fsolve(fun, B0[ix]))
return B return B
@ -1890,7 +1887,7 @@ class Spreading(object):
S : ndarray S : ndarray
spread parameter of COS2S functions spread parameter of COS2S functions
''' '''
if self.method == None: if self.method is None:
# no frequency dependent spreading, # no frequency dependent spreading,
# but possible frequency dependent direction # but possible frequency dependent direction
s = atleast_1d(self.s_a) s = atleast_1d(self.s_a)
@ -1923,12 +1920,12 @@ class Spreading(object):
# Banner parametrization for B in SECH-2 # Banner parametrization for B in SECH-2
s3m = spb * (wn_up) ** mb s3m = spb * (wn_up) ** mb
s3p = self._donelan(wn_up) s3p = self._donelan(wn_up)
# % Scale so that parametrization will be continous # Scale so that parametrization will be continous
scale = s3m / s3p scale = s3m / s3p
s[k] = scale * self.donelan(wn[k]) s[k] = scale * self.donelan(wn[k])
r1 = self.r1ofsech2(s) r1 = self.r1ofsech2(s)
#% Convert to S-paramater in COS-2S distribution # Convert to S-paramater in COS-2S distribution
s = r1 / (1. - r1) s = r1 / (1. - r1)
else: else:
s[k] = 0.0 s[k] = 0.0
@ -1994,7 +1991,7 @@ class Spreading(object):
theta = np.linspace(-pi, pi, nt) theta = np.linspace(-pi, pi, nt)
else: else:
L = abs(theta[-1] - theta[0]) L = abs(theta[-1] - theta[0])
if abs(L - pi) > eps: if abs(L - pi) > _EPS:
raise ValueError('theta must cover all angles -pi -> pi') raise ValueError('theta must cover all angles -pi -> pi')
nt = len(theta) nt = len(theta)
@ -2036,11 +2033,9 @@ def _test_some_spectra():
plb.show() plb.show()
plb.close('all') plb.close('all')
#import pylab as plb
#w = plb.linspace(0,3)
w, th = plb.ogrid[0:4, 0:6] w, th = plb.ogrid[0:4, 0:6]
k, k2 = w2k(w, th) k, k2 = w2k(w, th)
#k1, k12 = w2k(w, th, h=20)
plb.plot(w, k, w, k2) plb.plot(w, k, w, k2)
plb.show() plb.show()
@ -2080,8 +2075,8 @@ def _test_spreading():
pi = plb.pi pi = plb.pi
w = plb.linspace(0, 3, 257) w = plb.linspace(0, 3, 257)
theta = plb.linspace(-pi, pi, 129) theta = plb.linspace(-pi, pi, 129)
theta0 = lambda w: w * plb.pi / 6.0
D2 = Spreading('cos2s', theta0=theta0) D2 = Spreading('cos2s', theta0=lambda w: w * plb.pi / 6.0)
d1 = D2(theta, w)[0] d1 = D2(theta, w)[0]
_t = plb.contour(d1.squeeze()) _t = plb.contour(d1.squeeze())

1
wafo/spectrum/tests/__init__.py
Unescape Escape View File

@ -1 +0,0 @@

176
wafo/spectrum/tests/test_specdata1d.py
Unescape Escape View File

@ -3,7 +3,8 @@ import wafo.transform.models as wtm
import wafo.objects as wo import wafo.objects as wo
from wafo.spectrum import SpecData1D from wafo.spectrum import SpecData1D
import numpy as np import numpy as np
from numpy.testing import assert_array_almost_equal from numpy import NAN
from numpy.testing import assert_array_almost_equal, assert_array_equal
import unittest import unittest
@ -12,36 +13,32 @@ def slow(f):
return f return f
class TestSpectrum(unittest.TestCase): class TestSpectrumHs7(unittest.TestCase):
def setUp(self):
self.Sj = sm.Jonswap(Hm0=7.0, Tp=11)
self.S = self.Sj.tospecdata()
@slow
def test_tocovmatrix(self): def test_tocovmatrix(self):
Sj = sm.Jonswap() acfmat = self.S.tocov_matrix(nr=3, nt=256, dt=0.1)
S = Sj.tospecdata()
acfmat = S.tocov_matrix(nr=3, nt=256, dt=0.1)
vals = acfmat[:2, :] vals = acfmat[:2, :]
true_vals = np.array([[3.06073383, 0.0000000, -1.67748256, 0.], true_vals = np.array([[3.06073383, 0.0000000, -1.67748256, 0.],
[3.05235423, -0.1674357, -1.66811444, [3.05235423, -0.1674357, -1.66811444,
0.18693242]]) 0.18693242]])
self.assertTrue((np.abs(vals - true_vals) < 1e-7).all()) assert_array_almost_equal(vals, true_vals)
def test_tocovdata(self):
def test_tocovdata(): Nt = len(self.S.data) - 1
Sj = sm.Jonswap() acf = self.S.tocovdata(nr=0, nt=Nt)
S = Sj.tospecdata()
Nt = len(S.data) - 1
acf = S.tocovdata(nr=0, nt=Nt)
vals = acf.data[:5] vals = acf.data[:5]
true_vals = np.array( true_vals = np.array(
[3.06090339, 2.22658399, 0.45307391, -1.17495501, -2.05649042]) [3.06090339, 2.22658399, 0.45307391, -1.17495501, -2.05649042])
assert_array_almost_equal(vals, true_vals)
assert((np.abs(vals - true_vals) < 1e-6).all()) assert((np.abs(vals - true_vals) < 1e-6).all())
def test_to_t_pdf(self):
def test_to_t_pdf(): f = self.S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7, seed=100)
Sj = sm.Jonswap()
S = Sj.tospecdata()
f = S.to_t_pdf(pdef='Tc', paramt=(0, 10, 51), speed=7, seed=100)
vals = ['%2.3f' % val for val in f.data[:10]] vals = ['%2.3f' % val for val in f.data[:10]]
truevals = ['0.000', '0.014', '0.027', '0.040', truevals = ['0.000', '0.014', '0.027', '0.040',
'0.050', '0.059', '0.067', '0.073', '0.077', '0.082'] '0.050', '0.059', '0.067', '0.073', '0.077', '0.082']
@ -55,14 +52,9 @@ def test_to_t_pdf():
for t, v in zip(truevals, vals): for t, v in zip(truevals, vals):
assert(t == v) assert(t == v)
@slow @slow
def test_sim(): def test_sim(self):
Sj = sm.Jonswap() S = self.S
S = Sj.tospecdata()
#ns = 100
#dt = .2
#x1 = S.sim(ns, dt=dt)
import scipy.stats as st import scipy.stats as st
x2 = S.sim(20000, 20) x2 = S.sim(20000, 20)
@ -72,18 +64,12 @@ def test_sim():
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
assert(np.abs(m - trueval) < sa) assert(np.abs(m - trueval) < sa)
@slow @slow
def test_sim_nl(): def test_sim_nl(self):
S = self.S
Sj = sm.Jonswap()
S = Sj.tospecdata()
# ns = 100
# dt = .2
# x1 = S.sim_nl(ns, dt=dt)
import scipy.stats as st import scipy.stats as st
x2, _x1 = S.sim_nl(ns=20000, cases=40) x2, _x1 = S.sim_nl(ns=20000, cases=40)
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')
@ -100,27 +86,19 @@ def test_sim_nl():
# trueval, m, sa # trueval, m, sa
assert(np.abs(m - trueval) < 2 * sa) assert(np.abs(m - trueval) < 2 * sa)
def test_stats_nl(self):
def test_stats_nl(): S = self.S
Hs = 7.
Sj = sm.Jonswap(Hm0=Hs, Tp=11)
S = Sj.tospecdata()
me, va, sk, ku = S.stats_nl(moments='mvsk') me, va, sk, ku = S.stats_nl(moments='mvsk')
assert(me == 0.0) assert(me == 0.0)
assert_array_almost_equal(va, 3.0608203389019537) assert_array_almost_equal(va, 3.0608203389019537)
assert_array_almost_equal(sk, 0.18673120577479801) assert_array_almost_equal(sk, 0.18673120577479801)
assert_array_almost_equal(ku, 3.0619885212624176) assert_array_almost_equal(ku, 3.0619885212624176)
def test_testgaussian(self):
def test_testgaussian(): Hs = self.Sj.Hm0
S0 = self.S
Hs = 7
Sj = sm.Jonswap(Hm0=Hs)
S0 = Sj.tospecdata()
# ns =100; dt = .2 # ns =100; dt = .2
# x1 = S0.sim(ns, dt=dt) # x1 = S0.sim(ns, dt=dt)
S = S0.copy() S = S0.copy()
me, _va, sk, ku = S.stats_nl(moments='mvsk') me, _va, sk, ku = S.stats_nl(moments='mvsk')
S.tr = wtm.TrHermite( S.tr = wtm.TrHermite(
@ -132,93 +110,87 @@ def test_testgaussian():
assert(sum(t1 > t0) < 5) assert(sum(t1 > t0) < 5)
def test_moment(): class TestSpectrumHs5(unittest.TestCase):
Sj = sm.Jonswap(Hm0=5) def setUp(self):
S = Sj.tospecdata() # Make spectrum ob self.Sj = sm.Jonswap(Hm0=5.0)
self.S = self.Sj.tospecdata()
def test_moment(self):
S = self.S
vals, txt = S.moment() vals, txt = S.moment()
true_vals = [1.5614600345079888, 0.95567089481941048] true_vals = [1.5614600345079888, 0.95567089481941048]
true_txt = ['m0', 'm0tt'] true_txt = ['m0', 'm0tt']
for tv, v in zip(true_vals, vals):
assert_array_almost_equal(tv, v) assert_array_almost_equal(vals, true_vals)
for tv, v in zip(true_txt, txt): for tv, v in zip(true_txt, txt):
assert(tv == v) assert(tv == v)
def test_nyquist_freq(self):
S = self.S
assert_array_almost_equal(S.nyquist_freq(), 3.0)
def test_nyquist_freq(): def test_sampling_period(self):
Sj = sm.Jonswap(Hm0=5) S = self.S
S = Sj.tospecdata() # Make spectrum ob assert_array_almost_equal(S.sampling_period(), 1.0471975511965976)
assert(S.nyquist_freq() == 3.0)
def test_normalize(self):
def test_sampling_period(): S = self.S
Sj = sm.Jonswap(Hm0=5) mom, txt = S.moment(2)
S = Sj.tospecdata() # Make spectrum ob assert_array_almost_equal(mom,
assert(S.sampling_period() == 1.0471975511965976) [1.5614600345079888, 0.95567089481941048])
assert_array_equal(txt, ['m0', 'm0tt'])
def test_normalize():
Sj = sm.Jonswap(Hm0=5)
S = Sj.tospecdata() # Make spectrum ob
S.moment(2)
([1.5614600345079888, 0.95567089481941048], ['m0', 'm0tt'])
vals, _txt = S.moment(2) vals, _txt = S.moment(2)
true_vals = [1.5614600345079888, 0.95567089481941048] true_vals = [1.5614600345079888, 0.95567089481941048]
for tv, v in zip(true_vals, vals): assert_array_almost_equal(vals, true_vals)
assert_array_almost_equal(tv, v)
Sn = S.copy() Sn = S.copy()
Sn.normalize() Sn.normalize()
# Now the moments should be one # Now the moments should be one
new_vals, _txt = Sn.moment(2) new_vals, _txt = Sn.moment(2)
for v in new_vals: assert_array_almost_equal(new_vals, np.ones(2))
assert(np.abs(v - 1.0) < 1e-7)
def test_characteristic(self):
S = self.S
ch, R, txt = S.characteristic(1)
assert_array_almost_equal(ch, 8.59007646)
assert_array_almost_equal(R, 0.03040216)
self.assert_(txt == ['Tm01'])
def test_characteristic(): ch, R, txt = S.characteristic([1, 2, 3]) # fact a vector of integers
''' assert_array_almost_equal(ch, [8.59007646, 8.03139757, 5.62484314])
>>> import wafo.spectrum.models as sm assert_array_almost_equal(R,
>>> Sj = sm.Jonswap(Hm0=5) [[0.03040216, 0.02834263, NAN],
>>> S = Sj.tospecdata() #Make spectrum ob [0.02834263, 0.0274645, NAN],
>>> S.characteristic(1) [NAN, NAN, 0.01500249]])
(array([ 8.59007646]), array([[ 0.03040216]]), ['Tm01']) assert_array_equal(txt, ['Tm01', 'Tm02', 'Tm24'])
>>> [ch, R, txt] = S.characteristic([1,2,3]) # fact a vector of integers ch, R, txt = S.characteristic('Ss') # fact a string
>>> ch; R; txt assert_array_almost_equal(ch, [0.04963112])
array([ 8.59007646, 8.03139757, 5.62484314]) assert_array_almost_equal(R, [[2.63624782e-06]])
array([[ 0.03040216, 0.02834263, nan], assert_array_equal(txt, ['Ss'])
[ 0.02834263, 0.0274645 , nan],
[ nan, nan, 0.01500249]])
['Tm01', 'Tm02', 'Tm24']
>>> S.characteristic('Ss') # fact a string # fact a list of strings
(array([ 0.04963112]), array([[ 2.63624782e-06]]), ['Ss']) ch, R, txt = S.characteristic(['Hm0', 'Tm02'])
>>> S.characteristic(['Hm0','Tm02']) # fact a list of strings assert_array_almost_equal(ch,
(array([ 4.99833578, 8.03139757]), array([[ 0.05292989, 0.02511371], [4.99833578, 8.03139757])
[ 0.02511371, 0.0274645 ]]), ['Hm0', 'Tm02']) assert_array_almost_equal(R, [[0.05292989, 0.02511371],
''' [0.02511371, 0.0274645]])
assert_array_equal(txt, ['Hm0', 'Tm02'])
def test_bandwidth(): class TestSpectrumHs3(unittest.TestCase):
def test_bandwidth(self):
Sj = sm.Jonswap(Hm0=3, Tp=7) Sj = sm.Jonswap(Hm0=3, Tp=7)
w = np.linspace(0, 4, 256) w = np.linspace(0, 4, 256)
S = SpecData1D(Sj(w), w) # Make spectrum object from numerical values S = SpecData1D(Sj(w), w) # Make spectrum object from numerical values
vals = S.bandwidth([0, 1, 2, 3]) vals = S.bandwidth([0, 1, 2, 3])
true_vals = np.array([0.73062845, 0.34476034, 0.68277527, 2.90817052]) true_vals = [0.73062845, 0.34476034, 0.68277527, 2.90817052]
assert((np.abs(vals - true_vals) < 1e-7).all()) assert_array_almost_equal(vals, true_vals)
def test_docstrings():
import doctest
doctest.testmod()
if __name__ == '__main__': if __name__ == '__main__':
import nose import nose
nose.run() nose.run()
# test_docstrings()
# test_tocovdata()
# test_tocovmatrix()
# test_sim()
# test_bandwidth()

23
wafo/stats/_distn_infrastructure.py
Unescape Escape View File

@ -1,9 +1,11 @@
from scipy.stats._distn_infrastructure import * from __future__ import absolute_import
from scipy.stats._distn_infrastructure import (_skew, _kurtosis, # @UnresolvedImport from scipy.stats._distn_infrastructure import * # @UnusedWildImport
_lazywhere, _ncx2_log_pdf, _ncx2_pdf, _ncx2_cdf) from scipy.stats._distn_infrastructure import (_skew, # @UnusedImport
from wafo.stats.estimation import FitDistribution _kurtosis, _lazywhere, _ncx2_log_pdf, # @IgnorePep8 @UnusedImport
from wafo.stats._constants import _EPS, _XMAX _ncx2_pdf, _ncx2_cdf) # @UnusedImport @IgnorePep8
import numpy as np from .estimation import FitDistribution
from ._constants import _XMAX
_doc_default_example = """\ _doc_default_example = """\
Examples Examples
@ -326,10 +328,11 @@ def nlogps(self, theta, x):
return T return T
def _reduce_func(self, args, kwds): def _reduce_func(self, args, options):
# First of all, convert fshapes params to fnum: eg for stats.beta, # First of all, convert fshapes params to fnum: eg for stats.beta,
# shapes='a, b'. To fix `a`, can specify either `f1` or `fa`. # shapes='a, b'. To fix `a`, can specify either `f1` or `fa`.
# Convert the latter into the former. # Convert the latter into the former.
kwds = options.copy()
if self.shapes: if self.shapes:
shapes = self.shapes.replace(',', ' ').split() shapes = self.shapes.replace(',', ' ').split()
for j, s in enumerate(shapes): for j, s in enumerate(shapes):
@ -384,9 +387,9 @@ def _reduce_func(self, args, kwds):
return x0, func, restore, args return x0, func, restore, args
def fit(self, data, *args, **kwds): def fit(self, data, *args, **kwargs):
""" """
Return ML or MPS estimate for shape, location, and scale parameters from data. Return ML/MPS estimate for shape, location, and scale parameters from data.
ML and MPS stands for Maximum Likelihood and Maximum Product Spacing, ML and MPS stands for Maximum Likelihood and Maximum Product Spacing,
respectively. Starting estimates for respectively. Starting estimates for
@ -476,6 +479,7 @@ def fit(self, data, *args, **kwds):
if Narg > self.numargs: if Narg > self.numargs:
raise TypeError("Too many input arguments.") raise TypeError("Too many input arguments.")
kwds = kwargs.copy()
start = [None]*2 start = [None]*2
if (Narg < self.numargs) or not ('loc' in kwds and if (Narg < self.numargs) or not ('loc' in kwds and
'scale' in kwds): 'scale' in kwds):
@ -573,4 +577,3 @@ rv_continuous.nlogps = nlogps
rv_continuous._reduce_func = _reduce_func rv_continuous._reduce_func = _reduce_func
rv_continuous.fit = fit rv_continuous.fit = fit
rv_continuous.fit2 = fit2 rv_continuous.fit2 = fit2

1
wafo/tests/__init__.py
Unescape Escape View File

@ -1 +0,0 @@

2
wafo/tests/conftest.py
Unescape Escape View File

@ -9,4 +9,4 @@
""" """
from __future__ import print_function, absolute_import, division from __future__ import print_function, absolute_import, division
import pytest import pytest # @UnusedImport

63
wafo/tests/test_integrate.py
Unescape Escape View File

@ -7,10 +7,23 @@ import unittest
import numpy as np import numpy as np
from numpy import exp, Inf from numpy import exp, Inf
from numpy.testing import assert_array_almost_equal from numpy.testing import assert_array_almost_equal
from wafo.integrate import gaussq from wafo.integrate import gaussq, quadgr, clencurt, romberg
class Gaussq(unittest.TestCase): class TestIntegrators(unittest.TestCase):
def test_clencurt(self):
val, err = clencurt(np.exp, 0, 2)
assert_array_almost_equal(val, np.expm1(2))
self.assert_(err < 1e-10)
def test_romberg(self):
tol = 1e-7
q, err = romberg(np.sqrt, 0, 10, 0, abseps=tol)
assert_array_almost_equal(q, 2.0/3 * 10**(3./2))
self.assert_(err < tol)
class TestGaussq(unittest.TestCase):
''' '''
1 : p(x) = 1 a =-1, b = 1 Gauss-Legendre 1 : p(x) = 1 a =-1, b = 1 Gauss-Legendre
2 : p(x) = exp(-x^2) a =-inf, b = inf Hermite 2 : p(x) = exp(-x^2) a =-inf, b = inf Hermite
@ -60,6 +73,52 @@ class Gaussq(unittest.TestCase):
val, _err = gaussq(lambda x: x, 0, 1, wfun=9) val, _err = gaussq(lambda x: x, 0, 1, wfun=9)
assert_array_almost_equal(val, 0.26666667) assert_array_almost_equal(val, 0.26666667)
class TestQuadgr(unittest.TestCase):
def test_log(self):
Q, err = quadgr(np.log, 0, 1)
assert_array_almost_equal(Q, -1)
self.assert_(err < 1e-5)
def test_exp(self):
Q, err = quadgr(np.exp, 0, 9999*1j*np.pi)
assert_array_almost_equal(Q, -2.0000000000122662)
self.assert_(err < 1.0e-8)
def test_integral3(self):
tol = 1e-12
Q, err = quadgr(lambda x: np.sqrt(4-x**2), 0, 2, tol)
assert_array_almost_equal(Q, np.pi)
self.assert_(err < tol)
# (3.1415926535897811, 1.5809575870662229e-13)
def test_integral4(self):
Q, err = quadgr(lambda x: 1./x**0.75, 0, 1)
assert_array_almost_equal(Q, 4)
self.assert_(err < 1.0e-12)
def test_integrand4(self):
tol = 1e-10
Q, err = quadgr(lambda x: 1./np.sqrt(1-x**2), -1, 1, tol)
assert_array_almost_equal(Q, np.pi)
self.assert_(err < tol)
# (3.141596056985029, 6.2146261559092864e-06)
def test_integrand5(self):
tol = 1e-9
Q, err = quadgr(lambda x: np.exp(-x**2), -np.inf, np.inf, tol)
assert_array_almost_equal(Q, np.sqrt(np.pi))
self.assert_(err < tol)
# (1.7724538509055152, 1.9722334876348668e-11)
def test_integrand6(self):
tol = 1e-9
Q, err = quadgr(lambda x: np.cos(x)*np.exp(-x), 0, np.inf, tol)
assert_array_almost_equal(Q, 0.5)
self.assert_(err < tol)
# (0.50000000000000044, 7.3296813063450372e-11)
if __name__ == "__main__": if __name__ == "__main__":
# import sys;sys.argv = ['', 'Test.testName'] # import sys;sys.argv = ['', 'Test.testName']
unittest.main() unittest.main()

4
wafo/transform/models.py
Unescape Escape View File

@ -478,8 +478,8 @@ class TrOchi(TrCommon2):
''' '''
Returns ga, gb, sigma2, mean2 Returns ga, gb, sigma2, mean2
''' '''
if (self._phat is None or self.sigma != self._phat[0] if (self._phat is None or self.sigma != self._phat[0] or
or self.mean != self._phat[1]): self.mean != self._phat[1]):
self._par_from_stats() self._par_from_stats()
# sigma1 = self._phat[0] # sigma1 = self._phat[0]
# mean1 = self._phat[1] # mean1 = self._phat[1]

156
wafo/wave_theory/core.py
Unescape Escape View File

@ -7,9 +7,11 @@ from __future__ import absolute_import
import numpy as np import numpy as np
from numpy import exp, expm1, inf, nan, pi, hstack, where, atleast_1d, cos, sin from numpy import exp, expm1, inf, nan, pi, hstack, where, atleast_1d, cos, sin
from .dispersion_relation import w2k, k2w # @UnusedImport from .dispersion_relation import w2k, k2w # @UnusedImport
from ..misc import JITImport
__all__ = ['w2k', 'k2w', 'sensor_typeid', 'sensor_type', 'TransferFunction'] __all__ = ['w2k', 'k2w', 'sensor_typeid', 'sensor_type', 'TransferFunction']
_MODELS = JITImport('wafo.spectrum.models')
def hyperbolic_ratio(a, b, sa, sb): def hyperbolic_ratio(a, b, sa, sb):
''' '''
@ -349,7 +351,8 @@ class TransferFunction(object):
Gwt = -Gwt Gwt = -Gwt
return Hw, Gwt return Hw, Gwt
__call__ = tran __call__ = tran
#---Private member methods
# --- Private member methods ---
def _get_ee_cthxy(self, theta, kw): def _get_ee_cthxy(self, theta, kw):
# convert from angle in degrees to radians # convert from angle in degrees to radians
@ -387,7 +390,7 @@ class TransferFunction(object):
ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw) ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
return np.ones_like(w), ee return np.ones_like(w), ee
#---- Vertical surface velocity and acceleration----- # --- Vertical surface velocity and acceleration ---
def _n_t(self, w, theta, kw): def _n_t(self, w, theta, kw):
''' n_t = Eta_t ''' ''' n_t = Eta_t '''
ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw) ee, unused_cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
@ -434,7 +437,7 @@ class TransferFunction(object):
# hyperbolic_ratio = cosh(zk)/cosh(hk) # hyperbolic_ratio = cosh(zk)/cosh(hk)
return self.rho * self.g * hyperbolic_ratio(zk, hk, 1, 1), ee return self.rho * self.g * hyperbolic_ratio(zk, hk, 1, 1), ee
#---- Water particle velocities --- # --- Water particle velocities ---
def _u(self, w, theta, kw): def _u(self, w, theta, kw):
''' U = x-velocity''' ''' U = x-velocity'''
ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw) ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
@ -459,7 +462,7 @@ class TransferFunction(object):
# w*sinh(zk)/sinh(hk), -? # w*sinh(zk)/sinh(hk), -?
return w * hyperbolic_ratio(zk, hk, -1, -1), -1j * ee return w * hyperbolic_ratio(zk, hk, -1, -1), -1j * ee
#---- Water particle acceleration --- # --- Water particle acceleration ---
def _u_t(self, w, theta, kw): def _u_t(self, w, theta, kw):
''' U_t = x-acceleration''' ''' U_t = x-acceleration'''
ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw) ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
@ -484,7 +487,7 @@ class TransferFunction(object):
# w*sinh(zk)/sinh(hk), ? # w*sinh(zk)/sinh(hk), ?
return (w ** 2) * hyperbolic_ratio(zk, hk, -1, -1), -ee return (w ** 2) * hyperbolic_ratio(zk, hk, -1, -1), -ee
#---- Water particle displacement --- # --- Water particle displacement ---
def _x_p(self, w, theta, kw): def _x_p(self, w, theta, kw):
''' X_p = x-displacement''' ''' X_p = x-displacement'''
ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw) ee, cthx, unused_cthy = self._get_ee_cthxy(theta, kw)
@ -508,88 +511,73 @@ class TransferFunction(object):
zk = self._get_zk(kw) zk = self._get_zk(kw)
return hyperbolic_ratio(zk, hk, -1, -1), ee # sinh(zk)./sinh(hk), ee return hyperbolic_ratio(zk, hk, -1, -1), ee # sinh(zk)./sinh(hk), ee
# def wave_pressure(z, Hm0, h=10000, g=9.81, rho=1028):
# ''' def wave_pressure(z, Hm0, h=10000, g=9.81, rho=1028):
# Calculate pressure amplitude due to water waves. '''
# Calculate pressure amplitude due to water waves.
# Parameters
# ---------- Parameters
# z : array-like ----------
# depth where pressure is calculated [m] z : array-like
# Hm0 : array-like depth where pressure is calculated [m]
# significant wave height (same as the average of the 1/3'rd highest Hm0 : array-like
# waves in a seastate. [m] significant wave height (same as the average of the 1/3'rd highest
# h : real scalar waves in a seastate. [m]
# waterdepth (default 10000 [m]) h : real scalar
# g : real scalar waterdepth (default 10000 [m])
# acceleration of gravity (default 9.81 m/s**2) g : real scalar
# rho : real scalar acceleration of gravity (default 9.81 m/s**2)
# water density (default 1028 kg/m**3) rho : real scalar
# water density (default 1028 kg/m**3)
#
# Returns
# ------- Returns
# p : ndarray -------
# pressure amplitude due to water waves at water depth z. [Pa] p : ndarray
# pressure amplitude due to water waves at water depth z. [Pa]
# PRESSURE calculate pressure amplitude due to water waves according to
# linear theory. PRESSURE calculate pressure amplitude due to water waves according to
# linear theory.
# Example
# ----- Example
# >>> import pylab as plt -----
# >>> z = -np.linspace(10,20) >>> import pylab as plt
# >>> fh = plt.plot(z, wave_pressure(z, Hm0=1, h=20)) >>> z = -np.linspace(10,20)
# >>> plt.show() >>> fh = plt.plot(z, wave_pressure(z, Hm0=1, h=20))
# >>> plt.show()
# See also
# -------- See also
# w2k --------
# w2k
#
# u = psweep.Fn*sqrt(mgf.length*9.81) '''
# z = -10; h = inf;
# Hm0 = 1.5;Tp = 4*sqrt(Hm0);
# S = jonswap([],[Hm0,Tp]);
# Hw = tran(S.w,0,[0 0 -z],'P',h)
# Sm = S;
# Sm.S = Hw.'.*S.S;
# x1 = spec2sdat(Sm,1000);
# pwave = pressure(z,Hm0,h)
#
# plot(psweep.x{1}/u, psweep.f)
# hold on
# plot(x1(1:100,1)-30,x1(1:100,2),'r')
# '''
#
#
# Assume seastate with jonswap spectrum: # Assume seastate with jonswap spectrum:
# Tp = 4 * np.sqrt(Hm0)
# Tp = 4 * np.sqrt(Hm0) gam = _MODELS.jonswap_peakfact(Hm0, Tp)
# gam = jonswap_peakfact(Hm0, Tp) Tm02 = Tp / (1.30301 - 0.01698 * gam + 0.12102 / gam)
# Tm02 = Tp / (1.30301 - 0.01698 * gam + 0.12102 / gam) w = 2 * np.pi / Tm02
# w = 2 * np.pi / Tm02 kw, unused_kw2 = w2k(w, 0, h)
# kw, unused_kw2 = w2k(w, 0, h)
# hk = kw * h
# hk = kw * h zk1 = kw * z
# zk1 = kw * z zk = hk + zk1 # z measured positive upward from mean water level (default)
# zk = hk + zk1 # z measured positive upward from mean water level (default) # zk = hk-zk1 # z measured positive downward from mean water level
# zk = hk-zk1; % z measured positive downward from mean water level # zk1 = -zk1
# zk1 = -zk1; # zk = zk1 # z measured positive upward from sea floor
# zk = zk1; % z measured positive upward from sea floor
#
# cosh(zk)/cosh(hk) approx exp(zk) for large h # cosh(zk)/cosh(hk) approx exp(zk) for large h
# hyperbolic_ratio(zk,hk,1,1) = cosh(zk)/cosh(hk) # hyperbolic_ratio(zk,hk,1,1) = cosh(zk)/cosh(hk)
# pr = np.where(np.pi < hk, np.exp(zk1), hyperbolic_ratio(zk, hk, 1, 1)) # pr = np.where(np.pi < hk, np.exp(zk1), hyperbolic_ratio(zk, hk, 1, 1))
# pr = hyperbolic_ratio(zk, hk, 1, 1) pr = hyperbolic_ratio(zk, hk, 1, 1)
# pressure = (rho * g * Hm0 / 2) * pr pressure = (rho * g * Hm0 / 2) * pr
#
## pos = [np.zeros_like(z),np.zeros_like(z),z] # pos = [np.zeros_like(z),np.zeros_like(z),z]
## tf = TransferFunction(pos=pos, sensortype='p', h=h, rho=rho, g=g) # tf = TransferFunction(pos=pos, sensortype='p', h=h, rho=rho, g=g)
## Hw, Gwt = tf.tran(w,0) # Hw, Gwt = tf.tran(w,0)
## pressure2 = np.abs(Hw) * Hm0 / 2 # pressure2 = np.abs(Hw) * Hm0 / 2
#
# return pressure return pressure
def test_docstrings(): def test_docstrings():

Write Preview
Loading…
Cancel
Save