reorganized the fatigue module

README.md

@@ -25,6 +25,7 @@ General Time Series Data Format, a binary hdf5 data format for storing time seri
 
 ### [fatigue_tools](wetb/fatigue_tools)
 - [fatigue](wetb/fatigue_tools/fatigue.py): Rainflow counting, cycle matrix and equvivalent loads
+- [bearing_damage](wetb/fatigue_tools/bearing_damage.py): Calculate a comparable measure of bearing damage
 
 ### [wind](wetb/wind)
 - [shear](wetb/wind/shear.py): Calculate and fit wind shear 

wetb/fatigue_tools/fatigue.py

@@ -2,199 +2,23 @@
 Created on 04/03/2013
 @author: mmpe
 
-'rainflow_windap' and 'rainflow_astm' are two different methods to for rainflow counting
+
 'eq_load' calculate equivalent loads using one of the two rain flow counting methods
 'cycle_matrix' calculates a matrix of cycles (binned on amplitude and mean value)
 'eq_load_and_cycles' is used to calculate eq_loads of multiple time series (e.g. life time equivalent load)
 
+The methods uses the rainflow counting routines (See documentation in top of methods):
+- 'rainflow_windap': (Described in "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads",
+                      2. edition 1990, Appendix A)
+or
+- 'rainflow_astm' (based on the c-implementation by Adam Nieslony found at the MATLAB Central File Exchange
+                   http://www.mathworks.com/matlabcentral/fileexchange/3026)
 '''
-
-#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
+from wetb.fatigue_tools.rainflowcounting import rainflowcount
 
-
-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
+rainflow_windap = rainflowcount.rainflow_windap
+rainflow_astm = rainflowcount.rainflow_astm
 
 
 def eq_load(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1, rainflow_func=rainflow_windap):
@@ -239,14 +63,6 @@ def eq_load(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1, rainflow_func=ra
     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):
@@ -354,3 +170,18 @@ def cycle_matrix(signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_win
     return cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges
 
 
+if __name__ == "__main__":
+    signal1 = 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])
+    signal2 = signal1 * 1.1
+
+    # equivalent load for default whler slopes
+    print (eq_load(signal1, no_bins=50, neq=17, rainflow_func=rainflow_windap))
+    print (eq_load(signal1, no_bins=50, neq=17, rainflow_func=rainflow_astm))
+
+    # Cycle matrix with 4 amplitude bins and 4 mean value bins
+    print (cycle_matrix(signal1, 4, 4, rainflow_func=rainflow_windap))
+    print (cycle_matrix(signal1, 4, 4, rainflow_func=rainflow_astm))
+
+    # Cycle matrix where signal1 and signal2 contributes with 50% each
+    print (cycle_matrix([(.5, signal1), (.5, signal2)], 4, 8, rainflow_func=rainflow_astm))
+

wetb/fatigue_tools/rainflowcounting/compile.py

+import sys
+sys.path.append("../../../../MMPE/")
+from mmpe.cython_compile.cython_compile import cython_import
+
+cython_import('pair_range')
+cython_import('peak_trough')
+cython_import('rainflowcount_astm')

wetb/fatigue_tools/rainflowcounting/rainflowcount.py

+import numpy as np
+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.rainflowcounting.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.rainflowcounting.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.rainflowcounting.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
\ No newline at end of file

wetb/fatigue_tools/rainflowcounting/rfc_hist.py

+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
\ No newline at end of file