wetb/signal/fit/_bin_fit.py

 
 import numpy as np
-from scipy.interpolate.interpolate import interp1d
+from scipy.interpolate import interp1d
 
 
-
-def bin_fit(x,y, bins=10, kind='linear', bin_func=np.nanmean, bin_min_count=3, lower_upper='discard'):
+def bin_fit(x, y, bins=10, kind='linear', bin_func=np.nanmean, bin_min_count=3, lower_upper='discard'):
     """Fit observations based on bin statistics
-    
+
     Parameters
     ---------
     x : array_like
@@ -29,29 +28,29 @@ def bin_fit(x,y, bins=10, kind='linear', bin_func=np.nanmean, bin_min_count=3, l
         - "discard":
         - "extrapolate":
         - int: Set f(max(x)) to mean of first/last int observations
-       
-    
+
+
     Returns
     -------
     bin_x, fit_function
-        
+
     """
     x, y = np.array(x[:]), np.array(y[:])
     if isinstance(bins, int):
         bins = np.linspace(np.nanmin(x), np.nanmax(x) + 1e-10, bins + 1)
-    elif isinstance(bins, tuple) and len(bins)==2 and isinstance(bins[0], int) and isinstance(bins[1], int):
+    elif isinstance(bins, tuple) and len(bins) == 2 and isinstance(bins[0], int) and isinstance(bins[1], int):
         xbins, ybins = bins
-        if xbins>0:
+        if xbins > 0:
             xbinsx = np.linspace(np.nanmin(x), np.nanmax(x) + 1e-10, xbins + 1)
         else:
             xbinsx = []
-        if ybins>0:
-            x1, f1 = bin_fit(y,x, kind=1, bins=ybins)
+        if ybins > 0:
+            x1, f1 = bin_fit(y, x, kind=1, bins=ybins)
             xbinsy = f1(x1)
         else:
             xbinsy = []
         #x2, f2 = bin_fit(x,y, kind=1, bins=xbins)
-        bins = sorted(np.r_[xbinsx, xbinsy ])
+        bins = sorted(np.r_[xbinsx, xbinsy])
 
     digitized = np.digitize(x, bins)
     digitized[np.isnan(x) | np.isnan(y)] = -1
@@ -66,50 +65,50 @@ def bin_fit(x,y, bins=10, kind='linear', bin_func=np.nanmean, bin_min_count=3, l
     #bin_x_fit, bin_y = [b[bin_count >= bin_min_count] for b in [bin_x, bin_y]]
     bin_x_fit = bin_x
     m = np.isnan(bin_x_fit)
-    bin_x_fit[m] = ((bins[:-1]+bins[1:])/2)[m]
+    bin_x_fit[m] = ((bins[:-1] + bins[1:]) / 2)[m]
     bin_y_fit = bin_y.copy()
-    bin_y_fit[bin_count<bin_min_count]= np.nan
+    bin_y_fit[bin_count < bin_min_count] = np.nan
 
-    
     if isinstance(lower_upper, (str, int)):
         lower = upper = lower_upper
     else:
-        lower, upper = lower_upper 
-    
-    #Add value to min(x)
+        lower, upper = lower_upper
+
+    # Add value to min(x)
     if bin_x_fit[0] > np.nanmin(x) or np.isnan(bin_y_fit[0]):
-        if lower =='extrapolate':
-        
-            bin_y_fit = np.r_[bin_y_fit[0] - (bin_x_fit[0] - np.nanmin(x)) * (bin_y_fit[1] - bin_y_fit[0]) / (bin_x_fit[1] - bin_x_fit[0]), bin_y_fit]
+        if lower == 'extrapolate':
+
+            bin_y_fit = np.r_[bin_y_fit[0] - (bin_x_fit[0] - np.nanmin(x)) *
+                              (bin_y_fit[1] - bin_y_fit[0]) / (bin_x_fit[1] - bin_x_fit[0]), bin_y_fit]
             bin_x_fit = np.r_[np.nanmin(x), bin_x_fit]
-        elif lower=="discard":
+        elif lower == "discard":
             pass
         elif isinstance(lower, int):
             bin_y_fit = np.r_[np.mean(y[~np.isnan(x)][np.argsort(x[~np.isnan(x)])[:lower]]), bin_y_fit]
             bin_x_fit = np.r_[np.nanmin(x), bin_x_fit]
         else:
-            raise NotImplementedError("Argument for handling lower observations, %s, not implemented"%lower)
-        
-    #add value to max(x)
+            raise NotImplementedError("Argument for handling lower observations, %s, not implemented" % lower)
+
+    # add value to max(x)
     if bin_x_fit[-1] < np.nanmax(x) or np.isnan(bin_y_fit[-1]):
         if upper == 'extrapolate':
-            bin_y_fit = np.r_[bin_y_fit, bin_y_fit[-1] + (np.nanmax(x) - bin_x_fit[-1]) * (bin_y_fit[-1] - bin_y_fit[-2]) / (bin_x_fit[-1] - bin_x_fit[-2]) ]
+            bin_y_fit = np.r_[bin_y_fit, bin_y_fit[-1] + (np.nanmax(x) - bin_x_fit[-1])
+                              * (bin_y_fit[-1] - bin_y_fit[-2]) / (bin_x_fit[-1] - bin_x_fit[-2])]
             bin_x_fit = np.r_[bin_x_fit, np.nanmax(x)]
-        elif upper=="discard":
+        elif upper == "discard":
             pass
         elif isinstance(upper, int):
             bin_y_fit = np.r_[bin_y_fit, np.mean(y[~np.isnan(x)][np.argsort(x[~np.isnan(x)])[-upper:]])]
             bin_x_fit = np.r_[bin_x_fit, np.nanmax(x)]
         else:
-            raise NotImplementedError("Argument for handling upper observations, %s, not implemented"%upper)
-        
+            raise NotImplementedError("Argument for handling upper observations, %s, not implemented" % upper)
 
-    
     return bin_x_fit, _interpolate_fit(bin_x_fit, bin_y_fit, kind)
-    
-def perpendicular_bin_fit(x, y, bins = 30, fit_func=None, bin_min_count=3, plt=None):
+
+
+def perpendicular_bin_fit(x, y, bins=30, fit_func=None, bin_min_count=3, plt=None):
     """Fit a curve to the values, (x,y) using bins that are perpendicular to an initial fit
-    
+
     Parameters
     ---------
     x : array_like
@@ -124,77 +123,76 @@ def perpendicular_bin_fit(x, y, bins = 30, fit_func=None, bin_min_count=3, plt=N
     bin_min_count : int, optional
         Minimum number of observations in bins to include
         Default is 3
-    
+
     plt : pyplot or None
         If pyplot the fitting process is plotted on plt
-        
+
     Returns
     -------
     fit_x, fit_y
     """
-    
+
     if fit_func is None:
-        fit_func = lambda x,y : bin_fit(x, y, bins, bin_func=np.nanmean)
-    
-    x,y = [v[~np.isnan(x)&~np.isnan(y)] for v in [x,y]]
-    
-    bfx,f = fit_func(x, y)
+        def fit_func(x, y): return bin_fit(x, y, bins, bin_func=np.nanmean)
+
+    x, y = [v[~np.isnan(x) & ~np.isnan(y)] for v in [x, y]]
+
+    bfx, f = fit_func(x, y)
     bfy = f(bfx)
-    
-    
-    bfx, bfy = [v[~np.isnan(bfx)&~np.isnan(bfy)] for v in [bfx,bfy]]
+
+    bfx, bfy = [v[~np.isnan(bfx) & ~np.isnan(bfy)] for v in [bfx, bfy]]
     if plt:
-        x_range, y_range = [v.max()-v.min() for v in [x,y]]
-        plt.ylim([y.min()-y_range*.1, y.max()+y_range*.1])
-        plt.xlim([x.min()-x_range*.1, x.max()+x_range*.1])
-        
-    # divide curve into N segments of same normalized curve length 
-    xg, xo = np.nanmax(bfx)-np.nanmin(bfx), np.nanmin(bfx)
-    yg, yo = np.nanmax(bfy)-np.nanmin(bfy), np.nanmin(bfy)
-    nbfx = (bfx-xo)/xg
-    nbfy = (bfy-yo)/yg
-    l = np.cumsum(np.sqrt(np.diff(nbfx)**2+np.diff(nbfy)**2))
-    nx, ny = [np.interp(np.linspace(l[0], l[-1], bins+1), l, (xy[1:]+xy[:-1])/2) for xy in [nbfx,nbfy]]
-    
-    
-    last = (-1,0)
-    
+        x_range, y_range = [v.max() - v.min() for v in [x, y]]
+        plt.ylim([y.min() - y_range * .1, y.max() + y_range * .1])
+        plt.xlim([x.min() - x_range * .1, x.max() + x_range * .1])
+
+    # divide curve into N segments of same normalized curve length
+    xg, xo = np.nanmax(bfx) - np.nanmin(bfx), np.nanmin(bfx)
+    yg, yo = np.nanmax(bfy) - np.nanmin(bfy), np.nanmin(bfy)
+    nbfx = (bfx - xo) / xg
+    nbfy = (bfy - yo) / yg
+    l = np.cumsum(np.sqrt(np.diff(nbfx)**2 + np.diff(nbfy)**2))
+    nx, ny = [np.interp(np.linspace(l[0], l[-1], bins + 1), l, (xy[1:] + xy[:-1]) / 2) for xy in [nbfx, nbfy]]
+
+    last = (-1, 0)
+
     pc = []
-    used = np.zeros_like(x).astype(np.bool)
-    for i in range(0,len(nx)):
-        i1,i2 = max(0,i-1), min(len(nx)-1,i+1)
-        a =-(nx[i2]-nx[i1])/ (ny[i2]-ny[i1])
-        b = (ny[i]-(a*nx[i]))*yg+yo
-        a *=yg/xg
+    used = np.zeros_like(x).astype(bool)
+    for i in range(0, len(nx)):
+        i1, i2 = max(0, i - 1), min(len(nx) - 1, i + 1)
+        a = -(nx[i2] - nx[i1]) / (ny[i2] - ny[i1])
+        b = (ny[i] - (a * nx[i])) * yg + yo
+        a *= yg / xg
         x_ = [np.nanmin(x), np.nanmax(x)]
-        m1 = np.sign(last[0])*y < np.sign(last[0])*((x-xo)*last[0]+last[1])
-        m2 = np.sign(a)*y>np.sign(a)*(a*(x-xo)+b)
-        m = m1&m2&~used
+        m1 = np.sign(last[0]) * y < np.sign(last[0]) * ((x - xo) * last[0] + last[1])
+        m2 = np.sign(a) * y > np.sign(a) * (a * (x - xo) + b)
+        m = m1 & m2 & ~used
         if plt:
-            plt.plot(x_, ((a)*(x_-xo))+b)
-            plt.plot(x[m], y[m],'.')
-        
-        if np.sum(m)>=bin_min_count:
+            plt.plot(x_, ((a) * (x_ - xo)) + b)
+            plt.plot(x[m], y[m], '.')
+
+        if np.sum(m) >= bin_min_count:
             pc.append((np.median(x[m]), np.median(y[m])))
-            used = used|m
-        last = (a,b)
+            used = used | m
+        last = (a, b)
     #bfx,bfy = zip(*pc)
     if plt:
         pbfx, pbfy = np.array(pc).T
-        plt.plot(bfx,bfy, 'orange', label='initial_fit')
+        plt.plot(bfx, bfy, 'orange', label='initial_fit')
         plt.plot(pbfx, pbfy, 'gray', label="perpendicular fit")
         plt.legend()
     #PlotData(None, bfx,bfy)
-    bin_x_fit, bin_y_fit = np.array(pc).T 
+    bin_x_fit, bin_y_fit = np.array(pc).T
     return bin_x_fit, _interpolate_fit(bin_x_fit, bin_y_fit, kind="linear")
 
-#Create mean function
+# Create mean function
+
+
 def _interpolate_fit(bin_x_fit, bin_y_fit, kind='linear'):
     def fit(x):
-        x = np.atleast_1d(x)[:].copy().astype(np.float)
-        x[x<bin_x_fit[0]] = np.nan
-        x[x>bin_x_fit[-1]] = np.nan
-        m = ~(np.isnan(bin_x_fit)|np.isnan(bin_y_fit))
+        x = np.atleast_1d(x)[:].copy().astype(float)
+        x[x < bin_x_fit[0]] = np.nan
+        x[x > bin_x_fit[-1]] = np.nan
+        m = ~(np.isnan(bin_x_fit) | np.isnan(bin_y_fit))
         return interp1d(bin_x_fit[m], bin_y_fit[m], kind)(x[:])
     return fit
-

wetb/signal/fit/_fourier_fit.py

@@ -13,13 +13,13 @@ def fourier_fit(y, max_nfft, x=None):
     return d, lambda deg : np.interp(deg%360, d, F2x(x2F(y, max_nfft, x)))
 
 def fourier_fit_old(y, nfft):
-    F = np.zeros(len(y), dtype=np.complex)
+    F = np.zeros(len(y), dtype=complex)
     F[:nfft + 1] = x2F(y, nfft)[:nfft + 1]
     return np.fft.ifft(F) * len(F)
 
 def F2x(F_coefficients):
     """Compute signal from Fourier coefficients"""
-    F = np.zeros(360, dtype=np.complex)
+    F = np.zeros(360, dtype=complex)
     nfft = len(F_coefficients) // 2
     F[:nfft + 1] = np.conj(F_coefficients[:nfft + 1])
     F[1:nfft + 1] += (F_coefficients[-nfft:][::-1])
@@ -53,7 +53,7 @@ def x2F(y, max_nfft, x=None):
         b[r] = 2 * np.nansum(y * np.sin(i * theta))
     AB = np.linalg.solve(a, b)
 
-    F = np.zeros(n, dtype=np.complex)
+    F = np.zeros(n, dtype=complex)
 
     F = np.r_[AB[0], (AB[1:nfft + 1] + 1j * AB[nfft + 1:]), np.zeros(nfft) ]
     return F
@@ -73,5 +73,5 @@ def rF2x(rF):
     """Convert single sided Fourier components, that satisfies x(t) = sum(X(cos(iw)+sin(iw)), i=0..N) to non-complex signal"""
     rF = np.conj(rF)
     rF[1:] /= 2
-    rF = np.r_[rF, np.zeros(181 - len(rF), dtype=np.complex)]
+    rF = np.r_[rF, np.zeros(181 - len(rF), dtype=complex)]
     return np.fft.irfft(rF) * 360

wetb/signal/fix/_rotor_position.py

@@ -7,9 +7,10 @@ import numpy as np
 from wetb.signal.fit._spline_fit import spline_fit
 from wetb.signal.filters._differentiation import differentiation
 
+
 def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt=None):
     """Rotor position fitted with spline
-    
+
     Parameters
     ----------
     rotor_position : array_like
@@ -25,61 +26,63 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt
         Note that a significant speed up is achievable by specifying the parameter
     plt : PyPlot or None
         If PyPlot a visual interpretation is plotted
-            
+
     Returns
     -------
     y : nd_array
         Fitted rotor position
     """
