From 6257bc6a9636826309451c995ebe342c81ffea12 Mon Sep 17 00:00:00 2001
From: Per A Brodtkorb <per.andreas.brodtkorb@gmail.com>
Date: Fri, 26 Feb 2016 12:43:31 +0100
Subject: [PATCH] Simplified HampelFilter

---
 wafo/sg_filter.py | 117 +++++++++++++++++++++++++---------------------
 1 file changed, 64 insertions(+), 53 deletions(-)

diff --git a/wafo/sg_filter.py b/wafo/sg_filter.py
index 28b9f60..eca2050 100644
--- a/wafo/sg_filter.py
+++ b/wafo/sg_filter.py
@@ -438,8 +438,8 @@ def smoothn(data, s=None, weight=None, robust=False, z0=None, tolz=1e-3,
     W = W * IsFinite
     if (W < 0).any():
         raise ValueError('Weights must all be >=0')
-    else:
-        W = W / W.max()
+
+    W = W / W.max()
 
     isweighted = (W < 1).any()  # Weighted or missing data?
     isauto = s is None  # Automatic smoothing?
@@ -1217,70 +1217,81 @@ class HampelFilter(object):
         self.adaptive = adaptive
         self.fulloutput = fulloutput
 
-    def __call__(self, y, x=None):
-        Y = np.atleast_1d(y).ravel()
-        if x is None:
-            x = range(len(Y))
-        X = np.atleast_1d(x).ravel()
-
-        dx = self.dx
-        if dx is None:
-            dx = 3 * np.median(np.diff(X))
+    def _check(self, dx):
         if not np.isscalar(dx):
             raise ValueError('DX must be a scalar.')
-        elif dx < 0:
+        if dx < 0:
             raise ValueError('DX must be larger than zero.')
 
-        YY = Y
-        S0 = np.nan * np.zeros(YY.shape)
-        Y0 = np.nan * np.zeros(YY.shape)
+    @staticmethod
+    def localwindow(X, Y, DX, i):
+        mask = (X[i] - DX <= X) & (X <= X[i] + DX)
+        Y0 = np.median(Y[mask])
+        # Calculate Local Scale of Natural Variation
+        S0 = 1.4826 * np.median(np.abs(Y[mask] - Y0))
+        return Y0, S0
+
+    @staticmethod
+    def smgauss(X, V, DX):
+        Xj = X
+        Xk = np.atleast_2d(X).T
+        Wjk = np.exp(-((Xj - Xk) / (2 * DX)) ** 2)
+        G = np.dot(Wjk, V) / np.sum(Wjk, axis=0)
+        return G
+
+    def _adaptive(self, Y, X, dx):
+        localwindow = self.localwindow
+        Y0, S0, ADX = self._init(Y, dx)
+        Y0Tmp = np.nan * np.zeros(Y.shape)
+        S0Tmp = np.nan * np.zeros(Y.shape)
+        DXTmp = np.arange(1, len(S0) + 1) * dx
+        # Integer variation of Window Half Size
+        # Calculate Initial Guess of Optimal Parameters Y0, S0, ADX
+        for i in range(len(Y)):
+            j = 0
+            S0Rel = np.inf
+            while S0Rel > self.adaptive:
+                Y0Tmp[j], S0Tmp[j] = localwindow(X, Y, DXTmp[j], i)
+                if j > 0:
+                    S0Rel = abs((S0Tmp[j - 1] - S0Tmp[j]) /
+                                (S0Tmp[j - 1] + S0Tmp[j]) / 2)
+                j += 1
+
+            Y0[i] = Y0Tmp[j - 2]
+            S0[i] = S0Tmp[j - 2]
+            ADX[i] = DXTmp[j - 2] / dx
+
+    # Gaussian smoothing of relevant parameters
+        DX = 2 * np.median(np.diff(X))
+        ADX = self.smgauss(X, ADX, DX)
+        S0 = self.smgauss(X, S0, DX)
+        Y0 = self.smgauss(X, Y0, DX)
+        return Y0, S0, ADX
+
+    def _init(self, Y, dx):
+        S0 = np.nan * np.zeros(Y.shape)
+        Y0 = np.nan * np.zeros(Y.shape)
         ADX = dx * np.ones(Y.shape)
+        return Y0, S0, ADX
 
-        def localwindow(X, Y, DX, i):
-            mask = (X[i] - DX <= X) & (X <= X[i] + DX)
-            Y0 = np.median(Y[mask])
-            # Calculate Local Scale of Natural Variation
-            S0 = 1.4826 * np.median(np.abs(Y[mask] - Y0))
-            return Y0, S0
+    def __call__(self, y, x=None):
+        Y = np.atleast_1d(y).ravel()
+        if x is None:
+            x = range(len(Y))
+        X = np.atleast_1d(x).ravel()
 
-        def smgauss(X, V, DX):
-            Xj = X
-            Xk = np.atleast_2d(X).T
-            Wjk = np.exp(-((Xj - Xk) / (2 * DX)) ** 2)
-            G = np.dot(Wjk, V) / np.sum(Wjk, axis=0)
-            return G
+        dx = 3 * np.median(np.diff(X)) if self.dx is None else self.dx
+        self._check(dx)
 
         if len(X) > 1:
             if self.adaptive is None:
+                localwindow = self.localwindow
+                Y0, S0, ADX = self._init(Y, dx)
                 for i in range(len(Y)):
                     Y0[i], S0[i] = localwindow(X, Y, dx, i)
             else:  # 'adaptive'
-
-                Y0Tmp = np.nan * np.zeros(YY.shape)
-                S0Tmp = np.nan * np.zeros(YY.shape)
-                DXTmp = np.arange(1, len(S0) + 1) * dx
-                # Integer variation of Window Half Size
-
-                # Calculate Initial Guess of Optimal Parameters Y0, S0, ADX
-                for i in range(len(Y)):
-                    j = 0
-                    S0Rel = np.inf
-                    while S0Rel > self.adaptive:
-                        Y0Tmp[j], S0Tmp[j] = localwindow(X, Y, DXTmp[j], i)
-                        if j > 0:
-                            S0Rel = abs((S0Tmp[j - 1] - S0Tmp[j]) /
-                                        (S0Tmp[j - 1] + S0Tmp[j]) / 2)
-                        j += 1
-                    Y0[i] = Y0Tmp[j - 2]
-                    S0[i] = S0Tmp[j - 2]
-                    ADX[i] = DXTmp[j - 2] / dx
-
-                # Gaussian smoothing of relevant parameters
-                DX = 2 * np.median(np.diff(X))
-                ADX = smgauss(X, ADX, DX)
-                S0 = smgauss(X, S0, DX)
-                Y0 = smgauss(X, Y0, DX)
-
+                Y0, S0, ADX = self._adaptive(Y, X, dx)
+        YY = Y.copy()
         T = self.t
         # Prepare Output
         self.UB = Y0 + T * S0