fatigue + test

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

wetb/fatigue_tools/fatigue.py

+'''
+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
+
+

wetb/fatigue_tools/pair_range.py

+
+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

wetb/fatigue_tools/peak_trough.py

+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

wetb/fatigue_tools/rainflowcount_astm.py

+'''
+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

wetb/fatigue_tools/tests/test_fatigue.py

+'''
+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()

wetb/fatigue_tools/tests/test_files/test.sel

+________________________________________________________________________________________________________________________
+  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