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

Nieslony's ASTM rainflow counting algorithm. Partially integrated, no support for the CyclePairs object yet.

Browse Source
master
david.verelst 14 years ago
parent 79fa7ab190
commit c777ef2a43
5 changed files with 310 additions and 20 deletions
Show Stats Download Patch File Download Diff File

BIN
pywafo/src/wafo/c_library.so
View File

Binary file not shown.

78
pywafo/src/wafo/misc.py
Unescape Escape View File

@ -477,6 +477,36 @@ def findpeaks(data, n=2, min_h=None, min_p=0.0):
return ind
def findrfc_astm(tp):
"""
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)
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)
"""
y1 = atleast_1d(tp).ravel()
sig_rfc, cnr = clib.findrfc3_astm(y1)
# 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)
sig_rfc = sig_rfc.__getslice__(0,n-cnr[0])
# sig_rfc holds the actual rainflow counted cycles, not the indices
return sig_rfc
def findrfc(tp, hmin=0.0, method='clib'):
'''
@ -489,6 +519,10 @@ def findrfc(tp, hmin=0.0, method='clib'):
h : real scalar
rainflow threshold. If h>0, then all rainflow cycles with height
smaller than h are removed.
method : string, optional
'clib' 'None'
Specify 'clib' for calling the c_functions, otherwise fallback to
the Python implementation.
Returns
-------
@ -513,8 +547,9 @@ def findrfc(tp, hmin=0.0, method='clib'):
rfcfilter,
findtp.
'''
# TODO merge rfcfilter and findrfc
# TODO: merge rfcfilter and findrfc
y1 = atleast_1d(tp).ravel()
n = len(y1)
ind = zeros(0, dtype=np.int)
ix = 0
@ -539,7 +574,7 @@ def findrfc(tp, hmin=0.0, method='clib'):
warnings.warn('This is not a sequence of turningpoints, exit')
return ind
if clib is None or method!='clib':
if clib is None or method not in ('clib'):
ind = zeros(n, dtype=np.int)
NC = np.int(NC)
for i in xrange(NC):
@ -717,8 +752,8 @@ def findtp(x, h=0.0, kind=None):
if h>0, then all rainflow cycles with height smaller than
h are removed.
kind : string
defines the type of wave. Possible options are
'mw' 'Mw' or 'none'.
defines the type of wave or indicate the ASTM rainflow counting method.
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
@ -765,23 +800,32 @@ def findtp(x, h=0.0, kind=None):
if ind.size < 2:
return None
#% In order to get the exact up-crossing intensity from rfc by
#% mm2lc(tp2mm(rfc)) we have to add the indices
#% to the last value (and also the first if the
#% sequence of turning points does not start with a minimum).
if x[ind[0]] > x[ind[1]]:
#% adds indices to first and last value
ind = r_[0, ind, n - 1]
else: # adds index to the last value
ind = r_[ind, n - 1]
if kind == 'astm':
# the Nieslony approach always put the first loading point as the first
# turning point.
if x[ind[0]] != x[0]:
# add the first turning point is the first of the signal
ind = np.r_[0, ind, n - 1]
else:
# only add the last point of the signal
ind = np.r_[ind, n - 1]
else:
if x[ind[0]] > x[ind[1]]:
#% adds indices to first and last value
ind = r_[0, ind, n - 1]
else: # adds index to the last value
ind = r_[ind, n - 1]
if h > 0.0:
ind1 = findrfc(x[ind], h)
ind = ind[ind1]
if kind in ('mw', 'Mw'):
xor = lambda a, b: a ^ b
# make sure that the first is a Max if wdef == 'Mw'
@ -791,9 +835,9 @@ def findtp(x, h=0.0, kind=None):
remove_first = xor(first_is_max, kind.startswith('Mw'))
if remove_first:
ind = ind[1::]
# make sure the number of minima and Maxima are according to the wavedef.
# i.e., make sure Nm=length(ind) is odd
# make sure the number of minima and Maxima are according to the
# wavedef. i.e., make sure Nm=length(ind) is odd
if (mod(ind.size, 2)) != 1:
ind = ind[:-1]
return ind

