Added fit + test to wetb.signal

wetb/signal/error_measures.py

+'''
+Created on 30/06/2016
+
+@author: MMPE
+'''
+
+from scipy.interpolate.interpolate import interp1d
+
+import numpy as np
+from wetb.signal.fit import bin_fit
+
+
+def rms(a, b):
+    """Calculate the Root-Mean-Squared Error of two value sets
+
+    Parameters
+    ---------
+    a : array_like
+        First value set
+    b : array_like
+        Second value set
+
+    Returns
+    -------
+    y : float
+        Root mean squared error of a and b
+
+    """
+    a, b = [np.array(ab[:]) for ab in [a, b]]
+    if a.shape != b.shape:
+        raise ValueError("Dimensions differ: %s!=%s" % (a.shape, b.shape))
+    if len(a) == 0:
+        return np.nan
+    return np.sqrt(np.nanmean((a - b) ** 2))
+
+
+def rms2fit(x, y, fit_func=bin_fit):
+    """
+    Calculate the rms error of the points (xi, yi) relative to the mean curve
+
+    The mean curve is computed by:\n
+    - Divide x into bins + 1 bins\n
+    - Remove bins with less than 2 elements\n
+    - Calculate the mean of x and y in the bins\n
+    - Do a linear interpolation between the bin mean values\n
+    - Extrapolate to the minimum and maximum value of x using the slope of the first and last line segment\n
+
+    Usefull for calculating e.g. power curve scatter
+
+    Parameters
+    ---------
+    x : array_like
+        x values
+    y : array_like
+        y values
+    bins : int or array_like, optional
+        If int: Number of control points for the mean curve, default is 10\n
+        If array_like: Bin egdes
+    kind : str or int, optional
+        Specifies the kind of interpolation as a string ('linear', 'nearest', 'zero', 'slinear',
+        'quadratic','cubic' where 'slinear', 'quadratic' and 'cubic' refer to a spline interpolation
+        of first, second or third order) or as an integer specifying the order of the spline
+        interpolator to use. Default is 'cubic'.
+    fit_func : function, optional
+        Function to apply on each bin to find control points for fit
+    
+    Returns
+    -------
+    err : float
+        Mean error of points compared to mean curve
+    f : function
+        Interpolation function
+    """
+    x, y = np.array(x[:]), np.array(y[:])
+    _, fit = fit_func(x,y)
+    return rms(fit(x),y), fit
+    
+def rms2fit_old(x, y, bins=10, kind='cubic', fit_func=np.nanmean, normalize_with_slope=False):
+    """
+    Calculate the rms error of the points (xi, yi) relative to the mean curve
+
+    The mean curve is computed by:\n
+    - Divide x into bins + 1 bins\n
+    - Remove bins with less than 2 elements\n
+    - Calculate the mean of x and y in the bins\n
+    - Do a linear interpolation between the bin mean values\n
+    - Extrapolate to the minimum and maximum value of x using the slope of the first and last line segment\n
+
+    Usefull for calculating e.g. power curve scatter
+
+    Parameters
+    ---------
+    x : array_like
+        x values
+    y : array_like
+        y values
+    bins : int or array_like, optional
+        If int: Number of control points for the mean curve, default is 10\n
+        If array_like: Bin egdes
+    kind : str or int, optional
+        Specifies the kind of interpolation as a string ('linear', 'nearest', 'zero', 'slinear',
+        'quadratic','cubic' where 'slinear', 'quadratic' and 'cubic' refer to a spline interpolation
+        of first, second or third order) or as an integer specifying the order of the spline
+        interpolator to use. Default is 'cubic'.
+    fit_func : function, optional
+        Function to apply on each bin to find control points for fit
+    normalize_with_slope : boolean, optional
+        If True, the mean error in each bin is normalized with the slope of the corresponding line segment
+
+    Returns
+    -------
+    err : float
+        Mean error of points compared to mean curve
+    f : function
+        Interpolation 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)
+
+    digitized = np.digitize(x, bins)
+    digitized[np.isnan(x) | np.isnan(y)] = -1
+
+    masks = [digitized == i for i in range(1, len(bins))]
+    import warnings
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        bin_x = np.array([np.nanmean(x[mask]) for mask in masks])
+        bin_y = np.array([fit_func(y[mask]) for mask in masks])
+        bin_count = np.array([np.sum(mask) for mask in masks])
+    bin_x_fit, bin_y = [b[bin_count >= 1] for b in [bin_x, bin_y]]
+
+    #extrapolate to first and last value of x
+    if bin_x_fit[0] > np.nanmin(x):
+        bin_y = np.r_[bin_y[0] - (bin_x_fit[0] - np.nanmin(x)) * (bin_y[1] - bin_y[0]) / (bin_x_fit[1] - bin_x_fit[0]), bin_y]
+        bin_x_fit = np.r_[np.nanmin(x), bin_x_fit]
+
+    if bin_x_fit[-1] < np.nanmax(x):
+        bin_y = np.r_[bin_y, bin_y[-1] + (np.nanmax(x) - bin_x_fit[-1]) * (bin_y[-1] - bin_y[-2]) / (bin_x_fit[-1] - bin_x_fit[-2]) ]
+        bin_x_fit = np.r_[bin_x_fit, np.nanmax(x)]
+
+    #Create mean function
+    f = lambda x : interp1d(bin_x_fit, bin_y, kind)(x[:])
+
+    #calculate error of segment
+    digitized = np.digitize(x, bin_x[bin_count > 0])
+    bin_err = np.array([rms(y[digitized == i], f(x[digitized == i])) for i in range(1, len(bin_x_fit))])
+    if normalize_with_slope:
+        slopes = np.diff(bin_y) / np.diff(bin_x_fit)
+        return np.nanmean(bin_err / np.abs(slopes)), f
+    return np.sqrt(np.nanmean(bin_err ** 2)), f
+
+
+def rms2mean(x, y, bins=10, kind='cubic', normalize_with_slope=False):
+    """
+    Calculate the rms error of the points (xi, yi) relative to the mean curve
+
+    The mean curve is computed by:\n
+    - Divide x into bins + 1 bins\n
+    - Remove bins with less than 2 elements\n
+    - Calculate the mean of x and y in the bins\n
+    - Do a linear interpolation between the bin mean values\n
+    - Extrapolate to the minimum and maximum value of x using the slope of the first and last line segment\n
+
+    Usefull for calculating e.g. power curve scatter
+
+    Parameters
+    ---------
+    x : array_like
+        x values
+    y : array_like
+        y values
+    bins : int or array_like, optional
+        If int: Number of control points for the mean curve, default is 10\n
+        If array_like: Bin egdes
+    kind : str or int, optional
+        Specifies the kind of interpolation as a string ('linear', 'nearest', 'zero', 'slinear',
+        'quadratic','cubic' where 'slinear', 'quadratic' and 'cubic' refer to a spline interpolation
+        of first, second or third order) or as an integer specifying the order of the spline
+        interpolator to use. Default is 'cubic'.
+    normalize_with_slope : boolean, optional
+        If True, the mean error in each bin is normalized with the slope of the corresponding line segment
+
+    Returns
+    -------
+    err : float
+        Mean error of points compared to mean curve
+    f : function
+        Interpolation function
+    """
+    return rms2fit(x, y, lambda x,y : bin_fit(x, y, bins, kind))
+
+
+def bootstrap_comparison(x, y, kind=1, N=15, M=100):
+    f_lst = []
+    y_lst = []
+    x_min, x_max = max(np.percentile(x, 2), np.sort(x)[2]), min(np.percentile(x, 98), np.sort(x)[-2])
+
+    y_arr = np.empty((M, N * 10)) + np.NaN
+    inside = 0
+    for i in range(M):
+        indexes = np.random.randint(0, len(x) - 1, len(x))
+        while x[indexes].min() > x_min or x[indexes].max() < x_max:
+            indexes = np.random.randint(0, len(x) - 1, len(x))
+        #indexes = np.arange(i, len(x), M)
+
+        _, f = rms2fit(x[indexes], y[indexes], lambda x,y : bin_fit(x,y, kind=kind, bins=N))
+        x_ = np.linspace(x_min, x_max, N * 10)
+        y_ = (f(x_))
+        if i > 10:
+            if np.all(y_ < np.nanmax(y_arr, 0)) and np.all(y_ > np.nanmin(y_arr, 0)):
+                inside += 1
+                if inside == 5:
+                    #print ("break", i)
+                    #break
+                    pass
+        y_arr[i, :] = y_
+
+    return (np.mean(np.std(y_arr, 0)), x_, y_arr[~np.isnan(y_arr[:, 0])])

wetb/signal/fit/__init__.py

+d = None
+d = dir()
+
+from wetb.signal.fit._linear_fit import *
+from wetb.signal.fit._bin_fit import *
+from wetb.signal.fit._fourier_fit import *
+
+__all__ = sorted([m for m in set(dir()) - set(d)])

wetb/signal/fit/_bin_fit.py

+
+import numpy as np
+from scipy.interpolate.interpolate import interp1d
+
+
+
+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
+        x observations
+    y : array_like
+        y observations
+    bins : int, array_like or (int, int)
+        if int: <bins> binx evenly distributed on the x-axis
+        if (xbins,ybins): <xbins> and <ybins> evenly distributed on the x and y axis respectively\n
+        Note that ybins only make sense if every y-value maps to a single x-value
+    kind : int or string
+        degree of polynomial for fit (argument passed to scipy.interpolate.interpolate.interp1d)
+    bin_func : function, optional
+        Statistic function to apply on bins, default is nanmean
+    bin_min_count : int, optional
+        Minimum number of observations in bins to include
+        Default is 3
+    lower_upper : str, int, (str,str), (int,int)
+        How to handle observations below and above first and last bin values. Can be:\n
+        - "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):
+        xbins, ybins = bins
+        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)
+            xbinsy = f1(x1)
+        else:
+            xbinsy = []
+        #x2, f2 = bin_fit(x,y, kind=1, bins=xbins)
+        bins = sorted(np.r_[xbinsx, xbinsy ])
+
+    digitized = np.digitize(x, bins)
+    digitized[np.isnan(x) | np.isnan(y)] = -1
+
+    masks = [digitized == i for i in range(1, len(bins))]
+    import warnings
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        bin_x = np.array([np.nanmean(x[mask]) for mask in masks])
+        bin_y = np.array([bin_func(y[mask]) for mask in masks])
+        bin_count = np.array([np.sum(mask) for mask in masks])
+    #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_y_fit = bin_y.copy()
+    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)
+    if bin_x_fit[0] > np.nanmin(x):
+        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":
+            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)
+    if bin_x_fit[-1] < np.nanmax(x):
+        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_x_fit = np.r_[bin_x_fit, np.nanmax(x)]
+        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)
+        
+    #Create mean function
+    def fit(x):
+        x = x[:].copy().astype(np.float)
+        x[x<bin_x_fit[0]] = np.nan
+        x[x>bin_x_fit[-1]] = np.nan
+        return interp1d(bin_x_fit, bin_y_fit, kind)(x[:])
+    
+    return bin_x_fit, fit
+    
+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
+        x observations
+    y : array_like
+        y observations
+    bins : int
+        Number of perpendicular bins
+    fit_func : function(x,y) -> (x,y) or None
+        Initial fit function
+        If None, bin_fit with same number of bins are used
+    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)
+    bfy = f(bfx)
+    
+    
+    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), 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
+        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
+        if plt:
+            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)
+    #bfx,bfy = zip(*pc)
+    if plt:
+        pbfx, pbfy = np.array(pc).T
+        plt.plot(bfx,bfy, 'orange', label='initial_fit')
+        plt.plot(pbfx, pbfy, 'gray', label="perpendicular fit")
+        plt.legend()
+    #PlotData(None, bfx,bfy)
+    
+    return np.array(pc).T
+

wetb/signal/fit/_fourier_fit.py

+'''
+Created on 07/07/2015
+
+@author: MMPE
+'''
+import numpy as np
+from wetb.signal.fit import bin_fit
+
+
+def fourier_fit(y, max_nfft, x=None):
+    """Approximate a signal, y, with Fourier fit"""
+    d = np.arange(360)
+    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[: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)
+    nfft = len(F_coefficients) // 2
+    F[:nfft + 1] = np.conj(F_coefficients[:nfft + 1])
+    F[1:nfft + 1] += (F_coefficients[-nfft:][::-1])
+    return np.real(np.fft.ifft(F) * len(F))
+
+def x2F(y, max_nfft, x=None):
+    """Compute Fourier coefficients from signal (signal may contain NANs)"""
+    d = np.arange(360)
+    if x is not None:
+        x,fit = bin_fit(x,y, d)
+        y = fit(d)
+        
+    nfft = min(max_nfft, len(y) // 2 + 1)
+    n = len(y)
+    N = nfft * 2 + 1
+    theta = np.linspace(0, 2 * np.pi, n + 1)[:n]
+    theta[np.isnan(y)] = np.nan
+    a = np.empty((nfft * 2 + 1, nfft * 2 + 1))
+    b = np.empty(nfft * 2 + 1)
+    A0_lst = lambda dF : 2 * np.nansum(1 * dF)
+    A_lst = lambda dF : [2 * np.nansum(np.cos(i * theta) * dF) for i in range(1, nfft + 1)]
+    B_lst = lambda dF : [2 * np.nansum(np.sin(i * theta) * dF) for i in range(1, nfft + 1)]
+    row = lambda dF : np.r_[A0_lst(dF), A_lst(dF), B_lst(dF)]
+
+
+    for i in range(nfft + 1):
+        a[i, :] = row(np.cos(i * theta))
+        b[i] = 2 * np.nansum(y * np.cos(i * theta))
+    for i, r in enumerate(range(nfft + 1, nfft * 2 + 1), 1):
+        a[r, :] = row(np.sin(i * theta))
+        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.r_[AB[0], (AB[1:nfft + 1] + 1j * AB[nfft + 1:]), np.zeros(nfft) ]
+    return F
+
+def rx2F(y, max_fft):
+    """Convert non-complex signal, y, to single sided Fourier components, that satifies x(t) = sum(X(cos(iw)+sin(iw)), i=0..N)"""
+    F = np.fft.rfft(y) / len(y)
+    F[1:-1] *= 2  # add negative side
+    F = np.conj(F)
+    return F[:max_fft + 1]
+
+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)]
+    return np.fft.irfft(rF) * 360

wetb/signal/fit/_linear_fit.py

+'''
+Created on 22. mar. 2017
+
+@author: mmpe
+'''
+import numpy as np
+def linear_fit(x,y):
+    from scipy.stats import linregress
+    slope, intercept, r_value, p_value, std_err = linregress(x,y)
+    return np.array([x.min(), x.max()]), lambda x : np.array(x)*slope+intercept , (slope, intercept)
\ No newline at end of file

wetb/signal/tests/test_error_measures.py

+'''
+Created on 20/07/2016
+
+@author: MMPE
+'''
+import os
+import unittest
+
+import numpy as np
+from wetb.signal.error_measures import rms2fit_old, rms2fit
+from wetb.signal.fit import bin_fit
+
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
+class Test(unittest.TestCase):
+
+
+#    def test_rms2mean(self):
+#        data = np.load(tfp + "wsp_power.npy")
+#        print (data.shape)
+#        wsp = data[:, 1].flatten()
+#        power = data[:, 0].flatten()
+#
+#        import matplotlib.pyplot as plt
+#        plt.plot(wsp, power, '.')
+#        x = np.linspace(wsp.min(), wsp.max(), 100)
+#        err, f = rms2mean(wsp, power)
+#        plt.plot(x, f(x), label='rms2mean, err=%.1f' % err)
+#        err, f = rms2fit(wsp, power, bins=20, kind=3, fit_func=np.median)
+#        plt.plot(x, f(x), label='rms2median, err=%.1f' % err)
+#        print (list(x))
+#        print (list(f(x)))
+#        plt.legend()
+#        plt.show()
+
+    def test_rms2fit(self):
+        x = np.array([10.234302313156817, 13.98517783627376, 7.7902362498947921, 11.08597865379001, 8.430623529700588, 12.279982848438033, 33.89151260027775, 12.095047111211629, 13.731371675689642, 14.858309846006723, 15.185588405617654])
+        y = np.array([28.515665187174477, 46.285328159179684, 17.763652093098958, 32.949007991536462, 20.788106673177083, 38.819226477864589, 96.53278479817709, 38.479684539388025, 46.072654127604167, 51.875484233398439, 53.379342967122398])
+        err, fit = rms2fit_old(x, y, kind=1, bins=15)
+        err2, fit2 = rms2fit(x, y, fit_func=lambda x,y: bin_fit(x,y, kind=1, bins=15, bin_min_count=1, lower_upper='extrapolate'))
+        self.assertAlmostEqual(err, 0.306,2)
+        if 0:
+            import matplotlib.pyplot as plt
+            plt.plot(x,y, '.')
+            x_ = np.linspace(x.min(), x.max(), 100)
+            plt.plot(x_, fit(x_), label='rms2fit, err=%.5f' % err)
+            plt.plot(x_, fit2(x_), label='rms2fit, err=%.5f' % err2)
+            plt.legend()
+            plt.show()
+        
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.test_rms2mean']
+    unittest.main()

wetb/signal/tests/test_fit.py

+'''
+Created on 17. okt. 2016
+
+@author: mmpe
+'''
+from wetb import gtsdf
+from wetb.signal.fit import perpendicular_bin_fit, linear_fit
+import numpy as np
+import os
+import unittest
+from wetb.signal.fit import fourier_fit
+from wetb.signal.error_measures import rms
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
+
+class TestFit(unittest.TestCase):
+    pass
+
+#     def testBinFit(self):
+#         ds = gtsdf.Dataset(tfp+"binfit.hdf5")
+#         import matplotlib.pyplot as plt
+#         print (ds.info['attribute_names'])
+#         x,y = ds('Wsp_metmast'), ds('Power')
+#         m = ~np.isnan(x)
+#         x,y = x[m],y[m]
+#         
+#         
+#         x_, f = bin_fit(x,y, bins=(10,10))
+#         rms2fit(x,y)
+#         
+#         plt.plot(x,y,'.')
+#         plt.plot(x_, f(x_),'.-r')
+#         plt.show()
+
+# 
+#     def testBinFit2(self):
+#         X = np.r_[np.arange(90), np.arange(90,100,.01)]
+#         Y = np.r_[np.arange(90)/10, np.arange(90,100,.01)*10-890]
+#         import matplotlib.pyplot as plt
+#         x,fit = bin_fit(X,Y, (10,0), 1)
+#         plt.plot(X,Y,'k')
+#         plt.plot(x,fit(x),'r.', label='(10,0)')
+#         x,fit = bin_fit(X,Y, (0,10), 1)
+#         plt.plot(x,fit(x),'g.', label='(0,10)')
+#         x,fit = bin_fit(X,Y, (10,10), 1)
+#         plt.plot(x,fit(x),'bx', label="(10,10)")
+#         plt.legend()
+#         plt.show()
+#          
+
+#     def testBinFit_non_monotomic_y(self):
+#         X = np.arange(0,10,.01)
+#         Y = np.sin(X)
+#         import matplotlib.pyplot as plt
+#         plt.plot(X,Y,'k')
+# #        x,fit = bin_fit(X,Y, (10,0), 1)
+# #        plt.plot(x,fit(x),'r.-', label='(10,0)')
+#         x,fit = bin_fit(X,Y, (0,10), 1)
+#         plt.plot(x,fit(x),'g.', label='(0,10)')
+# #         x,fit = bin_fit(X,Y, (10,10), 1)
+# #         plt.plot(x,fit(x),'bx', label="(10,10)")
+# #         plt.legend()
+#         plt.show()
+
+
+#     def testBinFit_missing_segments(self):
+#         ds = gtsdf.Dataset(tfp+"binfit.hdf5")
+#         import matplotlib.pyplot as plt
+#         print (ds.info['attribute_names'])
+#         x,y = ds('Wsp_metmast'), ds('Power')
+#         
+#         m = ~np.isnan(x)
+#         x,y = [v[(x[m]<8)|(x[m]>15)] for v in [x[m],y[m]]]
+#          
+#          
+#         x_, f = bin_fit(x,y, bins=10)
+#         x_ = np.linspace(x.min(), x.max(), 100 ) 
+#         
+#          
+#         plt.plot(x,y,'.')
+#         plt.plot(x_, f(x_),'.-r')
+#         plt.show()
+
+
+    def test_perpendicular_fit(self):
+        ds = gtsdf.Dataset(tfp+"binfit.hdf5")
+        
+        
+        x,y = ds('Wsp_metmast'), ds('Power')
+        if 0:
+            import matplotlib.pyplot as plt
+            fx, fy = perpendicular_bin_fit(x,y,30,plt=plt)
+            plt.show()
+
+        else:
+            fx, fy = perpendicular_bin_fit(x,y,30)
+        self.assertEqual(len(fx), 30)
+            
+            
+    def test_linear_fit(self):
+        x = np.random.rand(1000)*10
+        y = x*.5+3+(np.random.rand(len(x))-.5)
+        x_, fit, (a,b) = linear_fit(x,y)
+        self.assertAlmostEqual(a,.5,delta=.1)
+        self.assertAlmostEqual(b,3,delta=.1)
+        if 0:
+            import matplotlib.pyplot as plt
+            plt.plot(x,y,'.')
+            
+            plt.plot(x_, fit(x_), label='y=%fx+%f'%(a,b))
+            plt.legend()
+            plt.show()
+        
+        
+         
+    
+    def test_fourier_fit(self):
+        from numpy import nan
+        import matplotlib.pyplot as plt
+        y = [nan, nan, nan, nan, -0.36773350834846497, -0.34342807531356812, -0.3105124831199646, -0.2949407696723938, nan, nan, nan, nan, nan, nan, nan, nan, -0.37076538801193237, -0.35946175456047058, -0.35366204380989075, -0.34812772274017334, -0.32674536108970642, -0.31197881698608398, -0.31780806183815002, -0.31430944800376892, -0.32355087995529175, -0.35628914833068848, -0.39329639077186584, -0.46684062480926514, -0.48477476835250854, -0.50368523597717285, -0.51693356037139893, -0.50966787338256836, -0.49876394867897034, -0.486896812915802, -0.48280572891235352, -0.4708983302116394, -0.46562659740447998, -0.4582551121711731, -0.46219301223754883, -0.46569514274597168, -0.4741971492767334, -0.48431938886642456, -0.49597686529159546, -0.50340372323989868, -0.50065416097640991]
+        #y = [-0.36884582, -0.36256081, -0.35047901, -0.33841938, -0.3289246, -0.32291612, -0.32149044, -0.32851833, -0.34011644, -0.35467893, -0.36627313, -0.37245053, -0.37924927, -0.39883283, -0.38590872, -0.39833149, -0.40406495, -0.4102158, -0.41886991, -0.42862922, -0.43947089, -0.45299602, -0.46831554, -0.48249167, -0.49108803, -0.500368, -0.50779951, -0.51360059, -0.51370221, -0.50541216, -0.49272588, -0.47430229, -0.45657015, -0.44043627, -0.4286592, -0.41741648, -0.41344571, -0.40986174, -0.40896985, -0.40939313, -0.40635225, -0.40435526, -0.40015101, -0.39243227, -0.38454708]
+        
+        x = np.linspace(0, 360, len(y) + 1)[:len(y)]
+        #plt.plot(, y)
+        
+#         x_fit = fourier_fit.fit_old(y, 5)
+#         plt.plot(x, x_fit[::-1], label='fit2')
+        x_,fit = fourier_fit(y, 5)
+        self.assertAlmostEqual(rms(fit(x), y), 0.0056, 3)
+        if 0:
+            plt.plot(x, y, label='Observations')
+            plt.plot(x_,fit(x_), label='fit')
+            plt.legend()
+            plt.show()
+
+#     def test_fourier_fit2(self):
+#         import matplotlib.pyplot as plt
+#         t = np.linspace(0, 2 * np.pi, 9 + 1)[:-1]
+#         y = np.cos(t) + np.sin(t) + 1
+#         nfft = 4
+#         print (fourier_fit.x2F(y, nfft))
+#         plt.plot(deg(t), y, label='cos')
+#         plt.plot(fourier_fit.fit(y, nfft), label='fit')
+#         plt.plot(fourier_fit.F2x(np.fft.fft(y) / len(y)), label='fft')
+#         plt.legend()
+#         plt.show()
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()
\ No newline at end of file

wetb/signal/tests/test_fouier_fit.py

+'''
+Created on 07/07/2015
+
+@author: MMPE
+'''
+import unittest
+
+import numpy as np
+from wetb.signal.error_measures import rms
+from wetb.signal.fit._fourier_fit import fourier_fit_old, fourier_fit, x2F, rx2F
+from wetb.utils.geometry import deg
+
+
+class Test(unittest.TestCase):
+
+    def test_fourier_fit(self):
+        from numpy import nan
+        import matplotlib.pyplot as plt
+        y = [nan, nan, nan, nan, -0.36773350834846497, -0.34342807531356812, -0.3105124831199646, -0.2949407696723938, nan, nan, nan, nan, nan, nan, nan, nan, -0.37076538801193237, -0.35946175456047058, -0.35366204380989075, -0.34812772274017334, -0.32674536108970642, -0.31197881698608398, -0.31780806183815002, -0.31430944800376892, -0.32355087995529175, -0.35628914833068848, -0.39329639077186584, -0.46684062480926514, -0.48477476835250854, -0.50368523597717285, -0.51693356037139893, -0.50966787338256836, -0.49876394867897034, -0.486896812915802, -0.48280572891235352, -0.4708983302116394, -0.46562659740447998, -0.4582551121711731, -0.46219301223754883, -0.46569514274597168, -0.4741971492767334, -0.48431938886642456, -0.49597686529159546, -0.50340372323989868, -0.50065416097640991]
+        #y = [-0.36884582, -0.36256081, -0.35047901, -0.33841938, -0.3289246, -0.32291612, -0.32149044, -0.32851833, -0.34011644, -0.35467893, -0.36627313, -0.37245053, -0.37924927, -0.39883283, -0.38590872, -0.39833149, -0.40406495, -0.4102158, -0.41886991, -0.42862922, -0.43947089, -0.45299602, -0.46831554, -0.48249167, -0.49108803, -0.500368, -0.50779951, -0.51360059, -0.51370221, -0.50541216, -0.49272588, -0.47430229, -0.45657015, -0.44043627, -0.4286592, -0.41741648, -0.41344571, -0.40986174, -0.40896985, -0.40939313, -0.40635225, -0.40435526, -0.40015101, -0.39243227, -0.38454708]
+        
+        x = np.linspace(0, 360, len(y) + 1)[:len(y)]
+        #plt.plot(, y)
+        
+        x_fit = fourier_fit_old(y, 5)[::-1]
+        x_fit = np.r_[x_fit[-1],x_fit[:-1]]
+        
+        x_,fit = fourier_fit(y, 5)
+        self.assertAlmostEqual(rms(fit(x), y), 0.0056, 3)
+        if 0:
+            plt.plot(x, y, label='Observations')
+            plt.plot(x, x_fit, label='fit_old')
+            plt.plot(x_,fit(x_), label='fit')
+            plt.legend()
+            plt.show()
+            
+    def test_x2F(self):
+        t = np.linspace(0, 2 * np.pi, 9 + 1)[:-1]
+        y = 6 + 5*np.cos(t) + 4*np.sin(t) + 3*np.cos(2*t)+ 2*np.sin(2*t)
+        nfft = 4
+        np.testing.assert_array_almost_equal(x2F(y, nfft)[:nfft+1], rx2F(y, nfft))
+        np.testing.assert_array_almost_equal(rx2F(y, nfft),[6, 5+4j, 3+2j,0,0])
+        
+    def test_fourier_fit2(self):
+        import matplotlib.pyplot as plt
+        t = np.linspace(0, 2 * np.pi, 9 + 1)[:-1]
+        y = np.cos(t) + np.sin(t) + 1
+        nfft = 4
+        x,fit = fourier_fit(y, nfft)
+        if 0:
+            plt.figure()
+            plt.plot(deg(t), y, label='cos+sin+1')
+            plt.plot(x,fit(x), label='fit')
+            plt.legend()
+            plt.show()
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()
+

wetb/signal/tests/test_subset_mean.py

+'''
+Created on 18/07/2016
+
+@author: MMPE
+'''
+import os
+import unittest
+
+import numpy as np
+from wetb import gtsdf
+from wetb.signal.subset_mean import time_trigger, subset_mean, \
+    non_nan_index_trigger, revolution_trigger
+from wetb.utils.geometry import rpm2rads
+
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
+class TestSubsetMean(unittest.TestCase):
+    def test_time_trigger(self):
+        time = np.arange(0, 99.5, .5)
+        np.testing.assert_array_equal(time[time_trigger(time, 20)], [0, 20, 40, 60, 80])
+        np.testing.assert_array_equal(time[time_trigger(time + .5, 20)], [0, 20, 40, 60, 80])
+        np.testing.assert_array_equal(time[time_trigger(time + 100000000.5, 20)], [0, 20, 40, 60, 80])
+        np.testing.assert_array_equal(time[time_trigger(time, 20, 20, 60)], [20, 40, 60])
+        np.testing.assert_array_equal(time_trigger(np.arange(101), 20), [0, 20, 40, 60, 80, 100])
+        time, data, info = gtsdf.load(tfp + "subset_mean_test.hdf5")
+        np.testing.assert_array_equal(time_trigger(time, 200), [0, 5000, 10000, 15000])
+
+    def test_subset_mean(self):
+
+        time, data, info = gtsdf.load(tfp + "subset_mean_test.hdf5")
+        triggers = time_trigger(time, 100)
+        t, p = subset_mean([time, data[:, 0]], triggers).T
+        self.assertEqual(t[1], time[2500:5000].mean())
+        self.assertEqual(p[1], data[2500:5000, 0].mean())
+
+        triggers[1] = 2501
+        t, p = subset_mean([time, data[:, 0]], triggers).T
+        self.assertEqual(t[1], time[2501:5000].mean())
+        self.assertEqual(p[1], data[2501:5000, 0].mean())
+
+    def test_non_nan_index_trigger(self):
+        sensor = np.arange(18).astype(np.float)
+        sensor[[5, 11]] = np.nan
+        triggers = non_nan_index_trigger(sensor, 3)
+        for i1, i2 in triggers:
+            self.assertFalse(np.any(np.isnan(sensor[i1:i2])))
+        self.assertEqual(len(triggers), 4)
+
+
+    def test_subset_mean_trigger_tuple(self):
+        sensor = np.arange(18).astype(np.float)
+        triggers = non_nan_index_trigger(sensor, 3)
+        np.testing.assert_array_equal(subset_mean(sensor, triggers), [ 1., 4., 7., 10., 13., 16])
+
+        #start with nan eq step, eq len
+        sensor[0] = np.nan
+        triggers = non_nan_index_trigger(sensor, 3)
+        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[[5, 11]] = np.nan
+        triggers = non_nan_index_trigger(sensor, 3)
+
+        np.testing.assert_array_equal(subset_mean(sensor, triggers), [1, 7, 13, 16])
+
+    def test_cycle_trigger(self):
+        ds = gtsdf.Dataset(tfp+'azi.hdf5')
+        azi, rpm, time = [ds(x)[8403:8803] for x in ['azi','Rot_cor','Time']]
+        
+        trigger = revolution_trigger(azi)
+        np.testing.assert_array_equal(trigger, [ 17, 128, 241, 354])
+        azi[64] = 358
+        trigger = revolution_trigger(azi, (ds('Rot_cor'), np.diff(time).mean()))
+        
+#         import matplotlib.pyplot as plt
+#         t = np.arange(len(azi))
+#         plt.plot(t, azi)
+#         for i1,i2 in trigger:
+#             plt.plot(t[i1:i2],azi[i1:i2],'.--')
+#         plt.show()
+        np.testing.assert_array_equal(trigger, [ (128,241),(241,354)])
+
+
+if __name__ == "__main__":
+    unittest.main()