diff --git a/pywafo/src/wafo/c_library.so b/pywafo/src/wafo/c_library.so index 3ddfc76..f382a86 100755 Binary files a/pywafo/src/wafo/c_library.so and b/pywafo/src/wafo/c_library.so differ diff --git a/pywafo/src/wafo/misc.py b/pywafo/src/wafo/misc.py index 01ddfd2..aeb5694 100644 --- a/pywafo/src/wafo/misc.py +++ b/pywafo/src/wafo/misc.py @@ -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 diff --git a/pywafo/src/wafo/objects.py b/pywafo/src/wafo/objects.py index d27f574..6d641a8 100644 --- a/pywafo/src/wafo/objects.py +++ b/pywafo/src/wafo/objects.py @@ -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 ------- diff --git a/pywafo/src/wafo/source/c_codes/c_functions.c b/pywafo/src/wafo/source/c_codes/c_functions.c index 7725902..a9bbdd5 100644 --- a/pywafo/src/wafo/source/c_codes/c_functions.c +++ b/pywafo/src/wafo/source/c_codes/c_functions.c @@ -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= 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 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= 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 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 */ diff --git a/pywafo/src/wafo/source/c_codes/c_library.pyf b/pywafo/src/wafo/source/c_codes/c_library.pyf index 694ec54..5e58404 100644 --- a/pywafo/src/wafo/source/c_codes/c_library.pyf +++ b/pywafo/src/wafo/source/c_codes/c_library.pyf @@ -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 \ No newline at end of file