-    
+
     from wetb.signal.subset_mean import revolution_trigger
-    
+
     t = np.arange(len(rotor_position))
     indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=4)
-    
+
     rp = rotor_position[:].copy()
-    
+
     for i in indexes:
         rp[i:] += 360
-        
+
     if fix_dt is None:
         fix_dt = find_fix_dt(rotor_position, sample_frq, rotor_speed)
-        
-    N = int(np.round(fix_dt*sample_frq))
-    N2 = N//2
+
+    N = int(np.round(fix_dt * sample_frq))
+    N2 = N // 2
     if plt:
-        a = (rp.max()-rp.min())/t.max()
-        plt.plot(t/sample_frq,rp-t*a,label='Continus rotor position (detrended)')
-    
+        a = (rp.max() - rp.min()) / t.max()
+        plt.plot(t / sample_frq, rp - t * a, label='Continuous rotor position (detrended)')
+
     points = []
-    for j, i in enumerate(range(0,len(rp), N)):
-        #indexes for subsets for overlapping a and b polynomials
-        i1 = max(0,i-N2)
-        i2 = min(len(rp),i+N2)
-        
-        #fit a polynomial
-        if i1<len(rp):
-            poly_coef = np.polyfit(t[i1:i2]-t[i1], rp[i1:i2], 1)
-            points.append((t[i],np.poly1d(poly_coef)(t[i]-t[i1])))
+    for j, i in enumerate(range(0, len(rp), N)):
+        # indexes for subsets for overlapping a and b polynomials
+        i1 = max(0, i - N2)
+        i2 = min(len(rp), i + N2)
+
+        # fit a polynomial
+        if i1 < len(rp):
+            poly_coef = np.polyfit(t[i1:i2] - t[i1], rp[i1:i2], 1)
+            points.append((t[i], np.poly1d(poly_coef)(t[i] - t[i1])))
             if plt:
-                plt.plot(t[i1:i2]/sample_frq, np.poly1d(poly_coef)(t[i1:i2]-t[i1])- t[i1:i2]*a, 'mc'[j%2], label=('', "Line fit for points (detrended)")[j<2])
-                
-    x,y = np.array(points).T
-    
+                plt.plot(t[i1:i2] / sample_frq, np.poly1d(poly_coef)(t[i1:i2] - t[i1]) - t[i1:i2]
+                         * a, 'mc'[j % 2], label=('', "Line fit for points (detrended)")[j < 2])
+
+    x, y = np.array(points).T
+
     if plt:
-        plt.plot(x/sample_frq,y-x*a,'.', label='Fit points (detrended)')
-        plt.plot(t/sample_frq, spline_fit(x,y)(t)-t*a, label='Spline (detrended)')
+        plt.plot(x / sample_frq, y - x * a, '.', label='Fit points (detrended)')
+        plt.plot(t / sample_frq, spline_fit(x, y)(t) - t * a, label='Spline (detrended)')
         plt.legend()
         plt.show()
-    fit = spline_fit(x,y)(t)
-    fit[differentiation(fit, "left")<0]=np.nan
-    return fit%360
+    fit = spline_fit(x, y)(t)
+    fit[differentiation(fit, "left") < 0] = np.nan
+    return fit % 360
+
 
 def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
     """Find the optimal fix_dt parameter for fix_rotor_position (function above).
     Optimal is defined as the value that minimizes the sum of squared differences 
     between differentiated rotor position and rotor speed 
-    
+
     Parameters
     ----------
     rotor_position : array_like
@@ -90,35 +93,37 @@ def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
         Rotor speed [RPM]
     plt : pyplot or None
         If pyplot, a visual interpretation is plotted
-        
+
     Returns
     -------
     y : int
         Optimal value for the fix_dt parameter for fix_rotor_position
-    
+
     """
     from wetb.signal.filters import differentiation
+
     def err(i):
         drp = differentiation(fix_rotor_position(rotor_position, sample_frq, rotor_speed, i))
-        rpm_pos = drp%180 / 360 * sample_frq * 60
+        rpm_pos = drp % 180 / 360 * sample_frq * 60
         return np.sum((rpm_pos - rotor_speed)**2)
-    
+
     best = 27
-    for step in [9,3,1]:
-        
-        x_lst = np.arange(-2,3)*step + best
+    for step in [9, 3, 1]:
+
+        x_lst = np.arange(-2, 3) * step + best
         res = [err(x) for x in x_lst]
         if plt is not None:
-            plt.plot(x_lst, res,'.-')
+            plt.plot(x_lst, res, '.-')
         best = x_lst[np.argmin(res)]
     if plt is not None:
         plt.show()
-    
-    return best 
+
+    return best
+
 
 def check_rotor_position(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None):
     """Rotor position fitted with spline
-    
+
     Parameters
     ----------
     rotor_position : array_like
@@ -134,50 +139,53 @@ def check_rotor_position(rotor_position, sample_frq, rotor_speed, max_rev_diff=1
         Note that a significant speed up is achievable by specifying the parameter
     plt : PyPlot or None
         If PyPlot a visual interpretation is plotted
-            
+
     Returns
     -------
     y : nd_array
         Fitted rotor position
     """
-    
+
     from wetb.signal.subset_mean import revolution_trigger
-    
-    t = np.arange(len(rotor_position))/sample_frq
+
+    t = np.arange(len(rotor_position)) / sample_frq
     indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=max_rev_diff)
-    
-    rotor_position = rotor_position[:]%360
-        
+
+    rotor_position = rotor_position[:] % 360
+
     rp_fit = fix_rotor_position(rotor_position, sample_frq, rotor_speed, None)
-        
+
     if plt:
         #a = (rp.max()-rp.min())/t.max()
         #plt.plot(t/sample_frq,rp-t*a,label='Continus rotor position (detrended)')
         f, axarr = plt.subplots(2)
-        print (rp_fit)
-        axarr[0].plot(t, rotor_position)
-        axarr[0].plot(indexes/sample_frq, rotor_position[indexes],'.')
-        axarr[0].plot(t, rp_fit)
-        axarr[1].plot(t, rotor_speed)
-        axarr[1].plot(indexes/sample_frq,60/( differentiation(indexes)/sample_frq))
-        print (t[170:200])
+        print(rp_fit)
+        axarr[0].plot(t, rotor_position, label='rotor position')
+        axarr[0].plot(indexes / sample_frq, rotor_position[indexes], '.', label='Revolution trigger')
+        axarr[0].plot(t, rp_fit, label='Fitted rotor position')
+        axarr[0].legend()
+        axarr[1].plot(t, rotor_speed, label='Rotor speed')
+        axarr[1].plot(indexes / sample_frq, 60 / (differentiation(indexes) /
+                                                  sample_frq), label='Rotor speed (revolution average)')
+        print(t[170:200])
         drp = differentiation(rp_fit)
         #drp[(drp<0)]= 0
