diff --git a/setup.py b/setup.py
index 8b4a3a679f9a4f930a02d198d08468dde76f21d9..368d369cda4e8a41f9348e793f5c8f2fdc436fbc 100644
--- a/setup.py
+++ b/setup.py
@@ -27,7 +27,7 @@ from Cython.Distutils import build_ext
def setup_package():
ex_info = [('wetb.fatigue_tools.rainflowcounting', ['pair_range', 'peak_trough', 'rainflowcount_astm']),
- ('wetb.signal_tools.filters', ['cy_filters'])]
+ ('wetb.signal.filters', ['cy_filters'])]
extlist = [Extension('%s.%s' % (module, n),
[os.path.join(module.replace(".","/"), n)+'.pyx'],
include_dirs=[np.get_include()]) for module, names in ex_info for n in names]
diff --git a/wetb/hawc2/htc_file.py b/wetb/hawc2/htc_file.py
index 0f559e2fb001c298c35e07d91c27087dd24c22ba..1c23189523a582029ac7679b47f75623a9d5e8f2 100644
--- a/wetb/hawc2/htc_file.py
+++ b/wetb/hawc2/htc_file.py
@@ -14,7 +14,7 @@ from io import open
from builtins import str
from future import standard_library
from wetb.utils.process_exec import pexec
-from wetb.utils.cluster_tools.cluster_resource import unix_path, unix_path_old
+from wetb.utils.cluster_tools.cluster_resource import unix_path_old
standard_library.install_aliases()
from collections import OrderedDict
diff --git a/wetb/hawc2/simulation.py b/wetb/hawc2/simulation.py
index 8920ecb1856a0a7ecb3bac3fd4b754cd8fa36972..387f2913c1dbbf81b514689451b182343303df4b 100755
--- a/wetb/hawc2/simulation.py
+++ b/wetb/hawc2/simulation.py
@@ -388,6 +388,8 @@ class Simulation(object):
def set_id(self, *args, **kwargs):
pass
+ def progress_callback(self,*args, **kwargs):
+ pass
class UpdateSimStatusThread(Thread):
def __init__(self, simulation, interval=1):
diff --git a/wetb/hawc2/simulation_resources.py b/wetb/hawc2/simulation_resources.py
index 85b4a497094f9c2f18b177452a5560878293ef6e..8904ebddaef98357858bc277ecd3eb3b49a3552c 100644
--- a/wetb/hawc2/simulation_resources.py
+++ b/wetb/hawc2/simulation_resources.py
@@ -16,12 +16,11 @@ import time
from wetb.hawc2 import log_file
from wetb.hawc2.log_file import LogInfo, LogFile
from wetb.hawc2.simulation import ERROR, ABORTED
-from wetb.utils.cluster_tools import pbsjob
+
from wetb.utils.cluster_tools.cluster_resource import LocalResource, \
- SSHPBSClusterResource, unix_path
+ SSHPBSClusterResource
+from wetb.utils.cluster_tools import pbsjob
from wetb.utils.cluster_tools.pbsjob import SSHPBSJob, NOT_SUBMITTED, DONE
-from wetb.utils.cluster_tools.ssh_client import SSHClient
-from wetb.utils.timing import print_time
from wetb.hawc2.htc_file import fmt_path
import numpy as np
@@ -249,7 +248,8 @@ class GormSimulationResource(PBSClusterSimulationResource):
def __init__(self, username, password, wine_cmd="WINEARCH=win32 WINEPREFIX=~/.wine32 wine"):
init_cmd = """export PATH=/home/python/miniconda3/bin:$PATH
source activate wetb_py3"""
- PBSClusterSimulationResource.__init__(self, "gorm.risoe.dk", username, password, 22, 25, 100, init_cmd, wine_cmd, "python")
+ from wetb.utils.cluster_tools.ssh_client import SSHClient
+ PBSClusterSimulationResource.__init__(self, SSHClient('gorm.risoe.dk', username, password, 22), 25, 100, init_cmd, wine_cmd, "python")
class PBSClusterSimulationHost(SimulationHost):
diff --git a/wetb/signal/error_measures.py b/wetb/signal/error_measures.py
new file mode 100644
index 0000000000000000000000000000000000000000..6fbcdddc2bfff3458016b04631992a14cb1f551f
--- /dev/null
+++ b/wetb/signal/error_measures.py
@@ -0,0 +1,219 @@
+'''
+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])])
diff --git a/wetb/signal_tools/filters/__init__.py b/wetb/signal/filters/__init__.py
similarity index 100%
rename from wetb/signal_tools/filters/__init__.py
rename to wetb/signal/filters/__init__.py
diff --git a/wetb/signal_tools/filters/cy_filters.py b/wetb/signal/filters/cy_filters.py
similarity index 100%
rename from wetb/signal_tools/filters/cy_filters.py
rename to wetb/signal/filters/cy_filters.py
diff --git a/wetb/signal_tools/filters/cy_filters.pyx b/wetb/signal/filters/cy_filters.pyx
similarity index 100%
rename from wetb/signal_tools/filters/cy_filters.pyx
rename to wetb/signal/filters/cy_filters.pyx
diff --git a/wetb/signal_tools/filters/despike.py b/wetb/signal/filters/despike.py
similarity index 95%
rename from wetb/signal_tools/filters/despike.py
rename to wetb/signal/filters/despike.py
index e65a61c9652bea27a0809cceba3302851558e94f..f7256c68b1bd499faf78987afa83b644ff8f4dee 100644
--- a/wetb/signal_tools/filters/despike.py
+++ b/wetb/signal/filters/despike.py
@@ -4,8 +4,8 @@ Created on 13/07/2016
@author: MMPE
'''
import numpy as np
-from wetb.signal_tools.filters.first_order import low_pass
-from wetb.signal_tools.filters import replacer
+from wetb.signal.filters.first_order import low_pass
+from wetb.signal.filters import replacer
replace_by_nan = replacer.replace_by_nan
diff --git a/wetb/signal_tools/filters/first_order.py b/wetb/signal/filters/first_order.py
similarity index 90%
rename from wetb/signal_tools/filters/first_order.py
rename to wetb/signal/filters/first_order.py
index 6bcfd04655633e423127d80a4537654aad9385cc..3b5824bd3117e24b4180a075e6d8195f5dcda68c 100644
--- a/wetb/signal_tools/filters/first_order.py
+++ b/wetb/signal/filters/first_order.py
@@ -4,7 +4,7 @@ Created on 10/01/2015
@author: mmpe
'''
import numpy as np
-from wetb.signal_tools.filters import cy_filters
+from wetb.signal.filters import cy_filters
def low_pass(input, delta_t, tau, method=1):
if isinstance(tau, (int, float)):
diff --git a/wetb/signal_tools/filters/replacer.py b/wetb/signal/filters/replacer.py
similarity index 100%
rename from wetb/signal_tools/filters/replacer.py
rename to wetb/signal/filters/replacer.py
diff --git a/wetb/signal/fit/__init__.py b/wetb/signal/fit/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..2578ab2dde2c387e6801224980c9c13848c1b548
--- /dev/null
+++ b/wetb/signal/fit/__init__.py
@@ -0,0 +1,8 @@
+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)])
diff --git a/wetb/signal/fit/_bin_fit.py b/wetb/signal/fit/_bin_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..f1124ff6fb0919a22f2ef8a1e45ad14a597ce9eb
--- /dev/null
+++ b/wetb/signal/fit/_bin_fit.py
@@ -0,0 +1,195 @@
+
+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
+
diff --git a/wetb/signal/fit/_fourier_fit.py b/wetb/signal/fit/_fourier_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..ae2e9b4ca21ef21eb92705a48ab6d1d19ee24e95
--- /dev/null
+++ b/wetb/signal/fit/_fourier_fit.py
@@ -0,0 +1,77 @@
+'''
+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_nfft, x=None):
+ """Convert non-complex signal, y, to single sided Fourier components, that satifies x(t) = sum(X(cos(iw)+sin(iw)), i=0..N)"""
+ d = np.arange(360)
+ if x is not None:
+ x,fit = bin_fit(x,y, d)
+ y = fit(d)
+ F = np.fft.rfft(y) / len(y)
+ F[1:-1] *= 2 # add negative side
+ F = np.conj(F)
+ return F[:max_nfft + 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
diff --git a/wetb/signal/fit/_linear_fit.py b/wetb/signal/fit/_linear_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..fb6d15813f101d3015b52cb53ba2a7f2e96dff48
--- /dev/null
+++ b/wetb/signal/fit/_linear_fit.py
@@ -0,0 +1,10 @@
+'''
+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
diff --git a/wetb/signal_tools/spectrum.py b/wetb/signal/spectrum.py
similarity index 100%
rename from wetb/signal_tools/spectrum.py
rename to wetb/signal/spectrum.py
diff --git a/wetb/utils/subset_mean.py b/wetb/signal/subset_mean.py
similarity index 100%
rename from wetb/utils/subset_mean.py
rename to wetb/signal/subset_mean.py
diff --git a/wetb/signal/tests/test_error_measures.py b/wetb/signal/tests/test_error_measures.py
new file mode 100644
index 0000000000000000000000000000000000000000..b83723250f6968d0d5d5400c6fb77167483ccf54
--- /dev/null
+++ b/wetb/signal/tests/test_error_measures.py
@@ -0,0 +1,54 @@
+'''
+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()
diff --git a/wetb/signal/tests/test_files/azi.hdf5 b/wetb/signal/tests/test_files/azi.hdf5
new file mode 100644
index 0000000000000000000000000000000000000000..a40195b540e2dc316a007861f7af3331112fbaf4
Binary files /dev/null and b/wetb/signal/tests/test_files/azi.hdf5 differ
diff --git a/wetb/signal/tests/test_files/binfit.hdf5 b/wetb/signal/tests/test_files/binfit.hdf5
new file mode 100644
index 0000000000000000000000000000000000000000..4e98d115f918c7324cec500cd4c590b703836516
Binary files /dev/null and b/wetb/signal/tests/test_files/binfit.hdf5 differ
diff --git a/wetb/signal/tests/test_files/subset_mean_test.hdf5 b/wetb/signal/tests/test_files/subset_mean_test.hdf5
new file mode 100644
index 0000000000000000000000000000000000000000..3a6f83e349f97c11652721bcb467551cc6bd8653
Binary files /dev/null and b/wetb/signal/tests/test_files/subset_mean_test.hdf5 differ
diff --git a/wetb/signal/tests/test_first_order.py b/wetb/signal/tests/test_first_order.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c138c4ed9a3f36c691f9a90b1954308eb31fa92
--- /dev/null
+++ b/wetb/signal/tests/test_first_order.py
@@ -0,0 +1,186 @@
+'''
+Created on 13. jan. 2017
+
+@author: mmpe
+'''
+import unittest
+
+from scipy import signal
+
+import numpy as np
+from wetb.signal.filters import first_order
+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)
+ b = first_order.low_pass(a, 1, 1)
+ self.assertLess(b.std(), a.std())
+ if 0:
+ import matplotlib.pyplot as plt
+ 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)
+ 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)
+# cutoff = .2
+# b, a = signal.butter(8,cutoff)
+# w, h = signal.freqs(b, a)
+# y = signal.filtfilt(b, a, x)
+# 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')
+# plt.xlabel('Frequency [radians / second]')
+# plt.ylabel('Amplitude [dB]')
+# plt.margins(0, 0.1)
+# 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]
+# hz2 = np.sin(t*2*2*np.pi) # frq = 2 hz or 2*2pi rad/s
+# hz10 = np.sin(t*10*2*np.pi) # frq = 10 hz or 10*2pi rad/s
+# hz20 = np.sin(t*100*2*np.pi) # frq = 100 hz or 100*2pi rad/s
+# #sig = np.sum([np.sin(t*x*2*np.pi) for x in range(1,200)],0)
+# sig = np.sum([np.sin(t*x*2*np.pi) for x in [1,5,10]],0)
+# 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)
+# sig_lp2 = signal.filtfilt(bb,ba, sig)
+# # print (sig_lp2)
+# # F_lp2 = rx2F(sig_lp2, max_nfft=len(t))
+# # bb, ba = signal.butter(10,10, 'low', analog=True)
+# # 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.xlim([0,205])
+# # plt.ylim([0,1])
+# # plt.show()
+# # plt.plot(w1, 20 * np.log10(abs(h1)))
+# # plt.plot(w2, 20 * np.log10(abs(h2)))
+# # plt.xscale('log')
+# # plt.title('Butterworth filter frequency response')
+# # plt.xlabel('Frequency [radians / second]')
+# # plt.ylabel('Amplitude [dB]')
+# # plt.margins(0, 0.1)
+# # 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)
+# # 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
+# # plt.plot(a)
+# # #plt.ylim([-2,2])
+# # plt.plot(b)
+# # #plt.plot(c)
+# # plt.show()
+# # print (h)
+# # plt.plot(w, abs(h))
+# # plt.xscale('log')
+# # plt.title('Butterworth filter frequency response')
+# # plt.xlabel('Frequency [radians / second]')
+# # plt.ylabel('Amplitude [dB]')
+# # plt.margins(0, 0.1)
+# # 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)
+ b = first_order.high_pass(a, 1, 1)
+ self.assertLess(b.mean(), a.mean())
+ if 0:
+ import matplotlib.pyplot as plt
+ 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
diff --git a/wetb/signal/tests/test_fit.py b/wetb/signal/tests/test_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..ecdfd2cc69796b268d7351a376552d2d401a23ef
--- /dev/null
+++ b/wetb/signal/tests/test_fit.py
@@ -0,0 +1,148 @@
+'''
+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
diff --git a/wetb/signal/tests/test_fouier_fit.py b/wetb/signal/tests/test_fouier_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..4992afafab0b40e0c86c016a0cc2938c88ddd792
--- /dev/null
+++ b/wetb/signal/tests/test_fouier_fit.py
@@ -0,0 +1,60 @@
+'''
+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()
+
diff --git a/wetb/signal_tools/tests/test_spectrum.py b/wetb/signal/tests/test_spectrum.py
similarity index 97%
rename from wetb/signal_tools/tests/test_spectrum.py
rename to wetb/signal/tests/test_spectrum.py
index fc20a3b8bb1177b401f601af12c8c5a39dc6f31c..8051c0407240956e18943ec5a4366247655fa827 100644
--- a/wetb/signal_tools/tests/test_spectrum.py
+++ b/wetb/signal/tests/test_spectrum.py
@@ -5,7 +5,7 @@ Created on 12/11/2015
'''
import unittest
import numpy as np
-from wetb.signal_tools.spectrum import psd
+from wetb.signal.spectrum import psd
class TestSpectrum(unittest.TestCase):
diff --git a/wetb/signal/tests/test_subset_mean.py b/wetb/signal/tests/test_subset_mean.py
new file mode 100644
index 0000000000000000000000000000000000000000..21f92445cba215d4e75357457e376fa85da1873a
--- /dev/null
+++ b/wetb/signal/tests/test_subset_mean.py
@@ -0,0 +1,86 @@
+'''
+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()
diff --git a/wetb/signal_tools/tests/test_first_order.py b/wetb/signal_tools/tests/test_first_order.py
deleted file mode 100644
index c49c464d36732c26bcf47a232ecb72ae75940b5c..0000000000000000000000000000000000000000
--- a/wetb/signal_tools/tests/test_first_order.py
+++ /dev/null
@@ -1,34 +0,0 @@
-'''
-Created on 13. jan. 2017
-
-@author: mmpe
-'''
-import unittest
-from wetb.signal_tools.filters import first_order
-import numpy as np
-
-class Test_first_order_filters(unittest.TestCase):
-
-
- def test_low_pass(self):
- a = np.random.randint(0,100,100).astype(np.float)
- b = first_order.low_pass(a, 1, 1)
- self.assertLess(b.std(), a.std())
- if 0:
- import matplotlib.pyplot as plt
- plt.plot(a)
- plt.plot(b)
- plt.show()
-
- def test_high_pass(self):
- a = np.random.randint(0,100,100).astype(np.float)
- b = first_order.high_pass(a, 1, 1)
- self.assertLess(b.mean(), a.mean())
- if 0:
- import matplotlib.pyplot as plt
- 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
diff --git a/wetb/utils/cluster_tools/cluster_resource.py b/wetb/utils/cluster_tools/cluster_resource.py
index 5c128cc59a8f84fafdd630cb3d31f0c1b17acfdc..9bec834ef61c5ef1097410cabc22b1ca81afeda6 100644
--- a/wetb/utils/cluster_tools/cluster_resource.py
+++ b/wetb/utils/cluster_tools/cluster_resource.py
@@ -10,9 +10,7 @@ import re
import threading
from wetb.utils.cluster_tools import pbswrap
-from wetb.utils.cluster_tools.ssh_client import SSHClient, SharedSSHClient
-from wetb.utils.timing import print_time
-import time
+
def unix_path(path, cwd=None, fail_on_missing=False):
@@ -119,9 +117,14 @@ class SSHPBSClusterResource(Resource):
@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
@@ -159,7 +162,11 @@ class SSHPBSClusterResource(Resource):
jobids = list(jobids)
self.ssh.execute("qdel %s" % (" ".join(jobids)))
-
+ 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""")
diff --git a/wetb/utils/cluster_tools/pbsjob.py b/wetb/utils/cluster_tools/pbsjob.py
index d1c1c99ec54246d467c42c94f17af921d6297dfb..0343b2de5a9a2d2bbf268f9bb86cf587f2f78ba2 100644
--- a/wetb/utils/cluster_tools/pbsjob.py
+++ b/wetb/utils/cluster_tools/pbsjob.py
@@ -4,7 +4,6 @@ Created on 04/12/2015
@author: mmpe
'''
import os
-from wetb.utils.cluster_tools.ssh_client import SSHClient
NOT_SUBMITTED = "Job not submitted"
PENDING = "Pending"
@@ -38,7 +37,6 @@ class SSHPBSJob(object):
@property
def status(self):
-
if self._status in [NOT_SUBMITTED, DONE]:
return self._status
with self.ssh:
diff --git a/wetb/utils/cluster_tools/ssh_client.py b/wetb/utils/cluster_tools/ssh_client.py
index 74ebd6b9f891371b3bde2fa2130b5b834415175b..66886fe362554a1f47e1b2d387821f64e06c95cc 100644
--- a/wetb/utils/cluster_tools/ssh_client.py
+++ b/wetb/utils/cluster_tools/ssh_client.py
@@ -5,9 +5,10 @@ Created on 27/11/2015
'''
from io import StringIO
+import sys
import paramiko
+
import os
-import sys
import threading
from _collections import deque
import time
diff --git a/wetb/wind/air_density.py b/wetb/wind/air_density.py
index 0af972420c99771bf88b5bf0d2f07b9aae196be7..13780781d4234b7047ce2da6ca5d4e1308f0f392 100644
--- a/wetb/wind/air_density.py
+++ b/wetb/wind/air_density.py
@@ -175,8 +175,8 @@ def R(rh=0, t=15, P=1014):
-------
Specific gas constant
"""
- assert np.all((900<P)&(P<1100)), "Pressure outside range 900 to 1100"
- assert np.all((-50<t)&(t<100)), "Temperature outside range -50 to 100"
+ #assert np.all((900<P)&(P<1100)), "Pressure outside range 900 to 1100"
+ #assert np.all((-50<t)&(t<100)), "Temperature outside range -50 to 100"
Tk = t + 273.15
return P * 100 / (air_density(P, t, rh) * Tk)