src/wafo/_misc_numba.py
"""
Created on 6. okt. 2016
@author: pab
"""
from numba import jit, float64, int64, int32, int8, void
import numpy as np
# @guvectorize(['void(int64[:], int8[:], int64[:])'], '(n),(n)->(), (), ()')
# def _find_first_cross(ind, y, ix, dcross, start):
# ix, dcross, start, v = 0, 0, 0, 0
# n = len(y)
# if y[0] < v:
# dcross = -1 # first is a up-crossing
# elif y[0] > v:
# dcross = 1 # first is a down-crossing
# elif y[0] == v:
# # Find out what type of crossing we have next time..
# for i in range(1, n):
# start = i
# if y[i] < v:
# ind[ix] = i - 1 # first crossing is a down crossing
# ix += 1
# dcross = -1 # The next crossing is a up-crossing
# break
# elif y[i] > v:
# ind[ix] = i - 1 # first crossing is a up-crossing
# ix += 1
# dcross = 1 # The next crossing is a down-crossing
# break
@jit(int64(int64[:], int8[:]), nopython=True)
def _findcross(ind, y):
"""Returns indices to zero level crossings of y vector
Notes
-----
Same implementation as findcross function found in c_functions.c.
"""
ix, dcross, start, v = 0, 0, 0, 0
n = len(y)
if y[0] < v:
dcross = -1 # first is a up-crossing
elif y[0] > v:
dcross = 1 # first is a down-crossing
elif y[0] == v:
# Find out what type of crossing we have next time..
for i in range(1, n):
start = i
if y[i] < v:
ind[ix] = i - 1 # first crossing is a down crossing
ix += 1
dcross = -1 # The next crossing is a up-crossing
break
elif y[i] > v:
ind[ix] = i - 1 # first crossing is a up-crossing
ix += 1
dcross = 1 # The next crossing is a down-crossing
break
for i in range(start, n - 1):
if ((dcross == -1 and y[i] <= v and v < y[i + 1])
or (dcross == 1 and v <= y[i] and y[i + 1] < v)):
ind[ix] = i
ix += 1
dcross = -dcross
return ix
def findcross(xn):
"""Returns indices to zero up- and down-crossings of a vector"""
ind = np.empty(len(xn), dtype=np.int64)
m = _findcross(ind, xn)
return ind[:m]
def _make_findrfc(cmp1, cmp2):
@jit(int64(int64[:], float64[:], float64), nopython=True)
def local_findrfc(t, y, h):
"""Returns indices, t, to RFC turningpoints of a vector y of turningpoints
Notes
-----
Same implementation as rfcfilter function found in rfcfilter.m in wafo matlab.
"""
# cmp1, cmp2 = (a_le_b, a_lt_b) if method==0 else (a_lt_b, a_le_b)
n = len(y)
j, t0, z0 = 0, 0, 0
y0 = y[t0]
# The rainflow filter
for ti in range(1, n):
fpi = y0 + h
fmi = y0 - h
yi = y[ti]
if z0 == 0:
if cmp1(yi, fmi):
z1 = -1
elif cmp1(fpi, yi):
z1 = +1
else:
z1 = 0
t1, y1 = (t0, y0) if z1 == 0 else (ti, yi)
else:
if (((z0 == +1) and cmp1(yi, fmi))
or ((z0 == -1) and cmp2(yi, fpi))):
z1 = -1
elif (((z0 == +1) and cmp2(fmi, yi)) or
((z0 == -1) and cmp1(fpi, yi))):
z1 = +1
else:
raise ValueError
# warnings.warn('Something wrong, i={}'.format(tim1))
# Update y1
if z1 != z0:
t1, y1 = ti, yi
elif z1 == -1:
t1, y1 = (t0, y0) if y0 < yi else (ti, yi)
elif z1 == +1:
t1, y1 = (t0, y0) if y0 > yi else (ti, yi)
# Update y if y0 is a turning point
if abs(z0 - z1) == 2:
j += 1
t[j] = t0
# Update t0, y0, z0
t0, y0, z0 = t1, y1, z1
# end
# Update y if last y0 is greater than (or equal) threshold
if cmp2(h, abs(y0 - y[t[j]])):
j += 1
t[j] = t0
return j + 1
return local_findrfc
@jit(int32(float64, float64), nopython=True)
def a_le_b(a, b):
return a <= b
@jit(int32(float64, float64), nopython=True)
def a_lt_b(a, b):
return a < b
_findrfc_le = _make_findrfc(a_le_b, a_lt_b)
_findrfc_lt = _make_findrfc(a_lt_b, a_le_b)
@jit(int64(int64[:], float64[:], float64), nopython=True)
def _findrfc(ind, y, h):
"""Returns indices to RFC turningpoints of a vector y of turningpoints
Notes
-----
Same implementation as findrfc function found in c_functions.c.
"""
n = len(y)
t_start = 0
nc = n // 2
ix = 0
for i in range(nc):
Tmi = t_start + 2 * i
Tpl = t_start + 2 * i + 2
xminus = y[2 * i]
xplus = y[2 * i + 2]
if(i != 0):
j = i - 1
while ((j >= 0) and (y[2 * j + 1] <= y[2 * i + 1])):
if (y[2 * j] < xminus):
xminus = y[2 * j]
Tmi = t_start + 2 * j
j -= 1
if (xminus >= xplus):
if (y[2 * i + 1] - xminus >= h):
ind[ix] = Tmi
ix += 1
ind[ix] = (t_start + 2 * i + 1)
ix += 1
# goto L180 continue
else:
j = i + 1
while (j < nc):
if (y[2 * j + 1] >= y[2 * i + 1]):
break # goto L170
if((y[2 * j + 2] <= xplus)):
xplus = y[2 * j + 2]
Tpl = (t_start + 2 * j + 2)
j += 1
else:
if ((y[2 * i + 1] - xminus) >= h):
ind[ix] = Tmi
ix += 1
ind[ix] = (t_start + 2 * i + 1)
ix += 1
# iy = i
continue
# goto L180
# L170:
if (xplus <= xminus):
if ((y[2 * i + 1] - xminus) >= h):
ind[ix] = Tmi
ix += 1
ind[ix] = (t_start + 2 * i + 1)
ix += 1
elif ((y[2 * i + 1] - xplus) >= h):
ind[ix] = (t_start + 2 * i + 1)
ix += 1
ind[ix] = Tpl
ix += 1
# L180:
# iy=i
# /* for i */
return ix
def findrfc(y, h, method=0):
"""Returns indices to RFC turningpoints of a vector y of turningpoints
Parameters
----------
method : scalar integer
0: removes cycles with range < h. (default)
1: removes cycles with range <= h.
2: removes cycles with range < h alternative version
Notes
-----
Method 0 and 1 is the same implementation as found in rfcfilter.m in matlab.
Method 2 is the same implementation as found in findrfc function found in c_functions.c.
"""
n = len(y)
t = np.zeros(n, dtype=np.int64)
findrfc_ = [_findrfc_le, _findrfc_lt, _findrfc][method]
m = findrfc_(t, y, h)
return t[:m]
@jit(void(float64[:], float64[:], float64[:], float64[:],
float64[:], float64[:], float64, float64,
int32, int32, int32, int32), nopython=True)
def _finite_water_disufq(rvec, ivec, rA, iA, w, kw, h, g, nmin, nmax, m, n):
# kfact is set to 2 in order to exploit the symmetry.
# If you set kfact to 1, you must uncomment all statements
# including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2].
kfact = 2.0
for ix in range(nmin - 1, nmax):
# for (ix = nmin-1;ix<nmax;ix++) {
kw1 = kw[ix]
w1 = w[ix]
tmp1 = np.tanh(kw1 * h)
# Cg, wave group velocity
Cg = 0.5 * g * (tmp1 + kw1 * h * (1.0 - tmp1 * tmp1)) / w1 # OK
tmp1 = 0.5 * g * (kw1 / w1) * (kw1 / w1)
tmp2 = 0.5 * w1 * w1 / g
tmp3 = g * kw1 / (w1 * Cg)
tmp4 = kw1 / np.sinh(2.0 * kw1 * h) if kw1 * h < 300.0 else 0.0
# Difference frequency effects finite water depth
Edij = (tmp1 - tmp2 + tmp3) / (1.0 - g * h / (Cg * Cg)) - tmp4 # OK
# Sum frequency effects finite water depth
Epij = (3.0 * (tmp1 - tmp2) /
(1.0 - tmp1 / kw1 * np.tanh(2.0 * kw1 * h)) +
3.0 * tmp2 - tmp1) # OK
# printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix * m
iz1 = 2 * ixi
# iz2 = n*m-ixi;
for i in range(m):
rrA = rA[ixi] * rA[ixi]
iiA = iA[ixi] * iA[ixi]
riA = rA[ixi] * iA[ixi]
# Sum frequency effects along the diagonal
rvec[iz1] += kfact * (rrA - iiA) * Epij
ivec[iz1] += kfact * 2.0 * riA * Epij
# rvec[iz2] += kfact*(rrA-iiA)*Epij;
# ivec[iz2] -= kfact*2.0*riA*Epij;
# iz2++;
# Difference frequency effects along the diagonal
# are only contributing to the mean
rvec[i] += 2.0 * (rrA + iiA) * Edij
ixi += 1
iz1 += 1
# }
for jy in range(ix + 1, nmax):
# w1 = w[ix];
# kw1 = kw[ix];
w2 = w[jy]
kw2 = kw[jy]
tmp1 = g * (kw1 / w1) * (kw2 / w2)
tmp2 = 0.5 / g * (w1 * w1 + w2 * w2 + w1 * w2)
tmp3 = 0.5 * g * \
(w1 * kw2 * kw2 + w2 * kw1 * kw1) / (w1 * w2 * (w1 + w2))
tmp4 = (1 - g * (kw1 + kw2) / (w1 + w2) / (w1 + w2) *
np.tanh((kw1 + kw2) * h))
Epij = (tmp1 - tmp2 + tmp3) / tmp4 + tmp2 - 0.5 * tmp1 # OK */
tmp2 = 0.5 / g * (w1 * w1 + w2 * w2 - w1 * w2) # OK*/
tmp3 = -0.5 * g * \
(w1 * kw2 * kw2 - w2 * kw1 * kw1) / (w1 * w2 * (w1 - w2))
tmp4 = (1.0 - g * (kw1 - kw2) / (w1 - w2) /
(w1 - w2) * np.tanh((kw1 - kw2) * h))
Edij = (tmp1 - tmp2 + tmp3) / tmp4 + tmp2 - 0.5 * tmp1 # OK */
# printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix * m
jyi = jy * m
iz1 = ixi + jyi
iv1 = jyi - ixi
# iz2 = (n*m-iz1)
# iv2 = n*m-iv1
for i in range(m):
# for (i=0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi] * rA[jyi] # rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi] * iA[jyi] # iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi] * iA[jyi] # riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi] * rA[jyi] # irA = iA[i][ix]*rA[i][jy];
# Sum frequency effects */
tmp1 = kfact * 2.0 * (rrA - iiA) * Epij
tmp2 = kfact * 2.0 * (riA + irA) * Epij
rvec[iz1] += tmp1 # rvec[i][jy+ix] += tmp1
ivec[iz1] += tmp2 # ivec[i][jy+ix] += tmp2
# rvec[iz2] += tmp1 # rvec[i][n*m-(jy+ix)] += tmp1
# ivec[iz2] -= tmp2 # ivec[i][n*m-(jy+ix)] -= tmp2
# iz2++
# Difference frequency effects */
tmp1 = kfact * 2.0 * (rrA + iiA) * Edij
tmp2 = kfact * 2.0 * (riA - irA) * Edij
rvec[iv1] += tmp1 # rvec[i][jy-ix] += tmp1
ivec[iv1] += tmp2 # ivec[i][jy-ix] -= tmp2
# rvec[iv2] += tmp1
# ivec[iv2] -= tmp2
# iv2 += 1
ixi += 1
jyi += 1
iz1 += 1
iv1 += 1
@jit(void(float64[:], float64[:], float64[:], float64[:],
float64[:], float64[:], float64, float64,
int32, int32, int32, int32), nopython=True)
def _deep_water_disufq(rvec, ivec, rA, iA, w, kw, h, g, nmin, nmax, m, n):
# kfact is set to 2 in order to exploit the symmetry.
# If you set kfact to 1, you must uncomment all statements
# including the expressions: rvec[iz2], rvec[iv2], ivec[iz2] and ivec[iv2].
kfact = 2.0
for ix in range(nmin - 1, nmax):
ixi = ix * m
iz1 = 2 * ixi
iz2 = n * m - ixi
kw1 = kw[ix]
Epij = kw1
for _i in range(m):
rrA = rA[ixi] * rA[ixi]
iiA = iA[ixi] * iA[ixi]
riA = rA[ixi] * iA[ixi]
# Sum frequency effects along the diagonal
tmp1 = kfact * (rrA - iiA) * Epij
tmp2 = kfact * 2.0 * riA * Epij
rvec[iz1] += tmp1
ivec[iz1] += tmp2
ixi += 1
iz1 += 1
# rvec[iz2] += tmp1
# ivec[iz2] -= tmp2
iz2 += 1
# Difference frequency effects are zero along the diagonal
# and are thus not contributing to the mean.
for jy in range(ix + 1, nmax):
kw2 = kw[jy]
Epij = 0.5 * (kw2 + kw1)
Edij = -0.5 * (kw2 - kw1)
# printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix * m
jyi = jy * m
iz1 = ixi + jyi
iv1 = jyi - ixi
iz2 = (n * m - iz1)
iv2 = (n * m - iv1)
for _i in range(m):
rrA = rA[ixi] * rA[jyi] # rrA = rA[i][ix]*rA[i][jy]
iiA = iA[ixi] * iA[jyi] # iiA = iA[i][ix]*iA[i][jy]
riA = rA[ixi] * iA[jyi] # riA = rA[i][ix]*iA[i][jy]
irA = iA[ixi] * rA[jyi] # irA = iA[i][ix]*rA[i][jy]
# Sum frequency effects
tmp1 = kfact * 2.0 * (rrA - iiA) * Epij
tmp2 = kfact * 2.0 * (riA + irA) * Epij
rvec[iz1] += tmp1 # rvec[i][ix+jy] += tmp1
ivec[iz1] += tmp2 # ivec[i][ix+jy] += tmp2
# rvec[iz2] += tmp1 # rvec[i][n*m-(ix+jy)] += tmp1
# ivec[iz2] -= tmp2 # ivec[i][n*m-(ix+jy)] -= tmp2
iz2 += 1
# Difference frequency effects */
tmp1 = kfact * 2.0 * (rrA + iiA) * Edij
tmp2 = kfact * 2.0 * (riA - irA) * Edij
rvec[iv1] += tmp1 # rvec[i][jy-ix] += tmp1
ivec[iv1] += tmp2 # ivec[i][jy-ix] += tmp2
# rvec[iv2] += tmp1 # rvec[i][n*m-(jy-ix)] += tmp1
# ivec[iv2] -= tmp2 # ivec[i][n*m-(jy-ix)] -= tmp2
iv2 += 1
ixi += 1
jyi += 1
iz1 += 1
iv1 += 1
def disufq(rA, iA, w, kw, h, g, nmin, nmax, m, n):
"""
DISUFQ Is an internal function to spec2nlsdat
Parameters
----------
rA, iA = real and imaginary parts of the amplitudes (size m X n).
w = vector with angular frequencies (w>=0)
kw = vector with wavenumbers (kw>=0)
h = water depth (h >=0)
g = constant acceleration of gravity
nmin = minimum index where rA(:,nmin) and iA(:,nmin) is
greater than zero.
nmax = maximum index where rA(:,nmax) and iA(:,nmax) is
greater than zero.
m = size(rA,1),size(iA,1)
n = size(rA,2),size(iA,2), or size(rvec,2),size(ivec,2)
returns
-------
rvec, ivec = real and imaginary parts of the resultant (size m X n).
DISUFQ returns the summation of difference frequency and sum
frequency effects in the vector vec = rvec +sqrt(-1)*ivec.
The 2'nd order contribution to the Stokes wave is then calculated by
a simple 1D Fourier transform, real(FFT(vec)).
"""
rvec = np.zeros(n * m)
ivec = np.zeros(n * m)
if h > 10000: # { /* deep water /Inifinite water depth */
_deep_water_disufq(rvec, ivec, rA, iA, w, kw, h, g, nmin, nmax, m, n)
else:
_finite_water_disufq(rvec, ivec, rA, iA, w, kw, h, g, nmin, nmax, m, n)
return rvec, ivec
@jit(int64(float64[:], float64[:], float64[:, :]), nopython=True)
def _findrfc3_astm(array_ext, a, array_out):
"""
Rain flow without time analysis
Return [ampl ampl_mean nr_of_cycle]
By Adam Nieslony
Visit the MATLAB Central File Exchange for latest version
http://www.mathworks.com/matlabcentral/fileexchange/3026
"""
n = len(array_ext)
po = 0
# The original rainflow counting by Nieslony, unchanged
j = -1
c_nr1 = 1
for i in range(n):
j += 1
a[j] = array_ext[i]
while j >= 2 and abs(a[j - 1] - a[j - 2]) <= abs(a[j] - a[j - 1]):
ampl = abs((a[j - 1] - a[j - 2]) / 2)
mean = (a[j - 1] + a[j - 2]) / 2
if j == 2:
a[0] = a[1]
a[1] = a[2]
j = 1
if (ampl > 0):
array_out[po, :] = (ampl, mean, 0.5)
po += 1
else:
a[j - 2] = a[j]
j = j - 2
if (ampl > 0):
array_out[po, :] = (ampl, mean, 1.0)
po += 1
c_nr1 += 1
c_nr2 = 1
for i in range(j):
ampl = abs(a[i] - a[i + 1]) / 2
mean = (a[i] + a[i + 1]) / 2
if (ampl > 0):
array_out[po, :] = (ampl, mean, 0.5)
po += 1
c_nr2 += 1
return c_nr1 #, c_nr2]
@jit(int64(float64[:], float64[:], float64[:], float64[:], float64[:, :]),
nopython=True)
def _findrfc5_astm(array_ext, array_t, a, t, array_out):
"""
Rain flow with time analysis
returns
[ampl ampl_mean nr_of_cycle cycle_begin_time cycle_period_time]
By Adam Nieslony
Visit the MATLAB Central File Exchange for latest version
http://www.mathworks.com/matlabcentral/fileexchange/3026
"""
n = len(array_ext)
po = 0
# The original rainflow counting by Nieslony, unchanged
j = -1
c_nr1 = 1
for i in range(n):
j += 1
a[j] = array_ext[i]
t[j] = array_t[i]
while (j >= 2) and (abs(a[j - 1] - a[j - 2]) <= abs(a[j] - a[j - 1])):
ampl = abs((a[j - 1] - a[j - 2]) / 2)
mean = (a[j - 1] + a[j - 2]) / 2
period = (t[j - 1] - t[j - 2]) * 2
atime = t[j - 2]
if j == 2:
a[0] = a[1]
a[1] = a[2]
t[0] = t[1]
t[1] = t[2]
j = 1
if (ampl > 0):
array_out[po, :] = (ampl, mean, 0.5, atime, period)
po += 1
else:
a[j - 2] = a[j]
t[j - 2] = t[j]
j = j - 2
if (ampl > 0):
array_out[po, :] = (ampl, mean, 1.0, atime, period)
po += 1
c_nr1 += 1
c_nr2 = 1
for i in range(j):
# for (i=0; i<j; i++) {
ampl = abs(a[i] - a[i + 1]) / 2
mean = (a[i] + a[i + 1]) / 2
period = (t[i + 1] - t[i]) * 2
atime = t[i]
if (ampl > 0):
array_out[po, :] = (ampl, mean, 0.5, atime, period)
po += 1
c_nr2 += 1
return c_nr1 #, c_nr2
def findrfc_astm(tp, t=None):
"""
Return rainflow counted cycles
Nieslony's Matlab implementation of the ASTM standard practice for rainflow
counting ported to a Python C module.
Parameters
----------
tp : array-like
vector of turningpoints (NB! Only values, not sampled times)
t : array-like
vector of sampled times
Returns
-------
sig_rfc : array-like
array of shape (n,3) or (n, 5) with:
sig_rfc[:,0] Cycles amplitude
sig_rfc[:,1] Cycles mean value
sig_rfc[:,2] Cycle type, half (=0.5) or full (=1.0)
sig_rfc[:,3] cycle_begin_time (only if t is given)
sig_rfc[:,4] cycle_period_time (only if t is given)
"""
y1 = np.atleast_1d(tp).ravel()
n = len(y1)
a = np.zeros(n)
if t is None:
sig_rfc = np.zeros((n, 3))
cnr = _findrfc3_astm(y1, a, sig_rfc)
else:
t1 = np.atleast_1d(t).ravel()
sig_rfc = np.zeros((n, 5))
t2 = np.zeros(n)
cnr = _findrfc5_astm(y1, t1, a, t2, sig_rfc)
# the sig_rfc was constructed too big in rainflow.rf3, so
# reduce the sig_rfc array as done originally by a matlab mex c function
# n = len(sig_rfc)
return sig_rfc[:n - cnr]
if __name__ == '__main__':
pass