47
pywafo/src/wafo/objects.py
Unescape Escape View File

@ -17,7 +17,7 @@ from wafo.transform.core import TrData
from wafo.transform.models import TrHermite, TrOchi, TrLinear
from wafo.stats import edf, distributions
from wafo.misc import (nextpow2, findtp, findrfc, findtc, findcross,
ecross, JITImport, DotDict, gravity)
ecross, JITImport, DotDict, gravity, findrfc_astm)
from wafodata import WafoData
from wafo.interpolate import SmoothSpline
from scipy.interpolate.interpolate import interp1d
@ -909,6 +909,10 @@ class TurningPoints(WafoData):
----------
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
------
@ -930,6 +934,7 @@ class TurningPoints(WafoData):
TurningPoints
SurvivalCycleCount
"""
if h>0:
ind = findrfc(self.data, h, method=method)
data = self.data[ind]
@ -951,8 +956,41 @@ class TurningPoints(WafoData):
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)
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)
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):
"""
@ -1543,11 +1581,14 @@ class TimeSeries(WafoData):
wavetype : string
defines the type of wave. Possible options are
'mw' 'Mw' or 'none'.
'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
-------

164
pywafo/src/wafo/source/c_codes/c_functions.c
Unescape Escape View File

@ -613,3 +613,167 @@ void disufq2(double *rsvec, double *isvec,
// return i;
}
/* ++++++++++ BEGIN RF3 [ampl ampl_mean nr_of_cycle] */
/* ++++++++++ Rain flow without time analysis */
//By Adam Nieslony
//Visit the MATLAB Central File Exchange for latest version
//http://www.mathworks.com/matlabcentral/fileexchange/3026
void findrfc3_astm(double *array_ext, double *array_out, int n, int *nout) {
double *pr, *po, a[16384], ampl, mean;
int tot_num, index, j, cNr1, cNr2;
tot_num = n;
// pointers to the first element of the arrays
pr = &array_ext[0];
po = &array_out[0];
// The original rainflow counting by Nieslony, unchanged
j = -1;
cNr1 = 1;
for (index=0; index<tot_num; index++) {
a[++j]=*pr++;
while ( (j >= 2) && (fabs(a[j-1]-a[j-2]) <= fabs(a[j]-a[j-1])) ) {
ampl=fabs( (a[j-1]-a[j-2])/2 );
switch(j) {
case 0: { break; }
case 1: { break; }
case 2: {
mean=(a[0]+a[1])/2;
a[0]=a[1];
a[1]=a[2];
j=1;
if (ampl > 0) {
*po++=ampl;
*po++=mean;
*po++=0.50;
}
break;
}
default: {
mean=(a[j-1]+a[j-2])/2;
a[j-2]=a[j];
j=j-2;
if (ampl > 0) {
*po++=ampl;
*po++=mean;
*po++=1.00;
cNr1++;
}
break;
}
}
}
}
cNr2 = 1;
for (index=0; index<j; index++) {
ampl=fabs(a[index]-a[index+1])/2;
mean=(a[index]+a[index+1])/2;
if (ampl > 0){
*po++=ampl;
*po++=mean;
*po++=0.50;
cNr2++;
}
}
// array of ints nout is outputted
nout[0] = cNr1;
nout[1] = cNr2;
}
/* ++++++++++ END RF3 */
// ++ BEGIN RF5 [ampl ampl_mean nr_of_cycle cycle_begin_time cycle_period_time]
/* ++++++++++ Rain flow with time analysis */
//By Adam Nieslony
//Visit the MATLAB Central File Exchange for latest version
//http://www.mathworks.com/matlabcentral/fileexchange/3026
void
findrfc5_astm(double *array_ext, double *array_t, double *array_out, int n, int *nout) {
double *pr, *pt, *po, a[16384], t[16384], ampl, mean, period, atime;
int tot_num, index, j, cNr1, cNr2;
// tot_num = mxGetM(array_ext) * mxGetN(array_ext);
tot_num = n;
// pointers to the first element of the arrays
pr = &array_ext[0];
pt = &array_t[0];
po = &array_out[0];
// array_out = mxCreateDoubleMatrix(5, tot_num-1, mxREAL);
// The original rainflow counting by Nieslony, unchanged
j = -1;
cNr1 = 1;
for (index=0; index<tot_num; index++) {
a[++j]=*pr++;
t[j]=*pt++;
while ( (j >= 2) && (fabs(a[j-1]-a[j-2]) <= fabs(a[j]-a[j-1])) ) {
ampl=fabs( (a[j-1]-a[j-2])/2 );
switch(j)
{
case 0: { break; }
case 1: { break; }
case 2: {
mean=(a[0]+a[1])/2;
period=(t[1]-t[0])*2;
atime=t[0];
a[0]=a[1];
a[1]=a[2];
t[0]=t[1];
t[1]=t[2];
j=1;
if (ampl > 0) {
*po++=ampl;
*po++=mean;
*po++=0.50;
*po++=atime;
*po++=period;
}
break;
}
default: {
mean=(a[j-1]+a[j-2])/2;
period=(t[j-1]-t[j-2])*2;
atime=t[j-2];
a[j-2]=a[j];
t[j-2]=t[j];
j=j-2;
if (ampl > 0) {
*po++=ampl;
*po++=mean;
*po++=1.00;
*po++=atime;
*po++=period;
cNr1++;
}
break;
}
}
}
}
cNr2 = 1;
for (index=0; index<j; index++) {
ampl=fabs(a[index]-a[index+1])/2;
mean=(a[index]+a[index+1])/2;
period=(t[index+1]-t[index])*2;
atime=t[index];
if (ampl > 0){
*po++=ampl;
*po++=mean;
*po++=0.50;
*po++=atime;
*po++=period;
cNr2++;
}
}
// /* free the memeory !!!*/
// mxSetN(array_out, tot_num - cNr);
nout[0] = cNr1;
nout[1] = cNr2;
}
/* ++++++++++ END RF5 */

