wetb/fast/tests/test_fast_io.py

@@ -35,6 +35,13 @@ class TestFastIO(unittest.TestCase):
         self.assertEqual(info['attribute_units'][7], 'deg/s^2')
         self.assertAlmostEqual(data[10, 4], 138.822277739535)
 
+    def test_load_output2(self):
+        data, info = load_output(testfilepath + 'DTU10MW.out')
+        self.assertEqual(info['name'], "DTU10MW")
+        self.assertEqual(info['attribute_names'][1], "RotPwr")
+        self.assertEqual(info['attribute_units'][1], "kW")
+
+
 
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testload_output']

wetb/fatigue_tools/fatigue.py

@@ -19,6 +19,7 @@ from __future__ import print_function
 from __future__ import unicode_literals
 from __future__ import absolute_import
 from future import standard_library
+import warnings
 standard_library.install_aliases()
 import numpy as np
 from wetb.fatigue_tools.rainflowcounting import rainflowcount
@@ -70,7 +71,6 @@ def eq_load(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1, rainflow_func=ra
         return [[np.nan] * len(np.atleast_1d(m))] * len(np.atleast_1d(neq))
 
 
-
 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
 
@@ -95,7 +95,7 @@ def eq_load_and_cycles(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=[10 ** 6
     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)
+        2d array with shape = (no_ampl_bins, 1)
     ampl_bin_mean : array_like
         mean amplitude of the bins
     ampl_bin_edges
@@ -104,12 +104,13 @@ def eq_load_and_cycles(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=[10 ** 6
     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)]
+    cycles, ampl_bin_mean = cycles.flatten(), ampl_bin_mean.flatten()
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")   
+        eq_loads = [[((np.nansum(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
 
@@ -118,8 +119,8 @@ def cycle_matrix(signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_win
     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)]
+        - if array-like, the raw signal\n
+        - if list of tuples, list of (weight, signal), e.g. [(0.1,sig1), (0.8,sig2), (.1,sig3)]\n
     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
@@ -132,7 +133,8 @@ def cycle_matrix(signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_win
     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.
+        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,)
@@ -150,28 +152,19 @@ def cycle_matrix(signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_win
     """
 
     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)]
+        weights, ampls, means = np.array([(np.zeros_like(ampl)+weight,ampl,mean) for weight, signal in signals for ampl,mean in rainflow_func(signal[:]).T], dtype=np.float64).T
     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()
+        ampl_bins = np.linspace(0, 1, num=ampl_bins + 1) * ampls[weights>0].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]
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")     
+        ampl_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * ampls)[0]
+        ampl_bin_mean = np.nanmean(ampl_bin_sum / np.where(cycles,cycles,np.nan),1)
+        mean_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * means)[0]
+        mean_bin_mean = np.nanmean(mean_bin_sum / np.where(cycles, cycles, np.nan), 1)
     cycles = cycles / 2  # to get full cycles
     return cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges
 

wetb/fatigue_tools/rainflowcounting/compile.py

-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import division
-from __future__ import absolute_import
-from future import standard_library
-standard_library.install_aliases()
-from wetb.utils.cython_compile.cython_compile import cython_import
-cython_import('pair_range')
-cython_import('peak_trough')
-cython_import('rainflowcount_astm')

wetb/fatigue_tools/rainflowcounting/pair_range.pyx

+import cython
+import numpy as np
+cimport numpy as np
+
+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/rainflowcounting/peak_trough.pyx

+import cython
+import numpy as np
+cimport numpy as np
+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)
+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/rainflowcounting/rainflowcount.py

@@ -6,7 +6,8 @@ from builtins import str
 from future import standard_library
 standard_library.install_aliases()
 import numpy as np
-from wetb.utils.cython_compile.cython_compile import cython_import
+from wetb.fatigue_tools.rainflowcounting import peak_trough
+from wetb.fatigue_tools.rainflowcounting import pair_range
 
 
 def check_signal(signal):
@@ -63,7 +64,8 @@ def rainflow_windap(signal, levels=255., thresshold=(255 / 50)):
     check_signal(signal)
     #type <double> is required by <find_extreme> and <rainflow>
     signal = signal.astype(np.double)
-
+    if np.all(np.isnan(signal)):
+        return None
     offset = np.nanmin(signal)
     signal -= offset
     if np.nanmax(signal) > 0:
@@ -73,19 +75,14 @@ def rainflow_windap(signal, levels=255., thresshold=(255 / 50)):
 
 
         # If possible the module is compiled using cython otherwise the python implementation is used
-        cython_import('wetb.fatigue_tools.rainflowcounting.peak_trough')
-        cython_import('wetb.fatigue_tools.rainflowcounting.pair_range')
-
-        from wetb.fatigue_tools.rainflowcounting.peak_trough import peak_trough
-        from wetb.fatigue_tools.rainflowcounting.pair_range import pair_range_amplitude_mean
 
 
         #Convert to list of local minima/maxima where difference > thresshold
-        sig_ext = peak_trough(signal, thresshold)
+        sig_ext = peak_trough.peak_trough(signal, thresshold)
 
 
         #rainflow count
-        ampl_mean = pair_range_amplitude_mean(sig_ext)
+        ampl_mean = pair_range.pair_range_amplitude_mean(sig_ext)
 
         ampl_mean = np.array(ampl_mean)
         ampl_mean = np.round(ampl_mean / thresshold) * gain * thresshold
@@ -129,7 +126,7 @@ def rainflow_astm(signal):
 
     # Import find extremes and rainflow.
     # If possible the module is compiled using cython otherwise the python implementation is used
-    cython_import('wetb.fatigue_tools.rainflowcounting.rainflowcount_astm')
+    
     from wetb.fatigue_tools.rainflowcounting.rainflowcount_astm import find_extremes, rainflowcount
 
     # Remove points which is not local minimum/maximum

wetb/fatigue_tools/rainflowcounting/rainflowcount_astm.pyx

+import cython
+import numpy as np
+cimport numpy as np
+'''
+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
+
+
+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]
+
+
+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/rainflowcounting/rfc_hist.py

-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import absolute_import
-from future import standard_library
-standard_library.install_aliases()
-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

wetb/fatigue_tools/tests/test_bearing_damage.py

@@ -8,6 +8,7 @@ from __future__ import print_function
 from __future__ import division
 from __future__ import absolute_import
 from future import standard_library
+from wetb import gtsdf
 standard_library.install_aliases()
 
 
@@ -26,8 +27,8 @@ class TestBearingDamage(unittest.TestCase):
         self.tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
 
     def test_bearing_damage_swp(self):
-        data = Hawc2io.ReadHawc2(self.tfp + "test_bearing_damage").ReadBinary((np.array([4, 26, 6, 32, 8, 38]) - 1).tolist())
-        self.assertAlmostEqual(bearing_damage([(data[:, i], data[:, i + 1]) for i in [0, 2, 4]]), 7.755595081475002e+13)
+        data = Hawc2io.ReadHawc2(self.tfp + "test_bearing_damage").ReadBinary((np.array([1, 4, 2, 5, 3, 6]) ).tolist())
+        self.assertAlmostEqual(bearing_damage([(data[:, i], data[:, i + 1]) for i in [0,2,4]]), 7.755595081475002e+13)
 
 
 if __name__ == "__main__":

wetb/fatigue_tools/tests/test_fatigue.py

@@ -8,6 +8,7 @@ from __future__ import unicode_literals
 from __future__ import print_function
 from __future__ import absolute_import
 from future import standard_library
+import sys
 standard_library.install_aliases()
 
 import unittest
@@ -22,6 +23,36 @@ testfilepath = os.path.join(os.path.dirname(__file__), 'test_files/')  # test fi
 
 class TestFatigueTools(unittest.TestCase):
 
+    def test_leq_1hz(self):
+        """Simple test of wetb.fatigue_tools.fatigue.eq_load using a sine
+        signal.
+        """
+        amplitude = 1
+        m = 1
+        point_per_deg = 100
+
+        # sine signal with 10 periods (20 peaks)
+        nr_periods = 10
+        time = np.linspace(0, nr_periods*2*np.pi, point_per_deg*180)
+        neq = time[-1]
+        # mean value of the signal shouldn't matter
+        signal = amplitude * np.sin(time) + 5
+        r_eq_1hz = eq_load(signal, no_bins=1, m=m, neq=neq)[0]
+        r_eq_1hz_expected = ((2*nr_periods*amplitude**m)/neq)**(1/m)
+        np.testing.assert_allclose(r_eq_1hz, r_eq_1hz_expected)
+
+        # sine signal with 20 periods (40 peaks)
+        nr_periods = 20
+        time = np.linspace(0, nr_periods*2*np.pi, point_per_deg*180)
+        neq = time[-1]
+        # mean value of the signal shouldn't matter
+        signal = amplitude * np.sin(time) + 9
+        r_eq_1hz2 = eq_load(signal, no_bins=1, m=m, neq=neq)[0]
+        r_eq_1hz_expected2 = ((2*nr_periods*amplitude**m)/neq)**(1/m)
+        np.testing.assert_allclose(r_eq_1hz2, r_eq_1hz_expected2)
+
+        # 1hz equivalent load should be independent of the length of the signal
+        np.testing.assert_allclose(r_eq_1hz, r_eq_1hz2)
 
     def test_astm1(self):
 
@@ -55,12 +86,15 @@ class TestFatigueTools(unittest.TestCase):
         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(testfilepath + "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_windap3(self):
+#         data = Hawc2io.ReadHawc2(testfilepath + "test").ReadBinary([2]).flatten()
+#         from wetb.fatigue_tools.rainflowcounting import peak_trough
+#         self.assertTrue(peak_trough.__file__.lower()[-4:] == ".pyd" or peak_trough.__file__.lower()[-3:] == ".so",
+#                         "not compiled, %s, %s\n%s"%(sys.executable, peak_trough.__file__, os.listdir(os.path.dirname(peak_trough.__file__))))
+#         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):
@@ -70,6 +104,13 @@ class TestFatigueTools(unittest.TestCase):
                                                                                                            [  0., 0., 0., 0.],
                                                                                                            [  0., 1., 2., 0.]]) / 2, 0.001)
 
+    def test_astm_weighted(self):
+        data = Hawc2io.ReadHawc2(testfilepath + "test").ReadBinary([2]).flatten()
+        np.testing.assert_allclose(cycle_matrix([(1, data),(1,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.]]) , 0.001)
+
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testName']
     unittest.main()

wetb/flex/__init__.py

+from ._io import load, read_sensor_info
+from wetb import gtsdf
+import numpy as np
+
+class Dataset(gtsdf.Dataset):
+    def __init__(self, filename, name_stop=8,dtype=np.float):
+        self.time, self.data, self.info = load(filename, name_stop,dtype)
\ No newline at end of file

wetb/flex/tests/test_files/test1/sensor

+Sensor list for C:/mmpe/programming/python/WindEnergyToolbox/wetb/flex/tests/test_files/test1/test.int
+ No   forst  offset  korr. c  Volt    Unit   Navn    Beskrivelse------------
+  1  1.000   0.000      1.00     0.00     m/s WSP_gl._ Free wind speed Vy, gl. coo, of gl. pos    0.75,   0.00, -40.75
+  2  1.000   0.000      1.00     0.00     m/s WSP_gl._ Free wind speed Vy, gl. coo, of gl. pos    1.00,   0.00, -39.00
+  3  1.000   0.000      1.00     0.00     m/s WSP_gl._ Free wind speed Vy, gl. coo, of gl. pos    7.00,   0.00, -33.00
+  4  1.000   0.000      1.00     0.00   rad/s Omega    Rotor speed
+  5  1.000   0.000      1.00     0.00       m Ae_pos_x Aero position x  of blade  1 at radius  20.50, global coo.
+  6  1.000   0.000      1.00     0.00       m Ae_pos_y Aero position y  of blade  1 at radius  20.50, global coo.
+  7  1.000   0.000      1.00     0.00       m Ae_pos_z Aero position z  of blade  1 at radius  20.50, global coo.

wetb/flex/tests/test_files/test_sensor_info/sensor

+  Version ID : HAWC2MB 12.3 (beta -rev_TJ6)
+ No   forst  offset  korr. c  Volt    Unit   Navn    Beskrivelse---------------
+   1    1.0    0.0    0.00      1.0    s      Time    Time
+   2    1.0    0.0    0.00      1.0    deg    bea1 an shaft_rot angle
+   3    1.0    0.0    0.00      1.0    rpm    bea1 an shaft_rot angle speed
+   4    1.0    0.0    0.00      1.0    deg    bea2 an pitch1 angle
+   5    1.0    0.0    0.00      1.0    deg/s  bea2 an pitch1 angle speed
+   6    1.0    0.0    0.00      1.0    deg    bea2 an pitch2 angle
+   7    1.0    0.0    0.00      1.0    deg/s  bea2 an pitch2 angle speed
+   8    1.0    0.0    0.00      1.0    deg    bea2 an pitch3 angle
+   9    1.0    0.0    0.00      1.0    deg/s  bea2 an pitch3 angle speed
+  10    1.0    0.0    0.00      1.0    kNm    Mx coo: MomentMx Mbdy:towertop nodenr:   1 coo: tower  tower top -1: below top mass
+  11    1.0    0.0    0.00      1.0    kNm    My coo: MomentMy Mbdy:towertop nodenr:   1 coo: tower  tower top -1: below top mass
+  12    1.0    0.0    0.00      1.0    kNm    Mz coo: MomentMz Mbdy:towertop nodenr:   1 coo: tower  tower top -1: below top mass
+  13    1.0    0.0    0.00      1.0    kN     Fx coo: Force  Fx Mbdy:towertop nodenr:   1 coo: tower  tower top -1: below top mass
+  14    1.0    0.0    0.00      1.0    kN     Fy coo: Force  Fy Mbdy:towertop nodenr:   1 coo: tower  tower top -1: below top mass
+  15    1.0    0.0    0.00      1.0    kN     Fz coo: Force  Fz Mbdy:towertop nodenr:   1 coo: tower  tower top -1: below top mass
+  16    1.0    0.0    0.00      1.0    kNm    Mx coo: MomentMx Mbdy:tower nodenr:   1 coo: tower  tower base flange
+  17    1.0    0.0    0.00      1.0    kNm    My coo: MomentMy Mbdy:tower nodenr:   1 coo: tower  tower base flange
+  18    1.0    0.0    0.00      1.0    kNm    Mz coo: MomentMz Mbdy:tower nodenr:   1 coo: tower  tower base flange
+  19    1.0    0.0    0.00      1.0    kNm    Mx coo: MomentMx Mbdy:towertop nodenr:   2 coo: towertop  yaw bearing
+  20    1.0    0.0    0.00      1.0    kNm    My coo: MomentMy Mbdy:towertop nodenr:   2 coo: towertop  yaw bearing
+  21    1.0    0.0    0.00      1.0    kNm    Mz coo: MomentMz Mbdy:towertop nodenr:   2 coo: towertop  yaw bearing
+  22    1.0    0.0    0.00      1.0    kN     Fx coo: Force  Fx Mbdy:towertop nodenr:   2 coo: towertop  yaw bering
+  23    1.0    0.0    0.00      1.0    kN     Fy coo: Force  Fy Mbdy:towertop nodenr:   2 coo: towertop  yaw bering
+  24    1.0    0.0    0.00      1.0    kN     Fz coo: Force  Fz Mbdy:towertop nodenr:   2 coo: towertop  yaw bering
+  25    1.0    0.0    0.00      1.0    kNm    Mx coo: MomentMx Mbdy:shaft nodenr:   4 coo: shaft  main bearing
+  26    1.0    0.0    0.00      1.0    kNm    My coo: MomentMy Mbdy:shaft nodenr:   4 coo: shaft  main bearing
+  27    1.0    0.0    0.00      1.0    kNm    Mz coo: MomentMz Mbdy:shaft nodenr:   4 coo: shaft  main bearing
+  28    1.0    0.0    0.00      1.0    kNm    Mx coo: MomentMx Mbdy:blade1 nodenr:   3 coo: blade1  blade 1 root
+  29    1.0    0.0    0.00      1.0    kNm    My coo: MomentMy Mbdy:blade1 nodenr:   3 coo: blade1  blade 1 root
+  30    1.0    0.0    0.00      1.0    kNm    Mz coo: MomentMz Mbdy:blade1 nodenr:   3 coo: blade1  blade 1 root
+  31    1.0    0.0    0.00      1.0    kNm    Mx coo: MomentMx Mbdy:blade1 nodenr:  10 coo: local  blade 1 50% local e coo
+  32    1.0    0.0    0.00      1.0    kNm    My coo: MomentMy Mbdy:blade1 nodenr:  10 coo: local  blade 1 50% local e coo
+  33    1.0    0.0    0.00      1.0    kNm    Mz coo: MomentMz Mbdy:blade1 nodenr:  10 coo: local  blade 1 50% local e coo
+  34    1.0    0.0    0.00      1.0    kN     Fx coo: Force  Fx Mbdy:blade1 nodenr:   1 coo: blade1  blade 1 root
+  35    1.0    0.0    0.00      1.0    kN     Fy coo: Force  Fy Mbdy:blade1 nodenr:   1 coo: blade1  blade 1 root
+  36    1.0    0.0    0.00      1.0    kN     Fz coo: Force  Fz Mbdy:blade1 nodenr:   1 coo: blade1  blade 1 root
+  37    1.0    0.0    0.00      1.0    kNm    Mx coo: MomentMx Mbdy:blade2 nodenr:   1 coo: blade2  blade 2 root
+  38    1.0    0.0    0.00      1.0    kNm    My coo: MomentMy Mbdy:blade2 nodenr:   1 coo: blade2  blade 2 root
+  39    1.0    0.0    0.00      1.0    kNm    Mz coo: MomentMz Mbdy:blade2 nodenr:   1 coo: blade2  blade 2 root
+  40    1.0    0.0    0.00      1.0    kNm    Mx coo: MomentMx Mbdy:blade3 nodenr:   1 coo: blade3  blade 3 root
+  41    1.0    0.0    0.00      1.0    kNm    My coo: MomentMy Mbdy:blade3 nodenr:   1 coo: blade3  blade 3 root
+  42    1.0    0.0    0.00      1.0    kNm    Mz coo: MomentMz Mbdy:blade3 nodenr:   1 coo: blade3  blade 3 root
+  43    1.0    0.0    0.00      1.0    m      State p State pos x  Mbdy:towertop E-nr:   1 Z-rel:1.00 coo: global  tower top flange position
+  44    1.0    0.0    0.00      1.0    m      State p State pos y  Mbdy:towertop E-nr:   1 Z-rel:1.00 coo: global  tower top flange position
+  45    1.0    0.0    0.00      1.0    m      State p State pos z  Mbdy:towertop E-nr:   1 Z-rel:1.00 coo: global  tower top flange position
+  46    1.0    0.0    0.00      1.0    m      State p State pos x  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: blade1  blade 1 tip pos
+  47    1.0    0.0    0.00      1.0    m      State p State pos y  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: blade1  blade 1 tip pos
+  48    1.0    0.0    0.00      1.0    m      State p State pos z  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: blade1  blade 1 tip pos
+  49    1.0    0.0    0.00      1.0    m      State p State pos x  Mbdy:blade2 E-nr:  18 Z-rel:1.00 coo: blade2  blade 2 tip pos
+  50    1.0    0.0    0.00      1.0    m      State p State pos y  Mbdy:blade2 E-nr:  18 Z-rel:1.00 coo: blade2  blade 2 tip pos
+  51    1.0    0.0    0.00      1.0    m      State p State pos z  Mbdy:blade2 E-nr:  18 Z-rel:1.00 coo: blade2  blade 2 tip pos
+  52    1.0    0.0    0.00      1.0    m      State p State pos x  Mbdy:blade3 E-nr:  18 Z-rel:1.00 coo: blade3  blade 3 tip pos
+  53    1.0    0.0    0.00      1.0    m      State p State pos y  Mbdy:blade3 E-nr:  18 Z-rel:1.00 coo: blade3  blade 3 tip pos
+  54    1.0    0.0    0.00      1.0    m      State p State pos z  Mbdy:blade3 E-nr:  18 Z-rel:1.00 coo: blade3  blade 3 tip pos
+  55    1.0    0.0    0.00      1.0    m      State p State pos x  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: global  blade 1 tip pos
+  56    1.0    0.0    0.00      1.0    m      State p State pos y  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: global  blade 1 tip pos
+  57    1.0    0.0    0.00      1.0    m      State p State pos z  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: global  blade 1 tip pos
+  58    1.0    0.0    0.00      1.0    m      State p State pos x  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: hub1  blade 1 tip pos hub coo
+  59    1.0    0.0    0.00      1.0    m      State p State pos y  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: hub1  blade 1 tip pos hub coo
+  60    1.0    0.0    0.00      1.0    m      State p State pos z  Mbdy:blade1 E-nr:  18 Z-rel:1.00 coo: hub1  blade 1 tip pos hub coo
+  61    1.0    0.0    0.00      1.0    m      State p State pos x  Mbdy:blade2 E-nr:  18 Z-rel:1.00 coo: hub2  blade 2 tip pos hub coo
+  62    1.0    0.0    0.00      1.0    m      State p State pos y  Mbdy:blade2 E-nr:  18 Z-rel:1.00 coo: hub2  blade 2 tip pos hub coo
+  63    1.0    0.0    0.00      1.0    m      State p State pos z  Mbdy:blade2 E-nr:  18 Z-rel:1.00 coo: hub2  blade 2 tip pos hub coo
+  64    1.0    0.0    0.00      1.0    m      State p State pos x  Mbdy:blade3 E-nr:  18 Z-rel:1.00 coo: hub3  blade 3 tip pos hub coo
+  65    1.0    0.0    0.00      1.0    m      State p State pos y  Mbdy:blade3 E-nr:  18 Z-rel:1.00 coo: hub3  blade 3 tip pos hub coo
+  66    1.0    0.0    0.00      1.0    m      State p State pos z  Mbdy:blade3 E-nr:  18 Z-rel:1.00 coo: hub3  blade 3 tip pos hub coo
+  67    1.0    0.0    0.00      1.0    rad/s  omega t State_rot omega tz Mbdy:shaft E-nr:   4 Z-rel:1.00 coo: shaft

wetb/flex/tests/test_flex_io.py

+'''
+Created on 11. apr. 2017
+
+@author: mmpe
+'''
+import os
+import unittest
+from wetb import flex
+
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')  # test file path
+class Test(unittest.TestCase):
+
+
+    def test_load(self):
+        time, data, info = flex.load(tfp+"test1/test.int")
+        self.assertEqual(data.shape, (800,7))
+        self.assertEqual(info['attribute_names'][0], "WSP_gl._")
+        self.assertEqual(info['attribute_units'][0], "m/s")
+        self.assertEqual(info['attribute_descriptions'][0], "Free wind speed Vy, gl. coo, of gl. pos 0.75, 0.00, -40.75")
+
+        self.assertAlmostEqual(data[0, 1], 12.037,3)
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()
\ No newline at end of file