danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Reduced complexity + pep8

Browse Source
master
Per A Brodtkorb 8 years ago
parent a54aa6a323
commit 65bfb085d3
2 changed files with 63 additions and 36 deletions
Show Stats Download Patch File Download Diff File

78
wafo/sg_filter/_core.py
Unescape Escape View File

@ -250,7 +250,6 @@ def evar(y):
S = np.zeros(sh0) S = np.zeros(sh0)
sh1 = np.ones((d,)) sh1 = np.ones((d,))
cos = np.cos cos = np.cos
pi = np.pi
for i in range(d): for i in range(d):
ni = sh0[i] ni = sh0[i]
sh1[i] = ni sh1[i] = ni
@ -314,7 +313,10 @@ class _Filter(object):
return (np.array(y.shape) != 1).sum() return (np.array(y.shape) != 1).sum()
@staticmethod @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 Return upper and lower bound for the smoothness parameter
@ -327,8 +329,8 @@ class _Filter(object):
h_min = 1e-6 ** (2. / n) h_min = 1e-6 ** (2. / n)
h_max = 0.99 ** (2. / n) h_max = 0.99 ** (2. / n)
s_min = (((1 + sqrt(1 + 8 * h_max)) / 4. / h_max) ** 2 - 1) / 16 s_min = self._smoothness_par(h_max)
s_max = (((1 + sqrt(1 + 8 * h_min)) / 4. / h_min) ** 2 - 1) / 16 s_max = self._smoothness_par(h_min)
return s_min, s_max return s_min, s_max
@staticmethod @staticmethod
@ -354,6 +356,8 @@ class _Filter(object):
Lambda = self._lambda_tensor(y) Lambda = self._lambda_tensor(y)
def gamma(s): def gamma(s):
if s is None:
return 1.0
return 1. / (1 + s * Lambda ** 2) return 1. / (1 + s * Lambda ** 2)
return gamma return gamma
@ -450,30 +454,28 @@ class _Filter(object):
def gcv(self, p, aow, DCTy, y, Wtot): def gcv(self, p, aow, DCTy, y, Wtot):
# Search the smoothing parameter s that minimizes the GCV score # Search the smoothing parameter s that minimizes the GCV score
s = 10.0 ** p s = 10.0 ** p
Gamma = self.gamma(s) gamma_s = self.gamma(s)
if aow > 0.9: if aow > 0.9:
# aow = 1 means that all of the data are equally weighted # aow = 1 means that all of the data are equally weighted
# very much faster: does not require any inverse DCT # very much faster: does not require any inverse DCT
residual = DCTy.ravel() * (Gamma.ravel() - 1) residual = DCTy.ravel() * (gamma_s.ravel() - 1)
else: else:
# take account of the weights to calculate RSS: # take account of the weights to calculate RSS:
is_finite = self.is_finite 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]) 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 RSS = linalg.norm(residual)**2 # Residual sum-of-squares
GCVscore = RSS / self.nof / (1.0 - TrH / y.size) ** 2 GCVscore = RSS / self.nof / (1.0 - TrH / y.size) ** 2
return GCVscore return GCVscore
def __call__(self, z, s): def _smooth(self, z, s):
auto_smooth = self.auto_smooth auto_smooth = self.auto_smooth
norm = linalg.norm norm = linalg.norm
y = self.y y = self.y
Wtot = self.Wtot Wtot = self.Wtot
Gamma = 1 gamma_s = self.gamma(s)
if s is not None:
Gamma = self.gamma(s)
# "amount" of weights (see the function GCVscore) # "amount" of weights (see the function GCVscore)
aow = Wtot.sum() / y.size # 0 < aow <= 1 aow = Wtot.sum() / y.size # 0 < aow <= 1
for nit in range(self.maxiter): for nit in range(self.maxiter):
@ -489,18 +491,23 @@ class _Filter(object):
args=(aow, DCTy, y, Wtot), args=(aow, DCTy, y, Wtot),
xtol=self.errp, full_output=False, disp=False) xtol=self.errp, full_output=False, disp=False)
s = 10 ** log10s s = 10 ** log10s
Gamma = self.gamma(s) gamma_s = self.gamma(s)
z0 = z 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) # if no weighted/missing data => tol=0 (no iteration)
tol = norm(z0.ravel() - z.ravel()) / norm(z.ravel()) tol = norm(z0.ravel() - z.ravel()) / norm(z.ravel())
converged = tol <= self.tolz or not self.is_weighted converged = tol <= self.tolz or not self.is_weighted
if converged: if converged:
break break
return z, s, converged
def __call__(self, z, s):
z, s, converged = self._smooth(z, s)
if self.robust: if self.robust:
# -- Robust Smoothing: iteratively re-weighted process # -- Robust Smoothing: iteratively re-weighted process
h = self._average_leverage(s, self.N) 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 # re-initialize for another iterative weighted process
self.is_weighted = True self.is_weighted = True
return z, s, converged return z, s, converged
@ -919,25 +926,58 @@ class Kalman(object):
def reset(self): def reset(self):
self._filter = self._filter_first 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): def _set_A(self, n):
if self.A is None: if self.A is None:
self.A = np.eye(n) 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): def _set_Q(self, n):
if self.Q is None: if self.Q is None:
self.Q = np.zeros((n, n)) 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): def _set_H(self, n):
if self.H is None: if self.H is None:
self.H = np.eye(n) 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): def _set_P(self, HI):
if self.P is None: if self.P is None:
self.P = np.dot(np.dot(HI, self.R), HI.T) self.P = np.dot(np.dot(HI, self.R), HI.T)
self.P = np.atleast_2d(self.P)
def _init_first(self, n): def _init_first(self, n):
self._set_A(n) self._set_A(n)

17
wafo/sg_filter/tests/test_sg_filter.py
Unescape Escape View File

@ -1,9 +1,3 @@
'''
Created on 26. feb. 2016
@author: pab
'''
import unittest import unittest
import numpy as np import numpy as np
from wafo.sg_filter._core import HampelFilter, SavitzkyGolay, Kalman from wafo.sg_filter._core import HampelFilter, SavitzkyGolay, Kalman
@ -42,19 +36,12 @@ class Test(unittest.TestCase):
def test_savitzky_golay(self): def test_savitzky_golay(self):
t = np.linspace(-4, 4, 100) t = np.linspace(-4, 4, 100)
# noise = np.random.normal(0, 0.05, t.shape) # noise = np.random.normal(0, 0.05, t.shape)
noise = np.array( 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
[-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
# print(', '.join(['{0}'.format(y) for y in noise])) # print(', '.join(['{0}'.format(y) for y in noise]))
y = np.exp(-t**2) + noise y = np.exp(-t**2) + noise
ysg = SavitzkyGolay(n=20, degree=2).smooth(y) ysg = SavitzkyGolay(n=20, degree=2).smooth(y)
# print(', '.join(['{0}'.format(y) for y in ysg])) # 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): def test_hampelfilter(self):

Write Preview
Loading…
Cancel
Save