Refactoring for #47

CHANGELOG.md

@@ -1,8 +1,10 @@
 # Changelog
 
-## v1.0.5 (04.06.2021)
+## v1.1.0.dev0
 
-* Introduced optional 'normalizedsmooth' argument to reduce dependence on xdata and weights
+_Unreleased_
+
+* Introduced optional 'normalizedsmooth' argument to reduce dependence on xdata and weights [#47](https://github.com/espdev/csaps/pull/47)
 
 ## v1.0.4 (04.05.2021)
 

csaps/_sspumv.py

@@ -240,6 +240,23 @@ def trace(m: sp.dia_matrix):
 
         return 1. / (1. + trace(a) / (6. * trace(b)))
 
+    @staticmethod
+    def _normalize_smooth(x: np.ndarray, w: np.ndarray, smooth: Optional[float]):
+        """
+        See the explanation here: https://github.com/espdev/csaps/pull/47
+        """
+
+        span = np.ptp(x)
+
+        eff_x = 1 + (span ** 2) / np.sum(np.diff(x) ** 2)
+        eff_w = np.sum(w) ** 2 / np.sum(w ** 2)
+        k = 80 * (span ** 3) * (x.size ** -2) * (eff_x ** -0.5) * (eff_w ** -0.5)
+
+        s = 0.5 if smooth is None else smooth
+        p = s / (s + (1 - s) * k)
+
+        return p
+
     @staticmethod
     def _make_spline(x, y, w, smooth, shape, normalizedsmooth):
         pcount = x.size
@@ -276,14 +293,9 @@ def _make_spline(x, y, w, smooth, shape, normalizedsmooth):
         qtw = qtw @ qtw.T
 
         p = smooth
+
         if normalizedsmooth:
-            span = np.ptp(x)
-            count = x.size
-            eff_x = 1+(np.ptp(x)**2)/np.sum(np.diff(x)**2)
-            eff_w = np.sum(w)**2 / np.sum(w**2)
-            K = (80)*(span**3)*(count**-2)*(eff_x**-0.5)*(eff_w**-0.5)
-            s = 0.5 if smooth is None else smooth
-            p = s/(s+(1-s)*K)
+            p = CubicSmoothingSpline._normalize_smooth(x, w, smooth)
         elif smooth is None:
             p = CubicSmoothingSpline._compute_smooth(r, qtw)
 

csaps/_version.py

@@ -1,3 +1,3 @@
 # -*- coding: utf-8 -*-
 
-__version__ = '1.0.4'
+__version__ = '1.1.0.dev0'

tests/test_csaps.py

@@ -73,20 +73,6 @@ def test_shortcut_output(data, tolist):
     sp = csaps(x, y, normalizedsmooth=True)
     assert isinstance(sp, sp_cls)
 
-    if tolist:
-        x2 = (2*np.array(x)).tolist()
-        xi2 = (2*np.array(xi)).tolist()
-    else:
-        x2 = ([2*np.array(xx, dtype=np.float64) for xx in x]
-              if isinstance(x, list)
-              else 2*np.array(x, dtype=np.float64))
-        xi2 = ([2*np.array(xx, dtype=np.float64) for xx in xi]
-               if isinstance(x, list)
-               else 2*np.array(xi, dtype=np.float64))
-    smoothed_dataA = csaps(x, y, xi, normalizedsmooth=True, smooth=0.25)
-    smoothed_dataB = csaps(x2, y, xi2, normalizedsmooth=True, smooth=0.25)
-    assert np.all(np.isclose(smoothed_dataA, smoothed_dataB))
-
 
 @pytest.mark.parametrize('smooth, cls', [
     (0.85, np.ndarray),
@@ -101,3 +87,28 @@ def test_shortcut_ndgrid_smooth_output(surface, smooth, cls):
 
     output = csaps(x, y, x, smooth=smooth)
     assert isinstance(output, cls)
+
+
+@pytest.mark.parametrize('data', ['univariate', 'ndgrid'], indirect=True)
+@pytest.mark.parametrize('smooth', [0.0, 0.1, 0.25, 0.5, 0.85, 1.0, None])
+@pytest.mark.parametrize('scale', [0.001, 0.005, 0.25, 0.5, 2.0, 5.0, 100.0, 1000.0])
+def test_normalized_smooth(data, smooth, scale):
+    x, y, xi, *_ = data
+
+    x2 = (
+        [scale * np.array(xx, dtype=np.float64) for xx in x]
+        if isinstance(x, list) else scale * np.array(x, dtype=np.float64)
+    )
+    xi2 = (
+        [scale * np.array(xx, dtype=np.float64) for xx in xi]
+        if isinstance(x, list) else scale * np.array(xi, dtype=np.float64)
+    )
+
+    smoothed_data_a = csaps(x, y, xi, smooth=smooth, normalizedsmooth=True)
+    smoothed_data_b = csaps(x2, y, xi2, smooth=smooth, normalizedsmooth=True)
+
+    if isinstance(smoothed_data_a, AutoSmoothingResult):
+        smoothed_data_a = smoothed_data_a.values
+        smoothed_data_b = smoothed_data_b.values
+
+    assert smoothed_data_a == pytest.approx(smoothed_data_b, rel=1e-2)