diff --git a/README.md b/README.md
index 0be201ee70b9c7e06f187461093b1dd72d92ed48..7d6778965c5e54d703e157525fd9133134184e00 100644
--- a/README.md
+++ b/README.md
@@ -10,13 +10,14 @@
- [at_time_file](wetb/hawc2/at_time_file.py): read at output_at_time files
- [ascii2bin](wetb/hawc2/ascii2bin): Compress HAWC2 ascii result files to binary
-
-
### [gtsdf](wetb/gtsdf)
General Time Series Data Format, a binary hdf5 data format for storing time series data.
- [gtsdf](wetb/gtsdf/gtsdf.py): read/write/append gtsdf files
- [unix_time](wetb/gtsdf/unix_time.py): convert between datetime and unix time (seconds since 1/1/1970)
+### [fatigue_tools](wetb/fatigue_tools)
+- [fatigue](wetb/fatigue_tools/fatigue.py): Rainflow counting, cycle matrix and equvivalent loads
+
### [FAST](wetb/fast)
Tools for working with NREL's FAST code (An aeroelastic computer-aided engineering (CAE) tool for horizontal axis wind turbines)
- [fast_io](wetb/fast/fast_io.py): Read binary and ascii result files
diff --git a/wetb/fatigue_tools/__init__.py b/wetb/fatigue_tools/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/wetb/fatigue_tools/fatigue.py b/wetb/fatigue_tools/fatigue.py
new file mode 100644
index 0000000000000000000000000000000000000000..e0b57e70deb6d692018ebd74d4b5af7a346758c6
--- /dev/null
+++ b/wetb/fatigue_tools/fatigue.py
@@ -0,0 +1,350 @@
+'''
+Created on 04/03/2013
+@author: mmpe
+'''
+
+#try:
+# """
+# The cython_import function compiles modules using cython.
+# It is found at: https://github.com/madsmpedersen/MMPE/blob/master/cython_compile/cython_compile.py
+# """
+# from mmpe.cython_compile.cython_compile import cython_import
+#except ImportError:
+# cython_import = __import__
+import numpy as np
+
+
+def rfc_hist(sig_rf, nrbins=46):
+ """
+ Histogram of rainflow counted cycles
+ ====================================
+
+ hist, bin_edges, bin_avg = rfc_hist(sig, nrbins=46)
+
+ Divide the rainflow counted cycles of a signal into equally spaced bins.
+
+ Created on Wed Feb 16 16:53:18 2011
+ @author: David Verelst
+ Modified 10.10.2011 by Mads M Pedersen to elimintate __copy__ and __eq__
+
+ Parameters
+ ----------
+ sig_rf : array-like
+ As output by rfc_astm or rainflow
+
+ nrbins : int, optional
+ Divide the rainflow counted amplitudes in a number of equally spaced
+ bins.
+
+ Returns
+ -------
+ hist : array-like
+ Counted rainflow cycles per bin, has nrbins elements
+
+ bin_edges : array-like
+ Edges of the bins, has nrbins+1 elements.
+
+ bin_avg : array-like
+ Average rainflow cycle amplitude per bin, has nrbins elements.
+ """
+
+ rf_half = sig_rf
+
+ # the Matlab approach is to divide into 46 bins
+ bin_edges = np.linspace(0, 1, num=nrbins + 1) * rf_half.max()
+ hist = np.histogram(rf_half, bins=bin_edges)[0]
+ # calculate the average per bin
+ hist_sum = np.histogram(rf_half, weights=rf_half, bins=bin_edges)[0]
+ # replace zeros with one, to avoid 0/0
+ hist_ = hist.copy()
+ hist_[(hist == 0).nonzero()] = 1.0
+ # since the sum is also 0, the avg remains zero for those whos hist is zero
+ bin_avg = hist_sum / hist_
+
+ return hist, bin_edges, bin_avg
+
+
+def check_signal(signal):
+ # check input data validity
+ if not type(signal).__name__ == 'ndarray':
+ raise TypeError('signal must be ndarray, not: ' + type(signal).__name__)
+
+ elif len(signal.shape) not in (1, 2):
+ raise TypeError('signal must be 1D or 2D, not: ' + str(len(signal.shape)))
+
+ if len(signal.shape) == 2:
+ if signal.shape[1] > 1:
+ raise TypeError('signal must have one column only, not: ' + str(signal.shape[1]))
+ if np.min(signal) == np.max(signal):
+ raise TypeError("Signal contains no variation")
+
+
+def rainflow_windap(signal, levels=255., thresshold=(255 / 50)):
+ """Windap equivalent rainflow counting
+
+
+ Calculate the amplitude and mean values of half cycles in signal
+
+ This algorithms used by this routine is implemented directly as described in
+ "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A
+
+ Parameters
+ ----------
+ Signal : array-like
+ The raw signal
+
+ levels : int, optional
+ The signal is discretize into this number of levels.
+ 255 is equivalent to the implementation in Windap
+
+ thresshold : int, optional
+ Cycles smaller than this thresshold are ignored
+ 255/50 is equivalent to the implementation in Windap
+
+ Returns
+ -------
+ ampl : array-like
+ Peak to peak amplitudes of the half cycles
+
+ mean : array-like
+ Mean values of the half cycles
+
+
+ Example
+ -------
+ >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+ >>> ampl, mean = rainflow_windap(signal)
+ """
+ check_signal(signal)
+ #type <double> is required by <find_extreme> and <rainflow>
+ signal = signal.astype(np.double)
+
+ offset = np.nanmin(signal)
+ signal -= offset
+ if np.nanmax(signal) > 0:
+ gain = np.nanmax(signal) / levels
+ signal = signal / gain
+ signal = np.round(signal).astype(np.int)
+
+ from wetb.fatigue_tools.peak_trough import peak_trough
+
+
+ #Convert to list of local minima/maxima where difference > thresshold
+ sig_ext = peak_trough(signal, thresshold)
+
+ from wetb.fatigue_tools.pair_range import pair_range_amplitude_mean
+
+ #rainflow count
+ ampl_mean = pair_range_amplitude_mean(sig_ext)
+
+ ampl_mean = np.array(ampl_mean)
+ ampl_mean = np.round(ampl_mean / thresshold) * gain * thresshold
+ ampl_mean[:, 1] += offset
+ return ampl_mean.T
+
+
+
+def rainflow_astm(signal):
+ """Matlab equivalent rainflow counting
+
+ Calculate the amplitude and mean values of half cycles in signal
+
+ This implemementation is based on the c-implementation by Adam Nieslony found at
+ the MATLAB Central File Exchange http://www.mathworks.com/matlabcentral/fileexchange/3026
+
+ Parameters
+ ----------
+ Signal : array-like
+ The raw signal
+
+ Returns
+ -------
+ ampl : array-like
+ peak to peak amplitudes of the half cycles (note that the matlab implementation
+ uses peak amplitude instead of peak to peak)
+
+ mean : array-like
+ Mean values of the half cycles
+
+
+ Examples
+ --------
+ >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+ >>> ampl, mean = rainflow_astm(signal)
+ """
+ check_signal(signal)
+
+ # type <double> is reuqired by <find_extreme> and <rainflow>
+ signal = signal.astype(np.double)
+
+ # Import find extremes and rainflow.
+ # If possible the module is compiled using cython otherwise the python implementation is used
+ #cython_import('fatigue_tools.rainflowcount_astm')
+ from wetb.fatigue_tools.rainflowcount_astm import find_extremes, rainflowcount
+
+ # Remove points which is not local minimum/maximum
+ sig_ext = find_extremes(signal)
+
+ # rainflow count
+ ampl_mean = np.array(rainflowcount(sig_ext))
+
+ return np.array(ampl_mean).T
+
+
+def eq_load(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1, rainflow_func=rainflow_windap):
+ """Equivalent load calculation
+
+ Calculate the equivalent loads for a list of Wohler exponent and number of equivalent loads
+
+ Parameters
+ ----------
+ signals : list of tuples or array_like
+ - if list of tuples: list must have format [(sig1_weight, sig1),(sig2_weight, sig1),...] where\n
+ - sigx_weight is the weight of signal x\n
+ - sigx is signal x\n
+ - if array_like: The signal
+ no_bins : int, optional
+ Number of bins in rainflow count histogram
+ m : int, float or array-like, optional
+ Wohler exponent (default is [3, 4, 6, 8, 10, 12])
+ neq : int, float or array-like, optional
+ The equivalent number of load cycles (default is 1, but normally the time duration in seconds is used)
+ rainflow_func : {rainflow_windap, rainflow_astm}, optional
+ The rainflow counting function to use (default is rainflow_windap)
+
+ Returns
+ -------
+ eq_loads : array-like
+ List of lists of equivalent loads for the corresponding equivalent number(s) and Wohler exponents
+
+ Examples
+ --------
+ >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+ >>> eq_load(signal, no_bins=50, neq=[1, 17], m=[3, 4, 6], rainflow_func=rainflow_windap)
+ [[10.311095426959747, 9.5942535021382174, 9.0789213365013932], # neq = 1, m=[3,4,6]
+ [4.010099657859783, 4.7249689509841746, 5.6618639965313005]], # neq = 17, m=[3,4,6]
+
+ eq_load([(.4, signal), (.6, signal)], no_bins=50, neq=[1, 17], m=[3, 4, 6], rainflow_func=rainflow_windap)
+ [[10.311095426959747, 9.5942535021382174, 9.0789213365013932], # neq = 1, m=[3,4,6]
+ [4.010099657859783, 4.7249689509841746, 5.6618639965313005]], # neq = 17, m=[3,4,6]
+ """
+ try:
+ return eq_load_and_cycles(signals, no_bins, m, neq, rainflow_func)[0]
+ except TypeError:
+ return [[np.nan] * len(np.atleast_1d(m))] * len(np.atleast_1d(neq))
+
+# ampl, _ = rainflow_func(signals)
+# if ampl is None:
+# return []
+# hist_data, x, bin_avg = rfc_hist(ampl, no_bins)
+#
+# m = np.atleast_1d(m)
+#
+# return np.array([np.power(np.sum(0.5 * hist_data * np.power(bin_avg, m[i])) / neq, 1. / m[i]) for i in range(len(m))])
+
+
+def eq_load_and_cycles(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=[10 ** 6, 10 ** 7, 10 ** 8], rainflow_func=rainflow_windap):
+ """Calculate combined fatigue equivalent load
+
+ Parameters
+ ----------
+ signals : list of tuples or array_like
+ - if list of tuples: list must have format [(sig1_weight, sig1),(sig2_weight, sig1),...] where\n
+ - sigx_weight is the weight of signal x\n
+ - sigx is signal x\n
+ - if array_like: The signal
+ no_bins : int, optional
+ Number of bins for rainflow counting
+ m : int, float or array-like, optional
+ Wohler exponent (default is [3, 4, 6, 8, 10, 12])
+ neq : int or array-like, optional
+ Equivalent number, default is [10^6, 10^7, 10^8]
+ rainflow_func : {rainflow_windap, rainflow_astm}, optional
+ The rainflow counting function to use (default is rainflow_windap)
+
+ Returns
+ -------
+ eq_loads : array-like
+ List of lists of equivalent loads for the corresponding equivalent number(s) and Wohler exponents
+ cycles : array_like
+ 2d array with shape = (no_ampl_bins, no_mean_bins)
+ ampl_bin_mean : array_like
+ mean amplitude of the bins
+ ampl_bin_edges
+ Edges of the amplitude bins
+ """
+ cycles, ampl_bin_mean, ampl_bin_edges, _, _ = cycle_matrix(signals, no_bins, 1, rainflow_func)
+ if 0: #to be similar to windap
+ ampl_bin_mean = (ampl_bin_edges[:-1] + ampl_bin_edges[1:]) / 2
+ cycles, ampl_bin_mean = cycles.flatten(), ampl_bin_mean.flatten()
+ eq_loads = [[((np.sum(cycles * ampl_bin_mean ** _m) / _neq) ** (1. / _m)) for _m in np.atleast_1d(m)] for _neq in np.atleast_1d(neq)]
+ return eq_loads, cycles, ampl_bin_mean, ampl_bin_edges
+
+
+
+def cycle_matrix(signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_windap):
+ """Markow load cycle matrix
+
+ Calculate the Markow load cycle matrix
+
+ Parameters
+ ----------
+ Signals : array-like or list of tuples
+ if array-like, the raw signal
+ if list of tuples, list of (weight, signal), e.g. [(0.1,sig1), (0.8,sig2), (.1,sig3)]
+ ampl_bins : int or array-like, optional
+ if int, Number of amplitude value bins (default is 10)
+ if array-like, the bin edges for amplitude
+ mean_bins : int or array-like, optional
+ if int, Number of mean value bins (default is 10)
+ if array-like, the bin edges for mea
+ rainflow_func : {rainflow_windap, rainflow_astm}, optional
+ The rainflow counting function to use (default is rainflow_windap)
+
+ Returns
+ -------
+ cycles : ndarray, shape(ampl_bins, mean_bins)
+ A bi-dimensional histogram of load cycles(full cycles). Amplitudes are histogrammed along the first dimension and mean values are histogrammed along the second dimension.
+ ampl_bin_mean : ndarray, shape(ampl_bins,)
+ The average cycle amplitude of the bins
+ ampl_edges : ndarray, shape(ampl_bins+1,)
+ The amplitude bin edges
+ mean_bin_mean : ndarray, shape(ampl_bins,)
+ The average cycle mean of the bins
+ mean_edges : ndarray, shape(mean_bins+1,)
+ The mean bin edges
+
+ Examples
+ --------
+ >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+ >>> cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges = cycle_matrix(signal)
+ >>> cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges = cycle_matrix([(.4, signal), (.6,signal)])
+ """
+
+ if isinstance(signals[0], tuple):
+ ampls = np.empty((0,), dtype=np.float64)
+ means = np.empty((0,), dtype=np.float64)
+ weights = np.empty((0,), dtype=np.float64)
+ for w, signal in signals:
+ a, m = rainflow_func(signal[:])
+ ampls = np.r_[ampls, a]
+ means = np.r_[means, m]
+ weights = np.r_[weights, (np.zeros_like(a) + w)]
+ else:
+ ampls, means = rainflow_func(signals[:])
+ weights = np.ones_like(ampls)
+ if isinstance(ampl_bins, int):
+ ampl_bins = np.linspace(0, 1, num=ampl_bins + 1) * ampls.max()
+ cycles, ampl_edges, mean_edges = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights)
+
+ ampl_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * ampls)[0]
+ ampl_bin_mean = np.zeros_like(cycles)
+ mask = (cycles > 0)
+ ampl_bin_mean[mask] = ampl_bin_sum[mask] / cycles[mask]
+ mean_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * means)[0]
+ mean_bin_mean = np.zeros_like(cycles)
+ mean_bin_mean[cycles > 0] = mean_bin_sum[cycles > 0] / cycles[cycles > 0]
+ cycles = cycles / 2 # to get full cycles
+ return cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges
+
+
diff --git a/wetb/fatigue_tools/pair_range.py b/wetb/fatigue_tools/pair_range.py
new file mode 100644
index 0000000000000000000000000000000000000000..6a493d701e5760682061369f6f2c1414f8d63ce1
--- /dev/null
+++ b/wetb/fatigue_tools/pair_range.py
@@ -0,0 +1,198 @@
+
+import cython
+import numpy as np
+# cimport numpy as np
+
+
+@cython.locals(p=cython.int, q=cython.int, f=cython.int, flow=list, k=cython.int, n=cython.int, ptr=cython.int)
+def pair_range_amplitude(x): # cpdef pair_range(np.ndarray[long,ndim=1] x):
+ """
+ Returns a list of half-cycle-amplitudes
+ x: Peak-Trough sequence (integer list of local minima and maxima)
+
+ This routine is implemented according to
+ "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A
+ except that a list of half-cycle-amplitudes are returned instead of a from_level-to_level-matrix
+ """
+
+ x = x - np.min(x)
+ k = np.max(x)
+ n = x.shape[0]
+ S = np.zeros(n + 1)
+
+ #A = np.zeros(k+1)
+ flow = []
+ S[1] = x[0]
+ ptr = 1
+ p = 1
+ q = 1
+ f = 0
+ # phase 1
+ while True:
+ p += 1
+ q += 1
+
+ # read
+ S[p] = x[ptr]
+ ptr += 1
+
+ if q == n:
+ f = 1
+ while p >= 4:
+ if (S[p - 2] > S[p - 3] and S[p - 1] >= S[p - 3] and S[p] >= S[p - 2]) \
+ or\
+ (S[p - 2] < S[p - 3] and S[p - 1] <= S[p - 3] and S[p] <= S[p - 2]):
+ ampl = abs(S[p - 2] - S[p - 1])
+ # A[ampl]+=2 #Two half cycles
+ flow.append(ampl)
+ flow.append(ampl)
+ S[p - 2] = S[p]
+
+ p -= 2
+ else:
+ break
+
+ if f == 0:
+ pass
+ else:
+ break
+ # phase 2
+ q = 0
+ while True:
+ q += 1
+ if p == q:
+ break
+ else:
+ ampl = abs(S[q + 1] - S[q])
+ # A[ampl]+=1
+ flow.append(ampl)
+ return flow
+
+
+
+
+
+@cython.locals(p=cython.int, q=cython.int, f=cython.int, flow=list, k=cython.int, n=cython.int, ptr=cython.int)
+def pair_range_from_to(x): # cpdef pair_range(np.ndarray[long,ndim=1] x):
+ """
+ Returns a list of half-cycle-amplitudes
+ x: Peak-Trough sequence (integer list of local minima and maxima)
+
+ This routine is implemented according to
+ "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A
+ except that a list of half-cycle-amplitudes are returned instead of a from_level-to_level-matrix
+ """
+
+ x = x - np.min(x)
+ k = np.max(x)
+ n = x.shape[0]
+ S = np.zeros(n + 1)
+
+ A = np.zeros((k + 1, k + 1))
+ S[1] = x[0]
+ ptr = 1
+ p = 1
+ q = 1
+ f = 0
+ # phase 1
+ while True:
+ p += 1
+ q += 1
+
+ # read
+ S[p] = x[ptr]
+ ptr += 1
+
+ if q == n:
+ f = 1
+ while p >= 4:
+ #print S[p - 3:p + 1]
+ #print S[p - 2], ">", S[p - 3], ", ", S[p - 1], ">=", S[p - 3], ", ", S[p], ">=", S[p - 2], (S[p - 2] > S[p - 3] and S[p - 1] >= S[p - 3] and S[p] >= S[p - 2])
+ #print S[p - 2], "<", S[p - 3], ", ", S[p - 1], "<=", S[p - 3], ", ", S[p], "<=", S[p - 2], (S[p - 2] < S[p - 3] and S[p - 1] <= S[p - 3] and S[p] <= S[p - 2])
+ #print (S[p - 2] > S[p - 3] and S[p - 1] >= S[p - 3] and S[p] >= S[p - 2]) or (S[p - 2] < S[p - 3] and S[p - 1] <= S[p - 3] and S[p] <= S[p - 2])
+ if (S[p - 2] > S[p - 3] and S[p - 1] >= S[p - 3] and S[p] >= S[p - 2]) or \
+ (S[p - 2] < S[p - 3] and S[p - 1] <= S[p - 3] and S[p] <= S[p - 2]):
+ A[S[p - 2], S[p - 1]] += 1
+ A[S[p - 1], S[p - 2]] += 1
+ S[p - 2] = S[p]
+ p -= 2
+ else:
+ break
+
+ if f == 1:
+ break # q==n
+ # phase 2
+ q = 0
+ while True:
+ q += 1
+ if p == q:
+ break
+ else:
+ #print S[q], "to", S[q + 1]
+ A[S[q], S[q + 1]] += 1
+ return A
+
+@cython.locals(p=cython.int, q=cython.int, f=cython.int, flow=list, k=cython.int, n=cython.int, ptr=cython.int)
+def pair_range_amplitude_mean(x): # cpdef pair_range(np.ndarray[long,ndim=1] x):
+ """
+ Returns a list of half-cycle-amplitudes
+ x: Peak-Trough sequence (integer list of local minima and maxima)
+
+ This routine is implemented according to
+ "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A
+ except that a list of half-cycle-amplitudes are returned instead of a from_level-to_level-matrix
+ """
+
+ x = x - np.min(x)
+ k = np.max(x)
+ n = x.shape[0]
+ S = np.zeros(n + 1)
+ ampl_mean = []
+ A = np.zeros((k + 1, k + 1))
+ S[1] = x[0]
+ ptr = 1
+ p = 1
+ q = 1
+ f = 0
+ # phase 1
+ while True:
+ p += 1
+ q += 1
+
+ # read
+ S[p] = x[ptr]
+ ptr += 1
+
+ if q == n:
+ f = 1
+ while p >= 4:
+ if (S[p - 2] > S[p - 3] and S[p - 1] >= S[p - 3] and S[p] >= S[p - 2]) \
+ or\
+ (S[p - 2] < S[p - 3] and S[p - 1] <= S[p - 3] and S[p] <= S[p - 2]):
+ # Extract two intermediate half cycles
+ ampl = abs(S[p - 2] - S[p - 1])
+ mean = (S[p - 2] + S[p - 1]) / 2
+ ampl_mean.append((ampl, mean))
+ ampl_mean.append((ampl, mean))
+
+ S[p - 2] = S[p]
+
+ p -= 2
+ else:
+ break
+
+ if f == 0:
+ pass
+ else:
+ break
+ # phase 2
+ q = 0
+ while True:
+ q += 1
+ if p == q:
+ break
+ else:
+ ampl = abs(S[q + 1] - S[q])
+ mean = (S[q + 1] + S[q]) / 2
+ ampl_mean.append((ampl, mean))
+ return ampl_mean
diff --git a/wetb/fatigue_tools/pair_range.pyd b/wetb/fatigue_tools/pair_range.pyd
new file mode 100644
index 0000000000000000000000000000000000000000..e3e28cb4e3bb7ab0fc866c3ab7d1ec242458abdb
Binary files /dev/null and b/wetb/fatigue_tools/pair_range.pyd differ
diff --git a/wetb/fatigue_tools/peak_trough.py b/wetb/fatigue_tools/peak_trough.py
new file mode 100644
index 0000000000000000000000000000000000000000..fc6bde699eee1cbd1406f58bca043a5f637bfa42
--- /dev/null
+++ b/wetb/fatigue_tools/peak_trough.py
@@ -0,0 +1,100 @@
+import cython
+import numpy as np
+#cimport numpy as np
+
+@cython.locals(BEGIN=cython.int, MINZO=cython.int, MAXZO=cython.int, ENDZO=cython.int, \
+ R=cython.int, L=cython.int, i=cython.int, p=cython.int, f=cython.int)
+def peak_trough(x, R): #cpdef np.ndarray[long,ndim=1] peak_trough(np.ndarray[long,ndim=1] x, int R):
+ """
+ Returns list of local maxima/minima.
+
+ x: 1-dimensional numpy array containing signal
+ R: Thresshold (minimum difference between succeeding min and max
+
+ This routine is implemented directly as described in
+ "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A
+ """
+
+ BEGIN = 0
+ MINZO = 1
+ MAXZO = 2
+ ENDZO = 3
+ S = np.zeros(x.shape[0] + 1, dtype=np.int)
+
+ L = x.shape[0]
+ goto = BEGIN
+
+ while 1:
+ if goto == BEGIN:
+ trough = x[0]
+ peak = x[0]
+
+ i = 0
+ p = 1
+ f = 0
+ while goto == BEGIN:
+ i += 1
+ if i == L:
+ goto = ENDZO
+ continue
+ else:
+ if x[i] > peak:
+ peak = x[i]
+ if peak - trough >= R:
+ S[p] = trough
+ goto = MAXZO
+ continue
+ elif x[i] < trough:
+ trough = x[i]
+ if peak - trough >= R:
+ S[p] = peak
+ goto = MINZO
+ continue
+
+ elif goto == MINZO:
+ f = -1
+
+ while goto == MINZO:
+ i += 1
+ if i == L:
+ goto = ENDZO
+ continue
+ else:
+ if x[i] < trough:
+ trough = x[i]
+ else:
+ if x[i] - trough >= R:
+ p += 1
+ S[p] = trough
+ peak = x[i]
+ goto = MAXZO
+ continue
+ elif goto == MAXZO:
+ f = 1
+ while goto == MAXZO:
+ i += 1
+ if i == L:
+ goto = ENDZO
+ continue
+ else:
+ if x[i] > peak:
+ peak = x[i]
+ else:
+ if peak - x[i] >= R:
+ p += 1
+ S[p] = peak
+ trough = x[i]
+ goto = MINZO
+ continue
+ elif goto == ENDZO:
+
+ n = p + 1
+ if abs(f) == 1:
+ if f == 1:
+ S[n] = peak
+ else:
+ S[n] = trough
+ else:
+ S[n] = (trough + peak) / 2
+ S = S[1:n + 1]
+ return S
diff --git a/wetb/fatigue_tools/peak_trough.pyd b/wetb/fatigue_tools/peak_trough.pyd
new file mode 100644
index 0000000000000000000000000000000000000000..db4087cd89878143c3de4f75a60b5026f5dceb1f
Binary files /dev/null and b/wetb/fatigue_tools/peak_trough.pyd differ
diff --git a/wetb/fatigue_tools/rainflowcount_astm.py b/wetb/fatigue_tools/rainflowcount_astm.py
new file mode 100644
index 0000000000000000000000000000000000000000..06924c526cb5007e3e50ed9a1e678c41b60a8f3c
--- /dev/null
+++ b/wetb/fatigue_tools/rainflowcount_astm.py
@@ -0,0 +1,102 @@
+'''
+Created on 27/02/2013
+
+@author: mmpe
+
+
+
+How to use:
+
+import_cython("cy_rainflowcount",'cy_rainflowcount.py','')
+from cy_rainflowcount import find_extremes,rainflow
+
+ext = find_extremes(np.array([-2,0,1,0,-3,0,5,0,-1,0,3,0,-4,0,4,0,-2]).astype(np.double))
+print rainflow(ext)
+'''
+
+import numpy as np
+
+
+def find_extremes(signal): #cpdef find_extremes(np.ndarray[double,ndim=1] signal):
+ """return indexes of local minima and maxima plus first and last element of signal"""
+
+ #cdef int pi, i
+ # sign of gradient
+ sign_grad = np.int8(np.sign(np.diff(signal)))
+
+ # remove plateaus(sign_grad==0) by sign_grad[plateau_index]=sign_grad[plateau_index-1]
+ plateau_indexes, = np.where(sign_grad == 0)
+ if len(plateau_indexes) > 0 and plateau_indexes[0] == 0:
+ # first element is a plateau
+ if len(plateau_indexes) == len(sign_grad):
+ # All values are equal to crossing level!
+ return np.array([0])
+
+ # set first element = first element which is not a plateau and delete plateau index
+ i = 0
+ while sign_grad[i] == 0:
+ i += 1
+ sign_grad[0] = sign_grad[i]
+
+ plateau_indexes = np.delete(plateau_indexes, 0)
+
+ for pi in plateau_indexes.tolist():
+ sign_grad[pi] = sign_grad[pi - 1]
+
+ extremes, = np.where(np.r_[1, (sign_grad[1:] * sign_grad[:-1] < 0), 1])
+
+ return signal[extremes]
+
+
+def rainflowcount(sig): #cpdef rainflowcount(np.ndarray[double,ndim=1] sig):
+ """Cython compilable rain ampl_mean count without time analysis
+
+
+ This implemementation is based on the c-implementation by Adam Nieslony found at
+ the MATLAB Central File Exchange http://www.mathworks.com/matlabcentral/fileexchange/3026
+
+ References
+ ----------
+ Adam Nieslony, "Determination of fragments of multiaxial service loading
+ strongly influencing the fatigue of machine components,"
+ Mechanical Systems and Signal Processing 23, no. 8 (2009): 2712-2721.
+
+ and is based on the following standard:
+ ASTM E 1049-85 (Reapproved 1997), Standard practices for cycle counting in
+ fatigue analysis, in: Annual Book of ASTM Standards, vol. 03.01, ASTM,
+ Philadelphia, 1999, pp. 710-718.
+
+ Copyright (c) 1999-2002 by Adam Nieslony
+
+ Ported to Cython compilable Python by Mads M Pedersen
+ In addition peak amplitude is changed to peak to peak amplitude
+
+
+ """
+
+ #cdef int sig_ptr, index
+ #cdef double ampl
+ a = []
+ sig_ptr = 0
+ ampl_mean = []
+ for _ in range(len(sig)):
+ a.append(sig[sig_ptr])
+ sig_ptr += 1
+ while len(a) > 2 and abs(a[-3] - a[-2]) <= abs(a[-2] - a[-1]):
+ ampl = abs(a[-3] - a[-2])
+ mean = (a[-3] + a[-2]) / 2;
+ if len(a) == 3:
+ del a[0]
+ if ampl > 0:
+ ampl_mean.append((ampl, mean))
+ elif len(a) > 3:
+ del a[-3:-1]
+ if ampl > 0:
+ ampl_mean.append((ampl, mean))
+ ampl_mean.append((ampl, mean))
+ for index in range(len(a) - 1):
+ ampl = abs(a[index] - a[index + 1])
+ mean = (a[index] + a[index + 1]) / 2;
+ if ampl > 0:
+ ampl_mean.append((ampl, mean))
+ return ampl_mean
diff --git a/wetb/fatigue_tools/rainflowcount_astm.pyd b/wetb/fatigue_tools/rainflowcount_astm.pyd
new file mode 100644
index 0000000000000000000000000000000000000000..6ef6cd427dfb019332137805ebe4db80d90dcbc8
Binary files /dev/null and b/wetb/fatigue_tools/rainflowcount_astm.pyd differ
diff --git a/wetb/fatigue_tools/tests/__init__.py b/wetb/fatigue_tools/tests/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/wetb/fatigue_tools/tests/test_fatigue.py b/wetb/fatigue_tools/tests/test_fatigue.py
new file mode 100644
index 0000000000000000000000000000000000000000..7949add8d12e7334b7a08730dfb0f3729d50dd47
--- /dev/null
+++ b/wetb/fatigue_tools/tests/test_fatigue.py
@@ -0,0 +1,67 @@
+'''
+Created on 16/07/2013
+
+@author: mmpe
+'''
+
+import unittest
+
+import numpy as np
+from wetb.fatigue_tools.fatigue import eq_load, rainflow_astm, rainflow_windap, \
+ cycle_matrix
+from wetb.hawc2 import Hawc2io
+
+
+class Test(unittest.TestCase):
+
+
+ def test_astm1(self):
+
+ signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+
+ ampl, mean = rainflow_astm(signal)
+ np.testing.assert_array_equal(np.histogram2d(ampl, mean, [6, 4])[0], np.array([[ 0., 1., 0., 0.],
+ [ 1., 0., 0., 2.],
+ [ 0., 0., 0., 0.],
+ [ 0., 0., 0., 1.],
+ [ 0., 0., 0., 0.],
+ [ 0., 0., 1., 2.]]))
+
+ def test_windap1(self):
+ signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
+ ampl, mean = rainflow_windap(signal, 18, 2)
+ np.testing.assert_array_equal(np.histogram2d(ampl, mean, [6, 4])[0], np.array([[ 0., 0., 1., 0.],
+ [ 1., 0., 0., 2.],
+ [ 0., 0., 0., 0.],
+ [ 0., 0., 0., 1.],
+ [ 0., 0., 0., 0.],
+ [ 0., 0., 2., 1.]]))
+
+ def test_windap2(self):
+ data = Hawc2io.ReadHawc2("test_files/test").ReadBinary([2]).flatten()
+ np.testing.assert_allclose(eq_load(data, neq=61), np.array([[1.356, 1.758, 2.370, 2.784, 3.077, 3.296]]), 0.01)
+
+
+ def test_astm2(self):
+ data = Hawc2io.ReadHawc2("test_files/test").ReadBinary([2]).flatten()
+ np.testing.assert_allclose(eq_load(data, neq=61, rainflow_func=rainflow_astm), np.array([[1.356, 1.758, 2.370, 2.784, 3.077, 3.296]]), 0.01)
+
+
+ def test_windap3(self):
+ data = Hawc2io.ReadHawc2("test_files/test").ReadBinary([2]).flatten()
+ np.testing.assert_array_equal(cycle_matrix(data, 4, 4, rainflow_func=rainflow_windap)[0], np.array([[ 14., 65., 39., 24.],
+ [ 0., 1., 4., 0.],
+ [ 0., 0., 0., 0.],
+ [ 0., 1., 2., 0.]]) / 2)
+
+
+ def test_astm3(self):
+ data = Hawc2io.ReadHawc2("test_files/test").ReadBinary([2]).flatten()
+ np.testing.assert_allclose(cycle_matrix(data, 4, 4, rainflow_func=rainflow_astm)[0], np.array([[ 24., 83., 53., 26.],
+ [ 0., 1., 4., 0.],
+ [ 0., 0., 0., 0.],
+ [ 0., 1., 2., 0.]]) / 2, 0.001)
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.testName']
+ unittest.main()
diff --git a/wetb/fatigue_tools/tests/test_files/test.dat b/wetb/fatigue_tools/tests/test_files/test.dat
new file mode 100644
index 0000000000000000000000000000000000000000..60c7efd914a0cf82ef78a7f21777d14c1bf40422
Binary files /dev/null and b/wetb/fatigue_tools/tests/test_files/test.dat differ
diff --git a/wetb/fatigue_tools/tests/test_files/test.sel b/wetb/fatigue_tools/tests/test_files/test.sel
new file mode 100644
index 0000000000000000000000000000000000000000..185e5253c0e5f1947b6a9d22079811714bc79206
--- /dev/null
+++ b/wetb/fatigue_tools/tests/test_files/test.sel
@@ -0,0 +1,114 @@
+________________________________________________________________________________________________________________________
+ Version ID : Pydap 0.1
+ Time : 15:04:16
+ Date : 01:08.2012
+________________________________________________________________________________________________________________________
+ Result file : unit_test/test_files/test.dat
+________________________________________________________________________________________________________________________
+ Scans Channels Time [sec] Format
+ 2440 50 61.000 BINARY
+
+ Channel Variable Description
+
+ 1 Time s Time
+ 2 WSP gl. coo.,Vy m/s Free wind speed Vy, gl. coo, of gl. pos 0.75, 0.00, -40.75
+ 3 WSP gl. coo.,Vy m/s Free wind speed Vy, gl. coo, of gl. pos 1.00, 0.00, -39.00
+ 4 WSP gl. coo.,Vy m/s Free wind speed Vy, gl. coo, of gl. pos 7.00, 0.00, -33.00
+ 5 Omega rad/s Rotor speed
+ 6 Ae pos x , glco, R= 20.5 m Aero position x of blade 1 at radius 20.50, global coo.
+ 7 Ae pos y , glco, R= 20.5 m Aero position y of blade 1 at radius 20.50, global coo.
+ 8 Ae pos z , glco, R= 20.5 m Aero position z of blade 1 at radius 20.50, global coo.
+ 9 Ae pos x , glco, R= 20.5 m Aero position x of blade 2 at radius 20.50, global coo.
+ 10 Ae pos y , glco, R= 20.5 m Aero position y of blade 2 at radius 20.50, global coo.
+ 11 Ae pos z , glco, R= 20.5 m Aero position z of blade 2 at radius 20.50, global coo.
+ 12 Ae pos x , glco, R= 20.5 m Aero position x of blade 3 at radius 20.50, global coo.
+ 13 Ae pos y , glco, R= 20.5 m Aero position y of blade 3 at radius 20.50, global coo.
+ 14 Ae pos z , glco, R= 20.5 m Aero position z of blade 3 at radius 20.50, global coo.
+ 15 Ae pos x , glco, R= 16.5 m Aero position x of blade 1 at radius 16.49, global coo.
+ 16 Ae pos y , glco, R= 16.5 m Aero position y of blade 1 at radius 16.49, global coo.
+ 17 Ae pos z , glco, R= 16.5 m Aero position z of blade 1 at radius 16.49, global coo.
+ 18 Ae pos x , glco, R= 16.5 m Aero position x of blade 2 at radius 16.49, global coo.
+ 19 Ae pos y , glco, R= 16.5 m Aero position y of blade 2 at radius 16.49, global coo.
+ 20 Ae pos z , glco, R= 16.5 m Aero position z of blade 2 at radius 16.49, global coo.
+ 21 Ae pos x , glco, R= 16.5 m Aero position x of blade 3 at radius 16.49, global coo.
+ 22 Ae pos y , glco, R= 16.5 m Aero position y of blade 3 at radius 16.49, global coo.
+ 23 Ae pos z , glco, R= 16.5 m Aero position z of blade 3 at radius 16.49, global coo.
+ 24 Ae pos x , glco, R= 11.6 m Aero position x of blade 1 at radius 11.59, global coo.
+ 25 Ae pos y , glco, R= 11.6 m Aero position y of blade 1 at radius 11.59, global coo.
+ 26 Ae pos z , glco, R= 11.6 m Aero position z of blade 1 at radius 11.59, global coo.
+ 27 Ae pos x , glco, R= 11.6 m Aero position x of blade 2 at radius 11.59, global coo.
+ 28 Ae pos y , glco, R= 11.6 m Aero position y of blade 2 at radius 11.59, global coo.
+ 29 Ae pos z , glco, R= 11.6 m Aero position z of blade 2 at radius 11.59, global coo.
+ 30 Ae pos x , glco, R= 11.6 m Aero position x of blade 3 at radius 11.59, global coo.
+ 31 Ae pos y , glco, R= 11.6 m Aero position y of blade 3 at radius 11.59, global coo.
+ 32 Ae pos z , glco, R= 11.6 m Aero position z of blade 3 at radius 11.59, global coo.
+ 33 Ae pos x , glco, R= 7.6 m Aero position x of blade 1 at radius 7.60, global coo.
+ 34 Ae pos y , glco, R= 7.6 m Aero position y of blade 1 at radius 7.60, global coo.
+ 35 Ae pos z , glco, R= 7.6 m Aero position z of blade 1 at radius 7.60, global coo.
+ 36 Ae pos x , glco, R= 7.6 m Aero position x of blade 2 at radius 7.60, global coo.
+ 37 Ae pos y , glco, R= 7.6 m Aero position y of blade 2 at radius 7.60, global coo.
+ 38 Ae pos z , glco, R= 7.6 m Aero position z of blade 2 at radius 7.60, global coo.
+ 39 Ae pos x , glco, R= 7.6 m Aero position x of blade 3 at radius 7.60, global coo.
+ 40 Ae pos y , glco, R= 7.6 m Aero position y of blade 3 at radius 7.60, global coo.
+ 41 Ae pos z , glco, R= 7.6 m Aero position z of blade 3 at radius 7.60, global coo.
+ 42 Ae pos x , glco, R= 4.0 m Aero position x of blade 1 at radius 4.01, global coo.
+ 43 Ae pos y , glco, R= 4.0 m Aero position y of blade 1 at radius 4.01, global coo.
+ 44 Ae pos z , glco, R= 4.0 m Aero position z of blade 1 at radius 4.01, global coo.
+ 45 Ae pos x , glco, R= 4.0 m Aero position x of blade 2 at radius 4.01, global coo.
+ 46 Ae pos y , glco, R= 4.0 m Aero position y of blade 2 at radius 4.01, global coo.
+ 47 Ae pos z , glco, R= 4.0 m Aero position z of blade 2 at radius 4.01, global coo.
+ 48 Ae pos x , glco, R= 4.0 m Aero position x of blade 3 at radius 4.01, global coo.
+ 49 Ae pos y , glco, R= 4.0 m Aero position y of blade 3 at radius 4.01, global coo.
+ 50 Ae pos z , glco, R= 4.0 m Aero position z of blade 3 at radius 4.01, global coo.
+________________________________________________________________________________________________________________________
+Scale factors:
+ 1.90625E-03
+ 3.86850E-04
+ 3.93534E-04
+ 3.63391E-04
+ 1.56250E-05
+ 6.40656E-04
+ 0.00000E+00
+ 1.89066E-03
+ 6.40656E-04
+ 0.00000E+00
+ 1.89066E-03
+ 6.40656E-04
+ 0.00000E+00
+ 1.89066E-03
+ 5.15428E-04
+ 6.87809E-07
+ 1.76543E-03
+ 5.15428E-04
+ 6.87809E-07
+ 1.76543E-03
+ 5.15428E-04
+ 6.87809E-07
+ 1.76543E-03
+ 3.62459E-04
+ 1.85211E-06
+ 1.61246E-03
+ 3.62459E-04
+ 1.85211E-06
+ 1.61246E-03
+ 3.62459E-04
+ 1.85211E-06
+ 1.61246E-03
+ 2.38552E-04
+ 5.23494E-06
+ 1.48855E-03
+ 2.38552E-04
+ 5.23494E-06
+ 1.48855E-03
+ 2.38551E-04
+ 5.23494E-06
+ 1.48855E-03
+ 1.29018E-04
+ 1.14916E-05
+ 1.37902E-03
+ 1.29018E-04
+ 1.14916E-05
+ 1.37902E-03
+ 1.29018E-04
+ 1.14916E-05
+ 1.37902E-03