diff --git a/wafo/sg_filter/_core.py b/wafo/sg_filter/_core.py index f23dea4..7e21ca1 100644 --- a/wafo/sg_filter/_core.py +++ b/wafo/sg_filter/_core.py @@ -137,7 +137,7 @@ class SavitzkyGolay(object): Cambridge University Press ISBN-13: 9780521880688 """ - def __init__(self, n, degree=1, diff_order=0, delta=1.0, axis=-1, + def __init__(self, n, degree=1, diff_order=0, delta=1.0, axis=-1, mode='interp', cval=0.0): self.n = n self.degree = degree @@ -250,7 +250,6 @@ def evar(y): S = np.zeros(sh0) sh1 = np.ones((d,)) cos = np.cos - pi = np.pi for i in range(d): ni = sh0[i] sh1[i] = ni @@ -314,7 +313,10 @@ class _Filter(object): return (np.array(y.shape) != 1).sum() @staticmethod - def _smoothness_limits(n): + def _smoothness_par(h): + return (((1 + sqrt(1 + 8 * h)) / 4. / h) ** 2 - 1) / 16 + + def _smoothness_limits(self, n): """ Return upper and lower bound for the smoothness parameter @@ -327,8 +329,8 @@ class _Filter(object): h_min = 1e-6 ** (2. / n) h_max = 0.99 ** (2. / n) - s_min = (((1 + sqrt(1 + 8 * h_max)) / 4. / h_max) ** 2 - 1) / 16 - s_max = (((1 + sqrt(1 + 8 * h_min)) / 4. / h_min) ** 2 - 1) / 16 + s_min = self._smoothness_par(h_max) + s_max = self._smoothness_par(h_min) return s_min, s_max @staticmethod @@ -354,6 +356,8 @@ class _Filter(object): Lambda = self._lambda_tensor(y) def gamma(s): + if s is None: + return 1.0 return 1. / (1 + s * Lambda ** 2) return gamma @@ -450,30 +454,28 @@ class _Filter(object): def gcv(self, p, aow, DCTy, y, Wtot): # Search the smoothing parameter s that minimizes the GCV score s = 10.0 ** p - Gamma = self.gamma(s) + gamma_s = self.gamma(s) if aow > 0.9: # aow = 1 means that all of the data are equally weighted # very much faster: does not require any inverse DCT - residual = DCTy.ravel() * (Gamma.ravel() - 1) + residual = DCTy.ravel() * (gamma_s.ravel() - 1) else: # take account of the weights to calculate RSS: is_finite = self.is_finite - yhat = idctn(Gamma * DCTy) + yhat = idctn(gamma_s * DCTy) residual = sqrt(Wtot[is_finite]) * (y[is_finite] - yhat[is_finite]) - TrH = Gamma.sum() + TrH = gamma_s.sum() RSS = linalg.norm(residual)**2 # Residual sum-of-squares GCVscore = RSS / self.nof / (1.0 - TrH / y.size) ** 2 return GCVscore - def __call__(self, z, s): + def _smooth(self, z, s): auto_smooth = self.auto_smooth norm = linalg.norm y = self.y Wtot = self.Wtot - Gamma = 1 - if s is not None: - Gamma = self.gamma(s) + gamma_s = self.gamma(s) # "amount" of weights (see the function GCVscore) aow = Wtot.sum() / y.size # 0 < aow <= 1 for nit in range(self.maxiter): @@ -489,18 +491,23 @@ class _Filter(object): args=(aow, DCTy, y, Wtot), xtol=self.errp, full_output=False, disp=False) s = 10 ** log10s - Gamma = self.gamma(s) + gamma_s = self.gamma(s) z0 = z - z = self.RF * idctn(Gamma * DCTy) + (1 - self.RF) * z + z = self.RF * idctn(gamma_s * DCTy) + (1 - self.RF) * z # if no weighted/missing data => tol=0 (no iteration) tol = norm(z0.ravel() - z.ravel()) / norm(z.ravel()) converged = tol <= self.tolz or not self.is_weighted if converged: break + return z, s, converged + + def __call__(self, z, s): + z, s, converged = self._smooth(z, s) if self.robust: # -- Robust Smoothing: iteratively re-weighted process h = self._average_leverage(s, self.N) - self.Wtot = self.W * self.robust_weights(y - z, self.is_finite, h) + self.Wtot = self.W * self.robust_weights(self.y - z, + self.is_finite, h) # re-initialize for another iterative weighted process self.is_weighted = True return z, s, converged @@ -919,25 +926,58 @@ class Kalman(object): def reset(self): self._filter = self._filter_first + def _none_or_atleast_2d(self, a): + if a is not None: + return np.atleast_2d(a) + return a + + @property + def A(self): + return self._a + + @A.setter + def A(self, a): + self._a = self._none_or_atleast_2d(a) + def _set_A(self, n): if self.A is None: self.A = np.eye(n) - self.A = np.atleast_2d(self.A) + + @property + def Q(self): + return self._q + + @Q.setter + def Q(self, q): + self._q = self._none_or_atleast_2d(q) def _set_Q(self, n): if self.Q is None: self.Q = np.zeros((n, n)) - self.Q = np.atleast_2d(self.Q) + + @property + def H(self): + return self._h + + @H.setter + def H(self, h): + self._h = self._none_or_atleast_2d(h) def _set_H(self, n): if self.H is None: self.H = np.eye(n) - self.H = np.atleast_2d(self.H) + + @property + def P(self): + return self._p + + @P.setter + def P(self, p): + self._p = self._none_or_atleast_2d(p) def _set_P(self, HI): if self.P is None: self.P = np.dot(np.dot(HI, self.R), HI.T) - self.P = np.atleast_2d(self.P) def _init_first(self, n): self._set_A(n) @@ -1165,7 +1205,7 @@ class HampelFilter(object): self.UB = Y0 + T * S0 self.LB = Y0 - T * S0 outliers = np.abs(Y - Y0) > T * S0 # possible outliers - np.putmask(YY, outliers, Y0) # YY[outliers] = Y0[outliers] + np.putmask(YY, outliers, Y0) # YY[outliers] = Y0[outliers] self.outliers = outliers self.num_outliers = outliers.sum() self.ADX = ADX diff --git a/wafo/sg_filter/tests/test_sg_filter.py b/wafo/sg_filter/tests/test_sg_filter.py index 51f6136..aa53c26 100644 --- a/wafo/sg_filter/tests/test_sg_filter.py +++ b/wafo/sg_filter/tests/test_sg_filter.py @@ -1,9 +1,3 @@ -''' -Created on 26. feb. 2016 - -@author: pab -''' - import unittest import numpy as np from wafo.sg_filter._core import HampelFilter, SavitzkyGolay, Kalman @@ -42,19 +36,12 @@ class Test(unittest.TestCase): def test_savitzky_golay(self): t = np.linspace(-4, 4, 100) # noise = np.random.normal(0, 0.05, t.shape) - noise = np.array( - [-0.0100738410948, 0.0434144660993, -0.00378151879776, - 0.0113356386249, 0.00435077505515, -0.0477602072563, - 0.0373494077292, -0.0133865881648, -0.058259106898, - 0.0661995082581, -0.0302004598725, 0.0251073570423, - 0.0504480752083, 0.00591197638959, -0.0936100824835, - -0.0475801028574, 0.164260739963, 0.000309393717298, - 0.0145265378376, 0.0451974209359, -0.0296464887607, 0.00113264461646, 0.0203243951278, 0.0250086489045, 0.0286663554031, 0.0162585459204, -0.0285237485127, 0.0841528003829, -0.0449728982279, 0.0549796633406, -0.0968901623365, 0.0306727572312, -0.06950957371, 0.00598396964893, 0.0583904237732, 0.0271856795988, 0.0503237614252, 0.179978222833, 0.0909965842309, -0.037723727536, 0.0766788834896, -0.135152061379, -0.0275445486266, 0.00437867535899, -0.0272885999582, -0.0167170993005, 0.0522077667657, -0.044618733434, -0.0323362657615, 0.0122726298648, 0.0311304227245, -0.106361849439, 0.0526578516171, 0.0600674983294, 0.0745276212563, -0.0488453979797, -0.000833273665862, 0.0232768529385, 0.0138869047145, -0.0581069396022, 0.0321001340163, -0.023975724861, -0.0281919023447, 0.0341550408589, -0.0128104984676, -0.0214928514613, -0.0930630211028, -0.086094936458, 0.0169055851217, -0.0106085935018, 0.0114499301052, 0.120823119533, 0.0401375561218, 0.0979979215035, -0.015930504967, 0.0730359376807, -0.0369380623571, 0.0628171477772, -0.0271592027456, -0.0805267873905, -0.047314352792, -0.00693573205681, -0.0276294819462, -0.00896110691771, 0.0248482543453, 0.0877066047642, -0.0325533337143, -0.0638923240986, 0.0512523528302, 0.0720952473332, 0.00555629231342, 0.0562912258368, 0.0518574742355, -0.0599288505845, 0.0129657160415, -0.0025435076237, -0.0136966498082, 0.0260013064028, 0.0315561663839, -0.0293194006299]) # @IgnorePep8 + noise = np.array([-0.0100738410948, 0.0434144660993, -0.00378151879776, 0.0113356386249, 0.00435077505515, -0.0477602072563, 0.0373494077292, -0.0133865881648, -0.058259106898, 0.0661995082581, -0.0302004598725, 0.0251073570423, 0.0504480752083, 0.00591197638959, -0.0936100824835, -0.0475801028574, 0.164260739963, 0.000309393717298, 0.0145265378376, 0.0451974209359, -0.0296464887607, 0.00113264461646, 0.0203243951278, 0.0250086489045, 0.0286663554031, 0.0162585459204, -0.0285237485127, 0.0841528003829, -0.0449728982279, 0.0549796633406, -0.0968901623365, 0.0306727572312, -0.06950957371, 0.00598396964893, 0.0583904237732, 0.0271856795988, 0.0503237614252, 0.179978222833, 0.0909965842309, -0.037723727536, 0.0766788834896, -0.135152061379, -0.0275445486266, 0.00437867535899, -0.0272885999582, -0.0167170993005, 0.0522077667657, -0.044618733434, -0.0323362657615, 0.0122726298648, 0.0311304227245, -0.106361849439, 0.0526578516171, 0.0600674983294, 0.0745276212563, -0.0488453979797, -0.000833273665862, 0.0232768529385, 0.0138869047145, -0.0581069396022, 0.0321001340163, -0.023975724861, -0.0281919023447, 0.0341550408589, -0.0128104984676, -0.0214928514613, -0.0930630211028, -0.086094936458, 0.0169055851217, -0.0106085935018, 0.0114499301052, 0.120823119533, 0.0401375561218, 0.0979979215035, -0.015930504967, 0.0730359376807, -0.0369380623571, 0.0628171477772, -0.0271592027456, -0.0805267873905, -0.047314352792, -0.00693573205681, -0.0276294819462, -0.00896110691771, 0.0248482543453, 0.0877066047642, -0.0325533337143, -0.0638923240986, 0.0512523528302, 0.0720952473332, 0.00555629231342, 0.0562912258368, 0.0518574742355, -0.0599288505845, 0.0129657160415, -0.0025435076237, -0.0136966498082, 0.0260013064028, 0.0315561663839, -0.0293194006299]) # nopep8 # print(', '.join(['{0}'.format(y) for y in noise])) y = np.exp(-t**2) + noise ysg = SavitzkyGolay(n=20, degree=2).smooth(y) # print(', '.join(['{0}'.format(y) for y in ysg])) - assert_array_almost_equal(ysg,[0.08136940307, 0.0628995512206, 0.0459066440291, 0.0303906814954, 0.0163516636196, 0.00378959040169, -0.00729553815838, -0.0169037220606, -0.0250349613049, -0.0316892558914, -0.03686660582, -0.0405670110908, -0.0427904717036, -0.0435369876586, -0.0428065589558, -0.0405991855951, -0.0369148675765, -0.0317536049001, -0.0251153975657, -0.0170002455736, -0.00740814892354, 0.00138883298422, 0.00903982318392, 0.0156981288342, 0.0262571138525, 0.038712235996, 0.0497401536595, 0.0708145615737, 0.0930139468092, 0.115399487253, 0.145859174356, 0.181612636791, 0.217133614374, 0.257717904903, 0.300137304571, 0.345750694351, 0.392370647909, 0.447698806473, 0.498669071812, 0.552242388661, 0.602481962016, 0.649961345551, 0.695325299394, 0.735374047294, 0.773010811455, 0.805919165272, 0.834865852047, 0.854796400219, 0.86842923705, 0.871717227059, 0.870973756546, 0.854091989437, 0.838398655214, 0.810673923439, 0.783634317013, 0.749908334726, 0.714315301461, 0.67175727743, 0.634614314389, 0.593943523363, 0.54497174327, 0.497539947485, 0.442464744684, 0.39039262486, 0.339641566641, 0.287962531039, 0.243889096129, 0.20650960823, 0.165038163827, 0.127423305227, 0.098217909507, 0.0723906956301, 0.0457840790019, 0.0325691800506, 0.0207065223098, 0.0137400799931, 0.00547944733578, -0.000268557810876, -0.00301934919753, -0.00213643706217, -0.0101846106787, -0.0170002602175, -0.0225833856784, -0.0269339870616, -0.030052064367, -0.0319376175946, -0.0325906467444, -0.0320111518164, -0.0301991328106, -0.0271545897271, -0.0228775225657, -0.0173679313266, -0.0106258160097, -0.00265117661497, 0.00655598685753, 0.0169956744078, 0.0286678860359, 0.0415726217417, 0.0557098815254, 0.0710796653868]) # @IgnorePep8 + assert_array_almost_equal(ysg, [0.08136940307, 0.0628995512206, 0.0459066440291, 0.0303906814954, 0.0163516636196, 0.00378959040169, -0.00729553815838, -0.0169037220606, -0.0250349613049, -0.0316892558914, -0.03686660582, -0.0405670110908, -0.0427904717036, -0.0435369876586, -0.0428065589558, -0.0405991855951, -0.0369148675765, -0.0317536049001, -0.0251153975657, -0.0170002455736, -0.00740814892354, 0.00138883298422, 0.00903982318392, 0.0156981288342, 0.0262571138525, 0.038712235996, 0.0497401536595, 0.0708145615737, 0.0930139468092, 0.115399487253, 0.145859174356, 0.181612636791, 0.217133614374, 0.257717904903, 0.300137304571, 0.345750694351, 0.392370647909, 0.447698806473, 0.498669071812, 0.552242388661, 0.602481962016, 0.649961345551, 0.695325299394, 0.735374047294, 0.773010811455, 0.805919165272, 0.834865852047, 0.854796400219, 0.86842923705, 0.871717227059, 0.870973756546, 0.854091989437, 0.838398655214, 0.810673923439, 0.783634317013, 0.749908334726, 0.714315301461, 0.67175727743, 0.634614314389, 0.593943523363, 0.54497174327, 0.497539947485, 0.442464744684, 0.39039262486, 0.339641566641, 0.287962531039, 0.243889096129, 0.20650960823, 0.165038163827, 0.127423305227, 0.098217909507, 0.0723906956301, 0.0457840790019, 0.0325691800506, 0.0207065223098, 0.0137400799931, 0.00547944733578, -0.000268557810876, -0.00301934919753, -0.00213643706217, -0.0101846106787, -0.0170002602175, -0.0225833856784, -0.0269339870616, -0.030052064367, -0.0319376175946, -0.0325906467444, -0.0320111518164, -0.0301991328106, -0.0271545897271, -0.0228775225657, -0.0173679313266, -0.0106258160097, -0.00265117661497, 0.00655598685753, 0.0169956744078, 0.0286678860359, 0.0415726217417, 0.0557098815254, 0.0710796653868]) # nopep8 def test_hampelfilter(self):