-        axarr[1].plot(t, drp%180 / 360 * sample_frq * 60, label='fix')
-        plt.legend()
-    
-    i1,i2 = indexes[0], indexes[-1]
-    print ("Rev from rotor position", np.sum(np.diff(rotor_position[i1:i2])%180)/360)
-    print ("Rev from rotor speed", np.mean(rotor_speed[i1:i2])*(i2-i1)/60/sample_frq)
-    print ("Rev from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360)
-    
-    print ("Mean RPM from rotor speed", np.mean(rotor_speed))
-    print ("Mean RPM from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360 / ((i2-i1)/60/sample_frq))
-    
-    spr1 = (np.diff(indexes)/sample_frq)
-    
+        axarr[1].plot(t, drp % 180 / 360 * sample_frq * 60, label='Fix(calc from fitted rotor position)')
+        axarr[1].legend()
+
+    i1, i2 = indexes[0], indexes[-1]
+    print("Rev from rotor position", np.sum(np.diff(rotor_position[i1:i2]) % 180) / 360)
+    print("Rev from rotor speed", np.mean(rotor_speed[i1:i2]) * (i2 - i1) / 60 / sample_frq)
+    print("Rev from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2]) % 180) / 360)
+
+    print("Mean RPM from rotor speed", np.mean(rotor_speed))
+    print("Mean RPM from fitted rotor position", np.sum(
+        np.diff(rp_fit[i1:i2]) % 180) / 360 / ((i2 - i1) / 60 / sample_frq))
+
+    spr1 = (np.diff(indexes) / sample_frq)
+
     #rs1 = 60/( np.diff(indexes)/sample_frq)
-    spr2 = np.array([60/rotor_speed[i1:i2].mean() for i1,i2 in zip(indexes[:-1], indexes[1:])])
-    
-    err = spr1-spr2
-    print (err.max())
\ No newline at end of file
+    spr2 = np.array([60 / rotor_speed[i1:i2].mean() for i1, i2 in zip(indexes[:-1], indexes[1:])])
+
+    err = spr1 - spr2
+    print(err.max())

wetb/signal/interpolation.py

@@ -43,8 +43,8 @@ def interpolate(x, xp, yp, max_xp_step=None, max_dydxp=None, cyclic_range=None,
     [359,nan,0,175,350]
     """
 
-    xp = np.array(xp, dtype=np.float)
-    yp = np.array(yp, dtype=np.float)
+    xp = np.array(xp, dtype=float)
+    yp = np.array(yp, dtype=float)
     assert xp.shape[0] == yp.shape[0], "xp and yp must have same length (%d!=%d)" % (xp.shape[0], yp.shape[0])
     non_nan = ~(np.isnan(xp) & np.isnan(yp))
     yp = yp[non_nan]

wetb/signal/subset_mean.py

@@ -5,14 +5,12 @@ Created on 24/06/2016
 '''
 
 import numpy as np
-import unittest
 from wetb.utils.geometry import rpm2rads
 from _collections import deque
-from tables.tests.test_index_backcompat import Indexes2_0TestCase
 from wetb.signal.filters._differentiation import differentiation
 
 
-def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1) :
+def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1):
     """Calculate the density normalized mean power, taking acceleration of the rotor into account
 
     Parameters
@@ -46,21 +44,24 @@ def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225,
         rs1 = rotor_speed[trigger_indexes[:-1]]
         rs2 = rotor_speed[trigger_indexes[1:] - 1]
     else:
-        rs = np.array([rotor_speed[max(i - rotor_speed_mean_samples, 0):i - 1 + rotor_speed_mean_samples].mean() for i in trigger_indexes])
+        rs = np.array([rotor_speed[max(i - rotor_speed_mean_samples, 0):i - 1 +
+                                   rotor_speed_mean_samples].mean() for i in trigger_indexes])
         rs1 = rs[:-1]
         rs2 = rs[1:]
 
-
-    power = np.array([np.nanmean(power[i1:i2], 0) for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
+    power = np.array([np.nanmean(power[i1:i2], 0)
+                      for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
     if isinstance(air_density, (int, float)):
         if air_density != 1.225:
             power = power / air_density * 1.225
     else:
-        air_density = np.array([np.nanmean(air_density[i1:i2], 0) for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
+        air_density = np.array([np.nanmean(air_density[i1:i2], 0)
+                                for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
         power = power / air_density * 1.225
     return power + 1 / 2 * I * (rs2 ** 2 - rs1 ** 2) / (time[trigger_indexes[1:] - 1] - time[trigger_indexes[:-1]])
 
-def power_mean_func_kW(I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1) :
+
+def power_mean_func_kW(I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1):
     """Return a power mean function [kW] used to Calculate the density normalized mean power, taking acceleration of the rotor into account
 
     Parameters
@@ -87,14 +88,14 @@ def power_mean_func_kW(I, rotor_speed, time, air_density=1.225, rotor_speed_mean
     wsp_mean, power_mean = subset_mean([wsp, power],trigger_indexes,mean_func={1:turbine_power_mean})
     """
     def mean_power(power, trigger_indexes):
-        return power_mean(power * 1000, trigger_indexes, I, rotor_speed, time , air_density, rotor_speed_mean_samples) / 1000
+        return power_mean(power * 1000, trigger_indexes, I, rotor_speed, time, air_density, rotor_speed_mean_samples) / 1000
     return mean_power
 
 
 def subset_mean(data, trigger_indexes, mean_func={}):
     if isinstance(data, list):
         data = np.array(data).T
-    if len(data.shape)==1:
+    if len(data.shape) == 1:
         no_sensors = 1
     else:
         no_sensors = data.shape[1]
@@ -103,7 +104,8 @@ def subset_mean(data, trigger_indexes, mean_func={}):
         steps = np.diff(triggers[:, 0])
         lengths = np.diff(triggers)[:, 0]
         if np.all(steps == steps[0]) and np.all(lengths == lengths[0]):
-            subset_mean = np.mean(np.r_[data.reshape(data.shape[0],no_sensors),np.empty((steps[0],no_sensors))+np.nan][triggers[0][0]:triggers.shape[0] * steps[0] + triggers[0][0]].reshape(triggers.shape[0], steps[0], no_sensors)[:, :lengths[0]], 1)
+            subset_mean = np.mean(np.r_[data.reshape(data.shape[0], no_sensors), np.empty((steps[0], no_sensors)) + np.nan][triggers[0]
+                                                                                                                            [0]:triggers.shape[0] * steps[0] + triggers[0][0]].reshape(triggers.shape[0], steps[0], no_sensors)[:, :lengths[0]], 1)
         else:
             subset_mean = np.array([np.mean(data[i1:i2], 0) for i1, i2 in trigger_indexes])
         for index, func in mean_func.items():
@@ -113,15 +115,17 @@ def subset_mean(data, trigger_indexes, mean_func={}):
         steps = np.diff(trigger_indexes)
 
         if np.all(steps == steps[0]):
-            #equal distance
-            subset_mean = np.mean(data[trigger_indexes[0]:trigger_indexes[-1]].reshape([ len(trigger_indexes) - 1, steps[0], data.shape[1]]), 1)
+            # equal distance
+            subset_mean = np.mean(data[trigger_indexes[0]:trigger_indexes[-1]
+                                       ].reshape([len(trigger_indexes) - 1, steps[0], data.shape[1]]), 1)
         else:
-            subset_mean = np.array([np.mean(data[i1:i2], 0) for i1, i2 in zip(trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
+            subset_mean = np.array([np.mean(data[i1:i2], 0) for i1, i2 in zip(
+                trigger_indexes[:-1].tolist(), trigger_indexes[1:].tolist())])
         for index, func in mean_func.items():
             att = data[:, index]
             subset_mean[:, index] = func(att, trigger_indexes)
-    if len(data.shape)==1 and len(subset_mean.shape)==2:
-        return subset_mean[:,0]
+    if len(data.shape) == 1 and len(subset_mean.shape) == 2:
+        return subset_mean[:, 0]
     else:
         return subset_mean
 
@@ -137,9 +141,10 @@ def cycle_trigger(values, trigger_value=None, step=1, ascending=True, tolerance=
     else:
         return np.where((values[1:] < trigger_value - tolerance) & (values[:-1] >= trigger_value + tolerance))[0][::step]
 
+
 def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None):
     """Returns one index per revolution (minimum rotor position)
-    
+
     Parameters
     ----------
     rotor_position : array_like
@@ -148,84 +153,75 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1,
         Sample frequency [Hz]
     rotor_speed : array_like
         Rotor speed [RPM]
-        
+
     Returns
     -------
     nd_array : Array of indexes
     """
     if isinstance(rotor_speed, (float, int)):
-        rotor_speed = np.ones_like(rotor_position)*rotor_speed
-    deg_per_sample = rotor_speed*360/60/sample_frq
-    thresshold = deg_per_sample.max()*3
-    drp = (np.diff(rotor_position)+thresshold)%360-thresshold
+        rotor_speed = np.ones_like(rotor_position) * rotor_speed
+    deg_per_sample = rotor_speed * 360 / 60 / sample_frq
+    thresshold = deg_per_sample.max() * 3
+    drp = (np.diff(rotor_position) + thresshold) % 360 - thresshold
     rp = rotor_position
 
-    
+    rp = np.array(rotor_position).copy() % 360
+    # filter degree increase > thresshold
+    rp[np.r_[False, np.diff(rp) > thresshold]] = 180
+
+    upper_indexes = np.where((rp[:-1] > (360 - thresshold)) & (rp[1:] < (360 - thresshold)))[0]
+    lower_indexes = np.where((rp[:-1] > thresshold) & (rp[1:] < thresshold))[0] + 1
 
-    
-    rp = np.array(rotor_position).copy()%360
-    #filter degree increase > thresshold
-    rp[np.r_[False, np.diff(rp)>thresshold]] = 180
-    
-    upper_indexes = np.where((rp[:-1]>(360-thresshold))&(rp[1:]<(360-thresshold)))[0]
-    lower_indexes = np.where((rp[:-1]>thresshold)&(rp[1:]<thresshold))[0] +1 
-    
     if plt:
         plt.plot(rp)
-        plt.plot(lower_indexes, rp[lower_indexes],'.')
-        plt.plot(upper_indexes, rp[upper_indexes],'.')
-        
+        plt.plot(lower_indexes, rp[lower_indexes], '.')
+        plt.plot(upper_indexes, rp[upper_indexes], '.')
 
     # Best lower is the first lower after upper
     best_lower = lower_indexes[np.searchsorted(lower_indexes, upper_indexes)]
     upper2lower = best_lower - upper_indexes
-    trigger_indexes = best_lower[upper2lower<upper2lower.mean()*2].tolist()
+    trigger_indexes = best_lower[upper2lower < upper2lower.mean() * 2].tolist()
+
+    if len(trigger_indexes) > 1:
 
-    if len(trigger_indexes)>1:
-        
-    
-        
         rpm_rs = np.array([rev.mean() for rev in np.split(rotor_speed, trigger_indexes)[1:-1]])
-        rpm_i = 1/np.diff(trigger_indexes)*60*sample_frq
-        spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
+        rpm_i = 1 / np.diff(trigger_indexes) * 60 * sample_frq
+        spr_rs = np.array([rev.mean() for rev in np.split(1 / rotor_speed * 60 * sample_frq, trigger_indexes)[1:-1]])
         spr_i = np.diff(trigger_indexes)
-        
 
-        while np.any(spr_rs>spr_i*1.9):
-            i = np.where(spr_rs>spr_i*1.9)[0][0]
-            if np.abs(spr_i[i-1] + spr_i[i] - spr_rs[i])< np.abs(spr_i[i] + spr_i[i+1] - spr_rs[i]):
+        while np.any(spr_rs > spr_i * 1.9):
+            i = np.where(spr_rs > spr_i * 1.9)[0][0]
+            if np.abs(spr_i[i - 1] + spr_i[i] - spr_rs[i]) < np.abs(spr_i[i] + spr_i[i + 1] - spr_rs[i]):
                 del trigger_indexes[i]
             else:
-                del trigger_indexes[i+1]
+                del trigger_indexes[i + 1]
             spr_i = np.diff(trigger_indexes)
-            spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
-        
-            
-        
+            spr_rs = np.array([rev.mean() for rev in np.split(
+                1 / rotor_speed * 60 * sample_frq, trigger_indexes)[1:-1]])
+
         # if a revolution is too long add trigger
-        if np.any(rpm_rs>rpm_i*2.1):
+        if np.any(rpm_rs > rpm_i * 2.1):
             # >1 missing triggers
             raise NotImplementedError
-        trigger_indexes.extend([np.mean(trigger_indexes[i:i+2]) for i in np.where(rpm_rs>rpm_i*1.9)[0]])
-        trigger_indexes = np.sort(trigger_indexes).astype(np.int)
-    
-    
-        i1,i2 = trigger_indexes[0], trigger_indexes[-1]
-        nround_rotor_speed = np.nansum(rotor_speed[i1:i2]/60/sample_frq)
-        #mod = [v for v in [5,10,30,60,90] if v>thresshold][0]    
-        nround_rotor_position = len(trigger_indexes)-1 #np.nansum(np.diff(rotor_position[i1:i2])%mod)/360
+        trigger_indexes.extend([np.mean(trigger_indexes[i:i + 2]) for i in np.where(rpm_rs > rpm_i * 1.9)[0]])
+        trigger_indexes = np.sort(trigger_indexes).astype(int)
+
+        i1, i2 = trigger_indexes[0], trigger_indexes[-1]
+        nround_rotor_speed = np.nansum(rotor_speed[i1:i2] / 60 / sample_frq)
+        #mod = [v for v in [5,10,30,60,90] if v>thresshold][0]
+        nround_rotor_position = len(trigger_indexes) - 1  # np.nansum(np.diff(rotor_position[i1:i2])%mod)/360
         if max_rev_diff is not None:
-            diff_pct = abs(nround_rotor_position-nround_rotor_speed)/nround_rotor_position*100
-            assert diff_pct<max_rev_diff, "No of revolution mismatch: rotor_position (%d), rotor_speed (%.1f), diff %.1f%%"%(nround_rotor_position, nround_rotor_speed, diff_pct)
-         
-    
+            diff_pct = abs(nround_rotor_position - nround_rotor_speed) / nround_rotor_position * 100
+            assert diff_pct < max_rev_diff, "No of revolution mismatch: rotor_position (%d), rotor_speed (%.1f), diff %.1f%%" % (
+                nround_rotor_position, nround_rotor_speed, diff_pct)
+
     return trigger_indexes
-    
+
 
 def revolution_trigger_old(values, rpm_dt=None, dmin=5, dmax=10, ):
     """Return indexes where values are > max(values)-dmin and decreases more than dmax
     If RPM and time step is provided, triggers steps < time of 1rpm is removed   
-    
+
     Parameters
     ---------
     values : array_like
@@ -237,40 +233,39 @@ def revolution_trigger_old(values, rpm_dt=None, dmin=5, dmax=10, ):
         Minimum normal position increase between samples
     dmax : float or int, optional
         Maximum normal position increase between samples 
-            
-    
+
+
     Returns
     -------
     trigger indexes 
         [i1,i2,...,in] if rpm_dt is not provided
         [(start1,stop1),(start2,stop2),...,(startn, stopn)] if rpm_dt is provided
     """
-    
+
     values = np.array(values)
-    indexes = np.where((np.diff(values)<-dmax)&(values[:-1]>values.max()-dmax))[0]
-    
+    indexes = np.where((np.diff(values) < -dmax) & (values[:-1] > values.max() - dmax))[0]
+
     if rpm_dt is None:
         return indexes
     else:
         index_pairs = []
         rpm, dt = rpm_dt
-        d_deg = rpm *360/60*dt
+        d_deg = rpm * 360 / 60 * dt
         cum_d_deg = np.cumsum(d_deg)
-        lbound, ubound = values.max()-dmax, values.max()+dmax
-        index_pairs = [(i1,i2) for i1,i2, deg in zip(indexes[:-1], indexes[1:], cum_d_deg[indexes[1:]-1]-cum_d_deg[indexes[:-1]]) 
-                       if deg > lbound and deg<ubound]
+        lbound, ubound = values.max() - dmax, values.max() + dmax
+        index_pairs = [(i1, i2) for i1, i2, deg in zip(indexes[:-1], indexes[1:], cum_d_deg[indexes[1:] - 1] - cum_d_deg[indexes[:-1]])
+                       if deg > lbound and deg < ubound]
         return index_pairs
-    
-    
+
+
 def time_trigger(time, step, start=None, stop=None):
     if start is None:
         start = time[0]
     decimals = int(np.ceil(np.log10(1 / np.nanmin(np.diff(time)))))
     time = np.round(time - start, decimals)
 
-
     steps = np.round(np.diff(time), decimals)
-    if np.sum(steps == steps[0])/len(time)>.99: #np.all(steps == steps[0]):
+    if np.sum(steps == steps[0]) / len(time) > .99:  # np.all(steps == steps[0]):
         # equal time step
         time = np.r_[time, time[-1] + max(set(steps), key=list(steps).count)]
     if stop is None:
@@ -286,10 +281,9 @@ def non_nan_index_trigger(sensor, step):
     i = 0
     nan_indexes = deque(np.where(np.isnan(sensor))[0].tolist() + [len(sensor)])
     while i + step <= sensor.shape[0]:
-        if i+step<=nan_indexes[0]:
-            trigger.append((i,i+step))
-            i+=step
+        if i + step <= nan_indexes[0]:
+            trigger.append((i, i + step))
+            i += step
         else:
-            i = nan_indexes.popleft()+1
+            i = nan_indexes.popleft() + 1
     return trigger
-

wetb/signal/tests/test_first_order.py

@@ -14,9 +14,8 @@ from wetb.signal.fit._fourier_fit import F2x, x2F, rx2F
 
 class Test_first_order_filters(unittest.TestCase):
 
-
     def test_low_pass(self):
-        a = np.random.randint(0,100,100).astype(np.float)
+        a = np.random.randint(0, 100, 100000).astype(float)
         b = first_order.low_pass(a, 1, 1)
         self.assertLess(b.std(), a.std())
         if 0:
@@ -24,30 +23,29 @@ class Test_first_order_filters(unittest.TestCase):
             plt.plot(a)
             plt.plot(b)
             plt.show()
-          
-          
+
     def test_low_pass2(self):
         t = np.linspace(0, 1.0, 2001)
         xlow = np.sin(2 * np.pi * 5 * t)
         xhigh = np.sin(2 * np.pi * 250 * t)
         x = xlow + xhigh
-                
+
         b, a = signal.butter(8, 0.125)
         y = signal.filtfilt(b, a, x, padlen=150)
-        self.assertAlmostEqual(np.abs(y - xlow).max(),0,4)
+        self.assertAlmostEqual(np.abs(y - xlow).max(), 0, 4)
         if 0:
             import matplotlib.pyplot as plt
             plt.plot(x)
             plt.plot(y)
             plt.show()
-        
+
 
 #     def test_low_pass3(self):
 #         t = np.linspace(0, 1.0, 2001)
 # #         xlow = np.sin(2 * np.pi * 5 * t)
 # #         xhigh = np.sin(2 * np.pi * 250 * t)
 # #         x = xlow + xhigh
-#         x = np.sum([np.sin(t*x*2*np.pi) for x in range(1,200)],0)                
+#         x = np.sum([np.sin(t*x*2*np.pi) for x in range(1,200)],0)
 #         cutoff = .2
 #         b, a = signal.butter(8,cutoff)
 #         w, h = signal.freqs(b, a)
@@ -55,19 +53,19 @@ class Test_first_order_filters(unittest.TestCase):
 #         F = rx2F(x, max_nfft=len(t))
 #         tF = np.linspace(0, len(F)/t[-1],len(F))
 #         Fy = rx2F(y, max_nfft=len(t))
-# 
+#
 #         if 1:
 #             import matplotlib.pyplot as plt
 # #             plt.plot(x)
 # #             plt.plot(y)
 # #             plt.show()
-# 
+#
 #             plt.plot(tF, np.abs(F))
 #             plt.plot(tF, np.abs(Fy))
 #             plt.xlim([0,260])
 #             plt.show()
 #             print (b,a)
-#             
+#
 #             plt.plot(w, 20 * np.log10(abs(h)))
 #             plt.xscale('log')
 #             plt.title('Butterworth filter frequency response')
@@ -77,7 +75,7 @@ class Test_first_order_filters(unittest.TestCase):
 #             plt.grid(which='both', axis='both')
 #             plt.axvline(cutoff, color='green') # cutoff frequency
 #             plt.show()
-#      
+#
 #     def test_low_pass2(self):
 #         t = np.arange(100000)/10000
 #         dt = t[1]-t[0]
@@ -89,11 +87,11 @@ class Test_first_order_filters(unittest.TestCase):
 #         print (sig.shape)
 #         F = rx2F(sig, max_nfft=len(t))
 #         tF = np.linspace(0, len(F)/t[-1],len(F))
-#           
-#           
+#
+#
 #         sig_lp1 = first_order.low_pass(sig, dt, 30*dt)
 #         F_lp1 = rx2F(sig_lp1, max_nfft=len(t))
-#    
+#
 #         bb, ba = signal.butter(1,5, 'low', analog=True)
 #         print (bb, ba)
 #         w1, h1 = signal.freqs(bb, ba)
@@ -104,24 +102,24 @@ class Test_first_order_filters(unittest.TestCase):
 # #         sig_lp3 = signal.lfilter(bb,ba, sig)
 # #         F_lp3 = rx2F(sig_lp3, max_nfft=len(t))
 # #         w2, h2 = signal.freqs(bb, ba)
-# # #          
+# # #
 # #         self.assertLess(b.std(), a.std())
 #         if 1:
 #             import matplotlib.pyplot as plt
 #             plt.plot(t, sig)
 # #             plt.plot(t, np.sum([1/x*np.sin(t*x*2*np.pi) for x in [1]],0))
-# #             
+# #
 #             plt.plot(t, sig_lp1)
 #             plt.plot(t, sig_lp2)
 #             #plt.plot(t, sig_lp3)
-#               
+#
 #             plt.show()
-#                
+#
 # #             plt.plot(tF, np.abs(F))
 # #             plt.plot(tF, np.abs(F_lp1))
 # # #             plt.plot(tF, np.abs(F_lp3))
 # # #             plt.plot(tF, np.abs(F_lp2))
-# #             plt.plot([0,1000],[0.708,.708])              
+# #             plt.plot([0,1000],[0.708,.708])
 # #             plt.xlim([0,205])
 # #             plt.ylim([0,1])
 # #             plt.show()
@@ -135,21 +133,20 @@ class Test_first_order_filters(unittest.TestCase):
 # #             plt.grid(which='both', axis='both')
 # #             plt.axvline(100, color='green') # cutoff frequency
 # #             plt.show()
- 
-          
- 
+
+
 # #     def test_low_pass2(self):
 # #         F = [0,1,0,0,0,0,0,0,0,0,0,0,0,0,.1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,.1]
 # #         F +=[0]*(len(F)-1)
 # #         x = F2x(F)
 # #         a = np.tile(x, 3)
 # #         print (x.shape)
-# # #         a = np.random.randint(0,100,100).astype(np.float)
+# # #         a = np.random.randint(0,100,100).astype(float)
 # #         b = first_order.low_pass(a, .1, 1)
 # #         bb, ba = signal.butter(10,100, 'low', analog=True)
 # #         #c = signal.lfilter(bb,ba, a)
 # #         w, h = signal.freqs(bb, ba)
-# #          
+# #
 # # #         self.assertLess(b.std(), a.std())
 # #         if 1:
 # #             import matplotlib.pyplot as plt
@@ -168,12 +165,11 @@ class Test_first_order_filters(unittest.TestCase):
 # #             plt.grid(which='both', axis='both')
 # #             plt.axvline(100, color='green') # cutoff frequency
 # #             plt.show()
-# #              
-# #             
-    
+# #
+# #
 
     def test_high_pass(self):
-        a = np.random.randint(0,100,100).astype(np.float)
+        a = np.random.randint(0, 100, 100).astype(float)
         b = first_order.high_pass(a, 1, 1)
         self.assertLess(b.mean(), a.mean())
         if 0:
@@ -181,6 +177,8 @@ class Test_first_order_filters(unittest.TestCase):
             plt.plot(a)
             plt.plot(b)
             plt.show()
+
+
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testName']
-    unittest.main()
\ No newline at end of file
+    unittest.main()

wetb/signal/tests/test_interpolation.py

@@ -27,10 +27,10 @@ class TestInterpolation(unittest.TestCase):
 
 
     def test_interpolate2(self):
-        x = np.arange(0, 100, 5, dtype=np.float)
-        xp = np.arange(0, 100, 10, dtype=np.float)
+        x = np.arange(0, 100, 5, dtype=float)
+        xp = np.arange(0, 100, 10, dtype=float)
         xp = np.r_[xp[:3], xp[5:]]
-        yp = np.arange(10, dtype=np.float)
+        yp = np.arange(10, dtype=float)
         yp[7:8] = np.nan
         yp = np.r_[yp[:3], yp[5:]]
 
@@ -39,10 +39,10 @@ class TestInterpolation(unittest.TestCase):
         #yp = [0.0, 1.0, 2.0, 5.0, 6.0, nan, 8.0, 9.0]
 
         y = signal.interpolate(x, xp, yp)
-        np.testing.assert_array_equal(y[~np.isnan(y)], [0., 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 8. , 8.5, 9. ])
+        np.testing.assert_array_equal(y[~np.isnan(y)], [0., 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6., 8. , 8.5, 9. ])
 
         y = signal.interpolate(x, xp, yp, 10)
-        np.testing.assert_array_equal(y[~np.isnan(y)], [ 0. , 0.5, 1. , 1.5, 2. , 5. , 5.5, 8. , 8.5, 9. ])
+        np.testing.assert_array_equal(y[~np.isnan(y)], [ 0. , 0.5, 1. , 1.5, 2. , 5. , 5.5, 6., 8. , 8.5, 9. ])
 
 
 

wetb/signal/tests/test_subset_mean.py

@@ -39,7 +39,7 @@ class TestSubsetMean(unittest.TestCase):
         self.assertEqual(p[1], data[2501:5000, 0].mean())
 
     def test_non_nan_index_trigger(self):
-        sensor = np.arange(18).astype(np.float)
+        sensor = np.arange(18).astype(float)
         sensor[[5, 11]] = np.nan
         triggers = non_nan_index_trigger(sensor, 3)
         for i1, i2 in triggers:
@@ -48,7 +48,7 @@ class TestSubsetMean(unittest.TestCase):
 
 
     def test_subset_mean_trigger_tuple(self):
-        sensor = np.arange(18).astype(np.float)
+        sensor = np.arange(18).astype(float)
         triggers = non_nan_index_trigger(sensor, 3)
         np.testing.assert_array_equal(subset_mean(sensor, triggers), [ 1., 4., 7., 10., 13., 16])
 
@@ -58,7 +58,7 @@ class TestSubsetMean(unittest.TestCase):
         np.testing.assert_array_equal(subset_mean(sensor, triggers), [ 2, 5, 8, 11, 14])
 
         #nan in the middle, noneq step and len
-        sensor = np.arange(18).astype(np.float)
+        sensor = np.arange(18).astype(float)
         sensor[[5, 11]] = np.nan
         triggers = non_nan_index_trigger(sensor, 3)
 

wetb/soil/KlinkvortModel.py

+# -*- coding: utf-8 -*-
+"""
+Created on Mon Jun 18 13:34:12 2018
+
+@author: shfe
+
+Note:
+This script is used for Klinkvort model used for cyclic accumulation calculation
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.close('all')
+
+class KlinkvortModel(object):
+    
+    
+    def __init__(self, xi_b, xi_c, N):
+        """
+        Initialize the input parameters
+        """
+        
+        self.xi_b = xi_b
+        self.xi_c = xi_c
+        self.N = N
+        
+    def alpha(self):
+        """
+        Calculate the parameter alpha
+        """
+        # Calculate the Tb, Tc
+        if not 0 <= self.xi_b <= 1:
+            raise ValueError('The maximum load ratio should between 0 to 1')
+        
+        if not -1 <= self.xi_c <= 1:
+            raise ValueError('The cyclic load ratio should between -1 to 1')
+            
+        if self.xi_b <= 0.022:
+            Tb = 0
+        else:
+            Tb = 0.61*self.xi_b - 0.013
+        
+        Tc = (self.xi_c + 0.63)*(self.xi_c - 1)*(self.xi_c - 1.64)
+        alpha = Tb * Tc
+        
+        return Tb, Tc, alpha 
+    
+    def evoluation(self, y1):
+        """
+        Evoluation model for cyclic accumulated rotation
+        """
+        
+        _, _, alpha = self.alpha()
+        yn = y1 * self.N**alpha
+        
+        return yn
+        
+    def equivalentnum(self, y1, yn, alpha):
+        """
+        Equivalent number for cylces
+        """
+        num = (yn/y1)**(1/alpha)
+        
+        return num
+        
+def superposition(xi_b, xi_c, N_cases, y1_cases, s0 = 0):
+
+    
+    # Check the array equal number of items
+    if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
+        raise ValueError('Size for input should be identical')
+    print(s0)
+    yn = np.zeros(np.size(xi_b))
+    for i, (b, c, N, y1) in enumerate(zip(xi_b, xi_c, N_cases, y1_cases)):
+        
+        if i==0:
+            
+            em = KlinkvortModel(b, c, N)
+            _, _, alpha = em.alpha()
+            if s0 == 0:
+                yn[i] = em.evoluation(y1)
+            else:
+                N_eq = (s0/y1)**(1/alpha)
+                if N_eq < 1e8:
+                    
+#                    N_eq = (s0/y1)**(1/alpha)
+                    yn[i] = y1*(N+N_eq)**alpha
+                else:
+                    yn[i] = s0
+                        
+                
+        
+        else:
+            em = KlinkvortModel(b, c, N)
+            _, _, alpha = em.alpha()
+            N_eq = (yn[i-1]/y1)**(1/alpha)
+            if N_eq < 1e8:
+#                N_eq = (yn[i-1]/y1)**(1/alpha)
+                yn[i] = y1*(N_eq + N)**alpha
+            else:
+                yn[i] = yn[i-1]
+
+    return yn
+    
+def alpha_super(xi_b, xi_c, N_cases):
+    # Check the alpha value
+    if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
+        raise ValueError('Size for input should be identical')
+    
+    alpha_list = np.zeros(np.size(xi_b))
+    for i, (b, c, N) in enumerate(zip(xi_b, xi_c, N_cases)):
+        
+        em = KlinkvortModel(b, c, N)
+        _, _, alpha = em.alpha()
+        alpha_list[i] = alpha
+        
+    return alpha_list
+    
+    
+    
+    
\ No newline at end of file

wetb/soil/LeblancModel.py

+# -*- coding: utf-8 -*-
+"""
+Created on Sun Mar 25 16:15:25 2018
+
+@author: shfe@dtu.dk
+
+The model proposed by Leblanc
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+plt.close('all')
+
+class LeblancModel(object):
+    
+    
+    def __init__(self, xi_b, xi_c, N):
+        """
+        Initialize the input parameters
+        """
+        
+        self.xi_b = xi_b
+        self.xi_c = xi_c
+        self.N = N
+        self.alpha = 0.31
+        
+    def coeff(self, rd = 0.38):
+        """
+        Calculate the parameter alpha
+        """
+        # Calculate the Tb, Tc
+        if not 0 <= self.xi_b <= 1:
+            raise ValueError('The maximum load ratio should between 0 to 1')
+        
+        if not -1 <= self.xi_c <= 1:
+            raise ValueError('The cyclic load ratio should between -1 to 1')
+        
+        if rd == 0.38:
+            if self.xi_b <= 0.051:
+                Tb = 0
+            else:
+                Tb = 0.4238*self.xi_b - 0.0217
+        
+        elif rd == 0.04:
+            if self.xi_b <= 0.1461:
+                Tb = 0
+            else:
+                Tb = 0.3087*self.xi_b - 0.0451
+                
+        else:
+            raise ValueError('The relative density value is not validated yet')
+                
+                
+        if -1 <= self.xi_c < -0.65:
+            Tc = 13.71*self.xi_c + 13.71
+            
+        elif -0.65 <= self.xi_c < 0:
+            Tc = -5.54*self.xi_c + 1.2
+            
+        else:
+            Tc = -1.2*self.xi_c + 1.2
+        
+        return Tb, Tc 
+    
+    def evoluation(self, y1):
+        """
+        Evoluation model for cyclic accumulated rotation
+        """
+        
+        Tb, Tc = self.coeff()
+        yn = y1 * Tb * Tc * self.N**0.31
+        
+        return yn
+        
+    def equivalentnum(self, y1, yn):
+        """
+        Equivalent number for cylces
+        """
+        Tb, Tc = self.coeff()
+        num = (yn/(y1*Tb*Tc))**(1/self.alpha)
+        
+        return num
+
+def superposition(xi_b, xi_c, N_cases, y1_cases, s0 = 0, y0 = 0):
+
+    
+    # Check the array equal number of items
+    if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
+        raise ValueError('Size for input should be identical')
+#    print(s0)
+    dyn = np.zeros(np.size(xi_b))
+    yn = np.zeros(np.size(xi_b))
+    y_static = np.zeros(np.size(xi_b))
+    for i, (b, c, N, y1) in enumerate(zip(xi_b, xi_c, N_cases, y1_cases)):
+        
+        if i==0:
+            
+            em = LeblancModel(b, c, N)
+            Tb, Tc = em.coeff()
+            if s0 == 0:
+                dyn[i] = em.evoluation(y1)
+            else:
+                N_eq = (s0/(y1*Tb*Tc))**(1/em.alpha)
+                if N_eq < 1e8:
+                    dyn[i] = (y1*Tb*Tc)*(N+N_eq)**em.alpha
+                else:
+                    dyn[i] = s0
+            
+            y_static[i] = max(y1, y0)
+                               
+        else:
+            em = LeblancModel(b, c, N)
+            Tb, Tc = em.coeff()
+            N_eq = (dyn[i-1]/(y1*Tb*Tc))**(1/em.alpha)
+            if N_eq < 1e8:
+                dyn[i] = (y1*Tb*Tc)*(N_eq + N)**em.alpha
+            else:
+                dyn[i] = dyn[i-1]
+                
+            y_static[i] = max(y_static[i-1], y1)
+        
+        yn[i] = dyn[i] + y_static[i]
+        
+    return dyn, y_static, yn
+    
+def coeff_super(xi_b, xi_c, N_cases):
+    # Check the alpha value
+    if not (np.size(xi_b) == np.size(xi_c) and np.size(xi_b) == np.size(N_cases)):
+        raise ValueError('Size for input should be identical')
+    
+    Tb_list = np.zeros(np.size(xi_b))
+    Tc_list = np.zeros(np.size(xi_c))
+    for i, (b, c, N) in enumerate(zip(xi_b, xi_c, N_cases)):
+        
+        em = LeblancModel(b, c, N)
+        Tb, Tc = em.coeff()
+        Tb_list[i] = Tb
+        Tc_list[i] = Tc
+        
+    return Tb_list, Tc_list
+    
+    
+def moment_ult(L, D, gamma):
+    """
+    Ultimate static moment
+    
+    Parameters
+    ----------------
+    L -- Pile length
+    D -- Pile diameter
+    gamma -- Submerged unit weight
+    
+    Output
+    ----------------
+    m_ult -- Ultimate static moment
+    """
+    
+    
+    m_ult = L**3*D*gamma*0.6
+    
+    return m_ult
+    
+def rot_moment(m, L, D, gamma, pa = 101, rd  = 0.04):
+    """
+    Ultimate static moment
+    
+    Parameters
+    ----------------
+    L -- Pile length
+    D -- Pile diameter
+    gamma -- Submerged unit weight
+    m -- Moment
+    
+    Output
+    ----------------
+    theta - Rotational angle in radians
+    """
+    
+    if rd == 0.04:
+        m_ult = L**3*D*gamma*0.6/1e3   # unit: MN
+        theta_norm = 0.038 * (m/m_ult)**2.33
+        theta = theta_norm*np.sqrt((L*gamma)/pa)
+    
+    else:
+        m_ult = L**3*D*gamma*1.2/1e3    # unit: MN
+        theta_norm = 0.042 * (m/m_ult)**1.92
+        theta = theta_norm*np.sqrt((L*gamma)/pa)
+
+    return theta
+

wetb/soil/hyperbolic_model.py

+# -*- coding: utf-8 -*-
+"""
+Created on Sun Apr 22 17:06:06 2018
+
+@author: shfe
+
+Hyperbolic model used for soil load displacement behaviour
+"""
+
+import numpy as np
+
+def hyperbolic_load(y, pu, a, b):
+    """
+    Static hyperbolic model developed by Lee
+    
+    Parameters
+    ----------------
+    y -- displacement
+    pu -- ultimate load
+    a, b -- coefficients
+    
+    Output
+    ----------------
+    p -- soil resistence pressure
+    """
+    
+    p = y/(a + b*y) * pu
+        
+    return p
+    
+def hyperbolic_disp(p, pu, a, b):
+    """
+    Static hyperbolic model developed by Lee
+    
+    Parameters
+    ----------------
+    p -- load
+    pu -- ultimate load
+    a, b -- coefficients
+    
+    Output
+    ----------------
+    y -- displacement
+    """
+    
+    y = p * a/(pu - p*b)
+        
+    return y
\ No newline at end of file

wetb/soil/py_model.py

+# -*- coding: utf-8 -*-
+"""
+Created on Sun Mar 25 16:15:25 2018
+
+@author: shfe@dtu.dk
+
+Static p-y soil model and cyclic degradation model
+"""
+
+import numpy as np
+
+def py_static(y, pu, yc):
+    """
+    Static p-y model developed by Matlock
+    
+    Parameters
+    ----------------
+    y -- displacement
+    pu -- ultimate stress
+    yc -- ultimate displacement
+    
+    Output
+    ----------------
+    p -- soil resistence pressure
+    """
+    
+    if y/yc <= 8:
+        p = 0.5*pu*(y/yc)**(1/3)
+        
+    else:
+        p = pu
+        
+    return p
+    
+def strain_static(p, pu, yc):
+    """
+    Static strain using p-y model developed by Matlock
+    
+    Parameters
+    ----------------
+    p -- stress
+    pu -- ultimate stress
+    yc -- ultimate displacement
+    
+    Output
+    ----------------
+    y -- soil strain
+    """
+    if p <= pu:
+        y = (p/(0.5*pu))**3*yc
+    
+    else:
+        raise ValueError("The applied stress is LARGER than ultimate strength")
+    
+    return y
+    
+def load_unload_py(p, pu, yc):
+    """
+    load/unload module using p-y model developed by Matlock
+    
+    Parameters
+    ----------------
+    p -- stress
+    pu -- ultimate stress
+    yc -- ultimate displacement
+    
+    Output
+    ----------------
+    y -- soil strain
+    y_cum -- cumulative soil displacement
+    """
+    
+    y = strain_static(p, pu, yc)
+    
+    if p <= 0.5*pu:  # elastic
+        y_cum = 0
+        
+    else:
+        y_cum = y - p/(0.5*pu/yc)
+        
+    return y, y_cum
+
+
+
+    
+def degra_coeff(N, y1, b):
+    """
+    Degradation factor for p-y model
+    
+    Parameters
+    ----------------
+    N -- number of cycle
+    y1 -- displacement predicted by static p-y curve
+    b -- pile diameter
+    
+    Output
+    ----------------
+    lambda -- Degradation factor
+    """
+    # check
+
+    lambda_N = np.abs(y1/(0.2*b)*np.log(N))
+    
+    try:
+        lambda_N <= 1
+    except:
+        print("soil must be degraded!!!")
+    
+    return lambda_N
+    
+def py_cyclic(y, pu, yc, N, y1, b):
+    """
+    Degradation p-y model with cyclic load
+    
+    Parameters
+    ----------------
+    y -- displacement
+    pu -- ultimate stress
+    yc -- ultimate displacement
+    N -- number of cycle
+    y1 -- displacement predicted by static p-y curve
+    b -- pile diameter
+    
+    Output
+    ----------------
+    p -- degraded soil resistence pressure
+    """
+    
+    lambda_N = degra_coeff(N, y1, b)
+    p = py_static(y, pu, yc) * (1 - lambda_N)
+
+    return p
+
+def py_parcel(load_parcel, cycle_parcel, pu, yc, b):
+    """
+    Degradation p-y model after a parcel of cyclic load
+    
+    Parameters
+    ----------------
+    load_parcel -- load parcel
+    cycle_parcel -- number of cycle of load parcel
+    yc -- ultimate displacement
+    b -- pile diameter
+    
+    Output
+    ----------------
+    p -- degraded soil resistence pressure
+    """
+    
+    pu_N = np.zeros(np.size(load_parcel)+1)
+    pu_N[0] = pu
+    lambda_N = np.zeros(np.size(load_parcel))
+    y_N = np.zeros(np.size(load_parcel))
+    
+    for i in np.arange(np.size(load_parcel)):
+        load_i = load_parcel[i]
+        cycle_i = cycle_parcel[i]
+        if np.isnan(load_i):
+            y_i = 0
+        else:
+            y_i = strain_static(load_i, pu_N[i], yc)
+        y_N[i] = y_i
+        lambda_i = degra_coeff(cycle_i, y_i, b)
+        lambda_N[i] = lambda_i
+        pu_N[i+1] = pu_N[i] * (1-lambda_i)
+    
+    return y_N, lambda_N, pu_N

wetb/standard_models/power_curve.py

+import numpy as np
+
+
+def standard_power_curve(power_norm, diameter, turbulence_intensity=.1, shear=(0, 100),
+                         rho=1.225, max_cp=.49,
+                         gear_loss_const=.01, gear_loss_var=.014, generator_loss=0.03, converter_loss=.03):
+    """Generate standard power curve
+
+    The method is extracted from an Excel spreadsheet made by Kenneth Thomsen, DTU Windenergy and
+    ported into python by Mads M. Pedersen, DTU Windenergy.
+
+    Parameters
+    ----------
+    power_norm : int or float
+        Nominal power [kW]
+    diameter : int or float
+        Rotor diameter [m]
+    turbulence_intensity : float
+        Turbulence intensity [%]
+    shear : (power shear coefficient, hub height)
+        Power shear arguments\n
+        - Power shear coeficient, alpha\n
+        - Hub height [m]
+    rho : float optional
+        Density of air [kg/m^3], defualt is 1.225
+    max_cp : float
+        Maximum power coefficient
+    gear_loss_const : float
+        Constant gear loss [%]
+    gear_loss_var : float
+        Variable gear loss [%]
+    generator_loss : float
+        Generator loss [%]
+    converter_loss : float
+
+    Examples
+    --------
+    wsp, power = standard_power_curve(10000, 178.3)
+    plot(wsp, power)
+    show()
+    """
+
+
+    area = (diameter / 2) ** 2 * np.pi
+    wsp_lst = np.arange(0.5, 25, .5)
+    sigma_lst = wsp_lst * turbulence_intensity
+
+    def norm_dist(x, my, sigma):
+        if turbulence_intensity > 0:
+            return 1 / np.sqrt(2 * np.pi * sigma ** 2) * np.exp(-(x - my) ** 2 / (2 * sigma ** 2))
+        else:
+            return x == my
+    p_aero = .5 * rho * area * wsp_lst ** 3 * max_cp / 1000
+
+    # calc power - gear, generator and conv loss
+#     gear_loss = gear_loss_const * power_norm + gear_loss_var * p_aero
+#     p_gear = p_aero - gear_loss
+#     p_gear[p_gear < 0] = 0
+#     gen_loss = generator_loss * p_gear
+#     p_gen = p_gear - gen_loss
+#     converter_loss = converter_loss * p_gen
+#     p_raw = p_gen - converter_loss
+#     p_raw[p_raw > power_norm] = power_norm
+    p_raw = p_aero - power_loss(p_aero, power_norm, gear_loss_const, gear_loss_var, generator_loss, converter_loss)
+
+    powers = []
+    shear_weighted_wsp = []
+    alpha, hub_height = shear
+    r = diameter / 2
+    z = np.linspace(hub_height - r, hub_height + r, 100, endpoint=True)
+    shear_factors = (z / hub_height) ** alpha
+    rotor_width = 2 * np.sqrt(r ** 2 - (hub_height - z) ** 2)
+
+    for wsp, sigma in zip(wsp_lst, sigma_lst):
+        shear_weighted_wsp.append(wsp * (np.trapz(shear_factors ** 3 * rotor_width, z) / (area)) ** (1 / 3))
+        ndist = norm_dist(wsp_lst, wsp, sigma)
+        powers.append((ndist * p_raw).sum() / ndist.sum())
+
+    return wsp_lst, lambda wsp : np.interp(wsp, wsp_lst, powers)
+    return wsp_lst, np.interp(shear_weighted_wsp, wsp_lst, powers)
+
+def power_loss(power_aero, power_norm, gear_loss_const=.01, gear_loss_var=.014, generator_loss=0.03, converter_loss=.03):
+    gear_loss = gear_loss_const * power_norm + gear_loss_var * power_aero
+    p_gear = power_aero - gear_loss
+    p_gear[p_gear < 0] = 0
+    gen_loss = generator_loss * p_gear
+    p_gen = p_gear - gen_loss
+    converter_loss = converter_loss * p_gen
+    p_electric = p_gen - converter_loss
+    p_electric[p_electric > power_norm] = power_norm
+    return power_aero - p_electric
+    
+
+if __name__ == '__main__':
+
+
+    from matplotlib.pyplot import plot, show
+    plot(*standard_power_curve(10000, 178.3, 0., (0, 119), gear_loss_const=.0, gear_loss_var=.0, generator_loss=0.0, converter_loss=.0))
+    plot(*standard_power_curve(10000, 178.3, 0.03, (0, 119), gear_loss_const=.0, gear_loss_var=.0, generator_loss=0.0, converter_loss=.0))
+
+
+
+
+    show()

wetb/utils/caching.py

@@ -3,15 +3,9 @@ Created on 07/02/2014
 
 @author: MMPE
 '''
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import division
-from __future__ import absolute_import
-from future import standard_library
 import sys
 from collections import OrderedDict
 import os
-standard_library.install_aliases()
 import inspect
 import numpy as np
 def set_cache_property(obj, name, get_func, set_func=None):

wetb/utils/cluster_tools/cluster_resource.py

@@ -34,15 +34,16 @@ def unix_path(path, cwd=None, fail_on_missing=False):
     if r:
         path = r[0]
     elif fail_on_missing:
-        raise FileExistsError("File or folder matching '%s' not found"%path)
+        raise FileExistsError("File or folder matching '%s' not found" % path)
     if cwd:
         path = os.path.relpath(path, cwd)
     if os.path.isdir(path):
-        path+="/"
-    return path.replace("\\","/")
+        path += "/"
+    return path.replace("\\", "/")
+
 
 def unix_path_old(filename):
-    filename = os.path.realpath(filename.replace("\\", "/")).replace("\\", "/")
+    filename = os.path.realpath(str(filename)).replace("\\", "/")
     ufn, rest = os.path.splitdrive(filename)
     ufn += "/"
     for f in rest[1:].split("/"):
@@ -51,22 +52,20 @@ def unix_path_old(filename):
             f_lst = [f_ for f_ in f_lst if f_ == f]
         elif len(f_lst) == 0:
             raise IOError("'%s' not found in '%s'" % (f, ufn))
-        else: # one match found
+        else:  # one match found
             ufn = os.path.join(ufn, f_lst[0])
     return ufn.replace("\\", "/")
 
 
-
-
 class Resource(object):
 
     def __init__(self, min_cpu, min_free):
         self.min_cpu = min_cpu
         self.min_free = min_free
-        self.cpu_free=0
+        self.cpu_free = 0
         self.acquired = 0
-        self.no_cpu="?"
-        self.used_by_user=0
+        self.no_cpu = "?"
+        self.used_by_user = 0
         self.resource_lock = threading.Lock()
 
     def ok2submit(self):
@@ -75,7 +74,7 @@ class Resource(object):
             #print ("ok2submit")
             total, free, user = self.check_resources()
             user = max(user, self.acquired)
-            user_max = max((self.min_cpu-user), self.cpu_free - self.min_free)
+            user_max = max((self.min_cpu - user), self.cpu_free - self.min_free)
             #print ("ok2submit", total, free, user, user_max)
             return user_max
         except:
@@ -95,7 +94,6 @@ class Resource(object):
         with self.resource_lock:
             self.acquired -= 1
 
-
     def update_resource_status(self):
         try:
             self.no_cpu, self.cpu_free, self.used_by_user = self.check_resources()
@@ -107,28 +105,27 @@ class SSHPBSClusterResource(Resource):
     finished = []
     loglines = {}
     is_executing = []
-    
+
     def __init__(self, sshclient, min_cpu, min_free):
         Resource.__init__(self, min_cpu, min_free)
         self.ssh = sshclient
         self.resource_lock = threading.Lock()
-        
+
     def glob(self, filepattern, cwd="", recursive=False):
         return self.ssh.glob(filepattern, cwd, recursive)
 
     @property
     def host(self):
         return self.ssh.host
-    
+
     @property
     def username(self):
         return self.ssh.username
 
-
     def new_ssh_connection(self):
         from wetb.utils.cluster_tools.ssh_client import SSHClient
         return SSHClient(self.host, self.ssh.username, self.ssh.password, self.ssh.port)
-        #return self.ssh
+        # return self.ssh
 
     def check_resources(self):
         with self.resource_lock:
@@ -153,22 +150,20 @@ class SSHPBSClusterResource(Resource):
             except Exception as e:
                 raise EnvironmentError(str(e))
 
-
     def jobids(self, jobname_prefix):
-            _, output, _ = self.ssh.execute('qstat -u %s' % self.username)
-            return [l.split()[0].split(".")[0] for l in output.split("\n")[5:] if l.strip() != "" and l.split()[3].startswith("h2l")]
+        _, output, _ = self.ssh.execute('qstat -u %s' % self.username)
+        return [l.split()[0].split(".")[0] for l in output.split("\n")[5:] if l.strip() != "" and l.split()[3].startswith("h2l")]
 
     def stop_pbsjobs(self, jobids):
         if not hasattr(jobids, "len"):
             jobids = list(jobids)
         self.ssh.execute("qdel %s" % (" ".join(jobids)))
-        
-    def setup_wine(self):    
+
+    def setup_wine(self):
         self.ssh.execute("""rm -f ./config-wine-hawc2.sh &&
 wget https://gitlab.windenergy.dtu.dk/toolbox/pbsutils/raw/master/config-wine-hawc2.sh &&
 chmod 777 config-wine-hawc2.sh &&
 ./config-wine-hawc2.sh""")
-   
 
 
 class LocalResource(Resource):
@@ -178,8 +173,6 @@ class LocalResource(Resource):
         #self.process_name = process_name
         self.host = 'Localhost'
 
-    
-    
     def check_resources(self):
         import psutil
         no_cpu = multiprocessing.cpu_count()

wetb/utils/cluster_tools/os_path.py

+
+import os
+
+
+def repl(path):
+    return path.replace("\\", "/")
+
+
+def abspath(path):
+    return repl(os.path.abspath(path))
+
+
+def relpath(path, start=None):
+    return repl(os.path.relpath(path, start))
+
+
+def realpath(path):
+    return repl(os.path.realpath(path))
+
+
+def pjoin(*path):
+    return repl(os.path.join(*path))
+
+
+def fixcase(path):
+    path = realpath(str(path)).replace("\\", "/")
+    p, rest = os.path.splitdrive(path)
+    p += "/"
+    for f in rest[1:].split("/"):
+        f_lst = [f_ for f_ in os.listdir(p) if f_.lower() == f.lower()]
+        if len(f_lst) > 1:
+            # use the case sensitive match
+            f_lst = [f_ for f_ in f_lst if f_ == f]
+        if len(f_lst) == 0:
+            raise IOError("'%s' not found in '%s'" % (f, p))
+        # Use matched folder
+        p = pjoin(p, f_lst[0])
+    return p
+
+
+def normpath(path):
+    return repl(os.path.normpath(path))
+
+
+def cluster_path(path):
+    try:
+        import win32wnet
+        drive, folder = os.path.splitdrive(abspath(path))
+        path = abspath(win32wnet.WNetGetUniversalName(drive, 1) + folder)
+    except Exception:
+        path = repl(path)
+    path = path.replace("//jess.dtu.dk", "/home")
+    path = path.replace("//mimer.risoe.dk/aiolos", "/mnt/aiolos")
+    path = path.replace("//mimer.risoe.dk", "/mnt/mimer")
+    return path

wetb/utils/cluster_tools/pbsfile.py

+import os
+import inspect
+from wetb.utils.cluster_tools.os_path import cluster_path, pjoin, normpath
+
+
+class Template():
+
+    def __init__(self, template):
+        self.template = template
+
+    def __call__(self, **kwargs):
+        s = self.template
+        found = True
+        while found:
+            found = False
+            for k, v in dict(kwargs).items():
+                if "[%s]" % k in s:
+                    found = True
+                s = s.replace("[%s]" % k, str(v))
+        return s
+
+
+pbs_template = Template('''### Jobid
+#PBS -N [jobname]
+### Standard Output
+#PBS -o [stdout_filename]
+[stderr]
+#PBS -W umask=0003
+### Maximum wallclock time format HOURS:MINUTES:SECONDS
+#PBS -l walltime=[walltime]
+#PBS -l nodes=[nodes]:ppn=[ppn]
+### Queue name
+#PBS -q [queue]
+cd "[workdir]"
+mkdir -p "stdout"
+if [ -z "$PBS_JOBID" ]; then echo "Run using qsub"; exit ; fi
+pwd
+[commands]
+exit
+''')
+
+
+class PBSFile():
+    _walltime = "00:30:00"
+
+    def __init__(self, workdir, jobname, commands, queue='workq', walltime='00:10:00', nodes=1, ppn=1, merge_std=True):
+        """Description
+
+        Parameters
+        ----------
+        walltime : int, str
+            wall time as string ("hh:mm:ss") or second (integer)
+
+        """
+        self.workdir = workdir
+        self.jobname = jobname
+        self.commands = commands
+        self.queue = queue
+        self.walltime = walltime
+        self.nodes = nodes
+        self.ppn = ppn
+        self.merge_std = merge_std
+        self.stdout_filename = normpath(pjoin(workdir, './stdout/%s.out' % jobname))
+        self.filename = "pbs_in/%s.in" % self.jobname
+
+    @property
+    def walltime(self):
+        return self._walltime
+
+    @walltime.setter
+    def walltime(self, walltime):
+        if isinstance(walltime, (float, int)):
+            from math import ceil
+            h = walltime // 3600
+            m = (walltime % 3600) // 60
+            s = ceil(walltime % 60)
+            self._walltime = "%02d:%02d:%02d" % (h, m, s)
+        else:
+            self._walltime = walltime
+
+    def __str__(self):
+        if self.merge_std:
+            stderr = "### merge stderr into stdout\n#PBS -j oe"
+        else:
+            stderr = '### Standard Error\n#PBS -e "./err/[jobname].err"'
+        if callable(self.commands):
+            commands = self.commands()
+        else:
+            commands = self.commands
+        return pbs_template(workdir=cluster_path(self.workdir),
+                            stdout_filename=cluster_path(self.stdout_filename),
+                            stderr=stderr,
+                            jobname=self.jobname,
+                            queue=self.queue,
+                            walltime=self.walltime,
+                            nodes=self.nodes,
+                            ppn=self.ppn,
+                            commands=commands)
+
+    def save(self, modelpath=None, filename=None):
+        modelpath = modelpath or self.workdir
+        self.filename = filename or self.filename
+        filename = os.path.join(modelpath, self.filename)
+        os.makedirs(os.path.dirname(filename), exist_ok=True)
+        with open(filename, 'w', newline='\n') as fid:
+            fid.write(str(self))
+        os.chmod(filename, 0o774)
+
+
+multirunner_template = Template("""echo "[make_dict]
+" | python
+
+for node in `cat $PBS_NODEFILE | sort | uniq`
+do
+
+     ssh -T $node << EOF &
+cd "[workdir]"
+python -c "[start_jobs]
+"
+EOF
+done
+wait
+""")
+
+
+class PBSMultiRunner(PBSFile):
+    def __init__(self, workdir, queue='workq', walltime='01:00:00', nodes=1, ppn=1, merge_std=True, pbsfiles=None,
+                 jobname='pbs_multirunner'):
+        if pbsfiles:
+            def fmt(pbsfile):
+                if isinstance(pbsfile, PBSFile):
+                    return pbsfile.filename
+                return pbsfile
+            self.PBS_FILES = "['" + "',\n             '".join([fmt(pbsfile) for pbsfile in pbsfiles]) + "']"
+        else:
+            self.PBS_FILES = """[os.path.join(root, f) for root, folders, f_lst in os.walk('.') for f in f_lst if f.endswith('.in')]"""
+
+        commands = multirunner_template(make_dict=self.get_src(self.make_dict),
+                                        start_jobs=self.get_src(self.start_jobs),
+                                        workdir=cluster_path(workdir)).replace("self.ppn", str(ppn))
+        PBSFile.__init__(self, workdir, jobname, commands, queue, walltime=walltime,
+                         nodes=nodes, ppn=ppn, merge_std=merge_std)
+        self.filename = "%s.%s" % (self.jobname, ("lst", "all")[pbsfiles is None])
+        self.pbsfiles = pbsfiles
+
+    def make_dict(self):
+        import os
+        import glob
+        import numpy as np
+        import re
+
+        # find available nodes
+        with open(os.environ['PBS_NODEFILE']) as fid:
+            nodes = set([f.strip() for f in fid.readlines() if f.strip() != ''])
+        pbs_files = eval(self.PBS_FILES)
+
+        # Make a list of [(pbs_in_filename, stdout_filename, walltime),...]
+        pat = re.compile(r'[\s\S]*#\s*PBS\s+-o\s+(.*)[\s\S]*(\d\d:\d\d:\d\d)[\s\S]*')
+
+        def get_info(f):
+            try:
+                with open(f) as fid:
+                    return (f,) + pat.match(fid.read()).groups()
+            except Exception:
+                return (f, f.replace('.in', '.out'), '00:30:00')
+        pbs_info_lst = map(get_info, pbs_files)
+
+        # sort wrt walltime
+        pbs_info_lst = sorted(pbs_info_lst, key=lambda fow: tuple(map(int, fow[2].split(':'))))[::-1]
+        # make dict {node1: pbs_info_lst1, ...} and save
+        d = dict([(f, pbs_info_lst[i::len(nodes)]) for i, f in enumerate(nodes)])
+        with open('pbs.dict', 'w') as fid:
+            fid.write(str(d))
+
+    def start_jobs(self):
+        import os
+        import multiprocessing
+        import platform
+        import time
+        with open('pbs.dict') as fid:
+            pbs_info_lst = eval(fid.read())[platform.node()]
+        arg_lst = ['echo starting %s && mkdir -p "%s" && env PBS_JOBID=$PBS_JOBID "%s" &> "%s" && echo finished %s' %
+                   (f, os.path.dirname(o), f, o, f) for f, o, _ in pbs_info_lst]
+        print(arg_lst[0])
+        print('Starting %d jobs on %s' % (len(arg_lst), platform.node()))
+        pool = multiprocessing.Pool(int('$PBS_NUM_PPN'))
+        res = pool.map_async(os.system, arg_lst)
+        t = time.time()
+        for (f, _, _), r in zip(pbs_info_lst, res.get()):
+            print('%-50s\t%s' % (f, ('Errorcode %d' % r, 'Done')[r == 0]))
+        print('Done %d jobs on %s in %ds' % (len(arg_lst), platform.node(), time.time() - t))
+
+    def get_src(self, func):
+        src_lines = inspect.getsource(func).split("\n")[1:]
+        indent = len(src_lines[0]) - len(src_lines[0].lstrip())
+        src = "\n".join([l[indent:] for l in src_lines])
+        if func == self.make_dict:
+            src = src.replace("eval(self.PBS_FILES)", self.PBS_FILES)
+        return src
+
+
+if __name__ == '__main__':
+    pbsmr = PBSMultiRunner(workdir="W:/simple1", queue='workq', nodes=2, ppn=10, pbsfiles=['/mnt/t1', "/mnt/t2"])
+    print(pbsmr.get_src(pbsmr.make_dict))
+    # pbsmr.save("c:/tmp/")

wetb/utils/cluster_tools/pbsfile_hawc2.py

+from wetb.hawc2.htc_file import HTCFile
+import os
+from os import path
+
+pbs_template = '''### Standard Output
+#PBS -N [jobname]
+#PBS -o [pbsoutdir]/[jobname].out
+### Standard Error
+#PBS -e [pbsoutdir]/[jobname].err
+#PBS -W umask=0003
+### Maximum wallclock time format HOURS:MINUTES:SECONDS
+#PBS -l walltime=[walltime]
+#PBS -l nodes=1:ppn=1
+### Queue name
+#PBS -q workq
+
+# ==============================================================================
+# single PBS mode: one case per PBS job
+# evaluates to true if LAUNCH_PBS_MODE is NOT set
+if [ -z ${LAUNCH_PBS_MODE+x} ] ; then
+  ### Create scratch directory and copy data to it
+  cd $PBS_O_WORKDIR
+  echo "current working dir (pwd):"
+  pwd
+  cp -R ./[modelzip] /scratch/$USER/$PBS_JOBID
+fi
+# ==============================================================================
+
+
+# ==============================================================================
+# single PBS mode: one case per PBS job
+# evaluates to true if LAUNCH_PBS_MODE is NOT set
+if [ -z ${LAUNCH_PBS_MODE+x} ] ; then
+  echo
+  echo 'Execute commands on scratch nodes'
+  cd /scratch/$USER/$PBS_JOBID
+  # create unique dir for each CPU
+  mkdir "1"; cd "1"
+  pwd
+  /usr/bin/unzip ../[modelzip]
+  mkdir -p [htcdir]
+  mkdir -p [resdir]
+  mkdir -p [logdir]
+  mkdir -p [turbdir]
+  cp -R $PBS_O_WORKDIR/[htcdir]/[jobname].htc ./[htcdir]
+  cp -R $PBS_O_WORKDIR/[turbdir][turbfileroot]*.bin [turbdir]
+  _HOSTNAME_=`hostname`
+  if [[ ${_HOSTNAME_:0:1} == "j" ]] ; then
+    WINEARCH=win64 WINEPREFIX=~/.wine winefix
+  fi
+# ==============================================================================
+
+# ------------------------------------------------------------------------------
+# find+xargs mode: 1 PBS job, multiple cases
+else
+  # with find+xargs we first browse to CPU folder
+  cd "$CPU_NR"
+fi
+# ------------------------------------------------------------------------------
+
+echo ""
+# ==============================================================================
+# single PBS mode: one case per PBS job
+# evaluates to true if LAUNCH_PBS_MODE is NOT set
+if [ -z ${LAUNCH_PBS_MODE+x} ] ; then
+  echo "execute HAWC2, fork to background"
+  time WINEARCH=win64 WINEPREFIX=~/.wine wine hawc2-latest ./[htcdir]/[jobname].htc &
+  wait
+# ==============================================================================
+
+# ------------------------------------------------------------------------------
+# find+xargs mode: 1 PBS job, multiple cases
+else
+  echo "execute HAWC2, do not fork and wait"
+  (time WINEARCH=win64 WINEPREFIX=~/.wine numactl --physcpubind=$CPU_NR wine hawc2-latest ./[htcdir]/[jobname].htc) 2>&1 | tee [pbsoutdir]/[jobname].err.out
+fi
+# ------------------------------------------------------------------------------
+
+
+### Epilogue
+# ==============================================================================
+# single PBS mode: one case per PBS job
+# evaluates to true if LAUNCH_PBS_MODE is NOT set
+if [ -z ${LAUNCH_PBS_MODE+x} ] ; then
+  ### wait for jobs to finish
+  wait
+  echo ""
+  echo "Copy back from scratch directory"
+  mkdir -p $PBS_O_WORKDIR/[resdir]
+  mkdir -p $PBS_O_WORKDIR/[logdir]
+  mkdir -p $PBS_O_WORKDIR/animation/
+  mkdir -p $PBS_O_WORKDIR/[turbdir]
+  cp -R [resdir]. $PBS_O_WORKDIR/[resdir].
+  cp -R [logdir]. $PBS_O_WORKDIR/[logdir].
+  cp -R animation/. $PBS_O_WORKDIR/animation/.
+
+  echo ""
+  echo "COPY BACK TURB IF APPLICABLE"
+  cd turb/
+  for i in `ls *.bin`; do  if [ -e $PBS_O_WORKDIR/[turbdir]$i ]; then echo "$i exists no copyback"; else echo "$i copyback"; cp $i $PBS_O_WORKDIR/[turbdir]; fi; done
+  cd /scratch/$USER/$PBS_JOBID/1/
+  echo "END COPY BACK TURB"
+  echo ""
+
+  echo "COPYBACK [copyback_files]/[copyback_frename]"
+  echo "END COPYBACK"
+  echo ""
+  echo ""
+  echo "following files are on node/cpu 1 (find .):"
+  find .
+# ==============================================================================
+# ------------------------------------------------------------------------------
+# find+xargs mode: 1 PBS job, multiple cases
+else
+  cd /scratch/$USER/$PBS_JOBID/$CPU_NR/
+  rsync -a --remove-source-files [resdir]. ../HAWC2SIM/[resdir].
+  rsync -a --remove-source-files [logdir]. ../HAWC2SIM/[logdir].
+  rsync -a --remove-source-files [pbsoutdir]/. ../HAWC2SIM/[pbsoutdir]/.
+  rsync -a --remove-source-files animation/. ../HAWC2SIM/animation/.
+
+  echo ""
+  echo "COPY BACK TURB IF APPLICABLE"
+  cd turb/
+  for i in `ls *.bin`; do  if [ -e $PBS_O_WORKDIR/[turbdir]$i ]; then echo "$i exists no copyback"; else echo "$i copyback"; cp $i $PBS_O_WORKDIR/[turbdir]; fi; done
+  cd /scratch/$USER/$PBS_JOBID/$CPU_NR/
+  echo "END COPY BACK TURB"
+  echo ""
+
+  echo "COPYBACK [copyback_files]/[copyback_frename]"
+  echo "END COPYBACK"
+  echo ""
+# ------------------------------------------------------------------------------
+fi
+exit
+'''
+
+def htc2pbs(htc_fn, walltime='00:40:00', zipfile=None):
+    """
+    Creates a PBS launch file (.p) based on a HAWC2 .htc file.
+    - Assumes htc files are within a htc/[casename]/ directory relative to current directory.
+    - Assumes there is a .zip file in the current directory which contains the turbine model.
+      If there is none, the zip file is set to 'model.zip' by default
+    - Will place a .p fine in pbs_in/[casename]/ directory relative to current directory.
+    -
+
+    Parameters
+    ----------
+    htc_fn : str
+        The file name and path to the .htc file relative to current directory.
+    walltime: str (default='00:40:00')
+        A string indicating the walltime of the job of the form 'HH:MM:SS'
+    zipfile: str (default=None)
+        The filename of the zipfile containing the wind turbine model files and
+        HAWC2 executable. if zipfile=None, searches the current directory for a
+        zip file. If none is found, sets zipfile to 'model.zip'
+
+
+    Returns
+    -------
+    str
+        The filename and path to the output .p file
+
+    Raises
+    ------
+    FileNotFoundError: If the file structure is not correct.
+    """
+
+
+
+    basename = path.relpath(path.dirname(htc_fn), 'htc')
+    jobname = path.splitext(path.basename(htc_fn))[0]
+    pbs_in_dir = path.join('pbs_in', basename)
+    if basename == '.':
+        raise FileNotFoundError('File structure is incorrect.')
+
+    if not zipfile:
+        try:
+            zipfile = [x for x in os.listdir() if x.lower().endswith('.zip')][0]
+        except:
+            print('No .zip file found in current directory. Set model zip to \'model.zip\'')
+            zipfile = 'model.zip'
+
+    #   get the required parameters for the pbs file from the htc file
+    htc = HTCFile(htc_fn) #modelpath='../..')
+    p = {
+        'walltime'      : walltime,
+        'modelzip'      : zipfile,
+        'jobname'       : jobname,
+        'htcdir'        : 'htc/' + basename,
+        'logdir'        :  path.dirname(htc.simulation.logfile.str_values())[2:] + '/',
+        'resdir'        : path.dirname(htc.output.filename.str_values())[2:] + '/',
+        'turbdir'       : path.dirname(htc.wind.mann.filename_u.str_values()) + '/',
+        'turbfileroot'  : path.basename(htc.wind.mann.filename_u.str_values()).split('u.')[0],
+        'pbsoutdir'     : 'pbs_out/' + basename
+        }
+
+
+    #Write pbs file based on template file and tags
+    if not os.path.exists(pbs_in_dir):
+        os.makedirs(pbs_in_dir)
+
+    template = str(pbs_template)
+
+    for key, value in p.items():
+        template = template.replace('[' + key + ']', value)
+
+    with open(os.path.join(pbs_in_dir, jobname + '.p'), 'w') as f:
+        f.write(template)

wetb/utils/cluster_tools/pbsjob.py

@@ -4,6 +4,9 @@ Created on 04/12/2015
 @author: mmpe
 '''
 import os
+import io
+from wetb.utils.cluster_tools.pbsfile import PBSFile
+from wetb.utils.cluster_tools.os_path import relpath
 
 NOT_SUBMITTED = "Job not submitted"
 PENDING = "Pending"
@@ -11,26 +14,35 @@ RUNNING = "Running"
 DONE = "Done"
 
 
+def pjoin(*args):
+    return os.path.join(*args).replace('\\', '/')
+
+
 class SSHPBSJob(object):
     _status = NOT_SUBMITTED
     nodeid = None
     jobid = None
 
-
     def __init__(self, sshClient):
         self.ssh = sshClient
 
-
-    def submit(self, job, cwd, pbs_out_file):
+    def submit(self, pbsfile, cwd=None, pbs_out_file=None):
         self.cwd = cwd
-        self.pbs_out_file = os.path.relpath(cwd + pbs_out_file).replace("\\", "/")
         self.nodeid = None
-        #self.execute()
-
+        if isinstance(pbsfile, PBSFile):
+            f = io.StringIO(str(pbsfile))
+            f.seek(0)
+            pbs_filename = pjoin(cwd, pbsfile.filename)
+            self.ssh.upload(f, pbs_filename)
+            self.pbs_out_file = pjoin(cwd, pbsfile.stdout_filename)
+            cwd = pbsfile.workdir
+            pbsfile = pbsfile.filename
+        else:
+            self.pbs_out_file = os.path.relpath(cwd + pbs_out_file).replace("\\", "/")
         cmds = ['rm -f %s' % self.pbs_out_file]
         if cwd != "":
             cmds.append("cd %s" % cwd)
-        cmds.append("qsub %s" % job)
+        cmds.append("qsub %s" % pbsfile)
         _, out, _ = self.ssh.execute(";".join(cmds))
         self.jobid = out.split(".")[0]
         self._status = PENDING
@@ -65,7 +77,6 @@ class SSHPBSJob(object):
                     return
                 raise e
 
-
     def is_executing(self):
         try:
             self.ssh.execute("qstat %s" % self.jobid)