41
pywafo/src/wafo/source/c_codes/c_library.pyf
Unescape Escape View File

@ -41,5 +41,46 @@ interface
integer intent(in) :: nmin, nmax
double precision dimension(n*m), intent(out) :: rsvec, isvec, rdvec, idvec ! output array,
end subroutine disufq2
! ===== START NIESLONY RAINFLOW FUNCTIONS
! RAINFLOW Revision: 1.1
! by Adam Nieslony, 2009
subroutine findrfc3_astm(array_ext, array_out, n, nout)
intent(c) findrfc3_astm ! rf3 is a C function
intent(c) ! all rf3 arguments are
! considered as C based
! n is the length of the input array array_ext
integer intent(hide), depend(array_ext) :: n=len(array_ext)
! of input array x
double precision intent(in) :: array_ext(n)
! the output array
double precision intent(out) :: array_out(n,3)
! nout array, to output additional ints
integer dimension(2), intent(out) :: nout
end subroutine findrfc3_astm
subroutine findrfc5_astm(array_ext, array_t, array_out, n, nout)
intent(c) findrfc5_astm ! rf5 is a C function
intent(c) ! all rf5 arguments are
! considered as C based
! n is the length of the input array array_ext
integer intent(hide), depend(array_ext) :: n=len(array_ext)
! of input array x
double precision intent(in) :: array_ext(n), array_t(n)
! the output array
double precision intent(out) :: array_out(n,5)
! nout array, to output additional ints
integer dimension(2), intent(out) :: nout
end subroutine findrfc5_astm
! ===== END NIESLONY RAINFLOW FUNCTIONS
end interface
end python module c_library
Write Preview
Loading…
Cancel
Save