Making CubicSmoothingSpline smoothing parameter independent of xdata scale

CHANGELOG.md

@@ -1,5 +1,9 @@
 # Changelog
 
+## v1.0.5 (04.06.2021)
+
+* Introduced optional 'normalizedsmooth' argument to reduce dependence on xdata and weights
+
 ## v1.0.4 (04.05.2021)
 
 * Bump numpy dependency version

csaps/_shortcut.py

@@ -33,7 +33,8 @@ def csaps(xdata: UnivariateDataType,
           *,
           weights: Optional[UnivariateDataType] = None,
           smooth: Optional[float] = None,
-          axis: Optional[int] = None) -> ISmoothingSpline: ...
+          axis: Optional[int] = None,
+          normalizedsmooth: bool = False) -> ISmoothingSpline: ...
 
 
 @overload
@@ -42,7 +43,8 @@ def csaps(xdata: UnivariateDataType,
           xidata: UnivariateDataType,
           *,
           weights: Optional[UnivariateDataType] = None,
-          axis: Optional[int] = None) -> AutoSmoothingResult: ...
+          axis: Optional[int] = None,
+          normalizedsmooth: bool = False) -> AutoSmoothingResult: ...
 
 
 @overload
@@ -52,7 +54,8 @@ def csaps(xdata: UnivariateDataType,
           *,
           smooth: float,
           weights: Optional[UnivariateDataType] = None,
-          axis: Optional[int] = None) -> MultivariateDataType: ...
+          axis: Optional[int] = None,
+          normalizedsmooth: bool = False) -> MultivariateDataType: ...
 
 
 @overload
@@ -61,7 +64,8 @@ def csaps(xdata: NdGridDataType,
           *,
           weights: Optional[NdGridDataType] = None,
           smooth: Optional[Sequence[float]] = None,
-          axis: Optional[int] = None) -> ISmoothingSpline: ...
+          axis: Optional[int] = None,
+          normalizedsmooth: bool = False) -> ISmoothingSpline: ...
 
 
 @overload
@@ -70,7 +74,8 @@ def csaps(xdata: NdGridDataType,
           xidata: NdGridDataType,
           *,
           weights: Optional[NdGridDataType] = None,
-          axis: Optional[int] = None) -> AutoSmoothingResult: ...
+          axis: Optional[int] = None,
+          normalizedsmooth: bool = False) -> AutoSmoothingResult: ...
 
 
 @overload
@@ -80,7 +85,8 @@ def csaps(xdata: NdGridDataType,
           *,
           smooth: Sequence[float],
           weights: Optional[NdGridDataType] = None,
-          axis: Optional[int] = None) -> MultivariateDataType: ...
+          axis: Optional[int] = None,
+          normalizedsmooth: bool = False) -> MultivariateDataType: ...
 #
 # csaps signatures
 # **************************************
@@ -92,7 +98,8 @@ def csaps(xdata: Union[UnivariateDataType, NdGridDataType],
           *,
           weights: Optional[Union[UnivariateDataType, NdGridDataType]] = None,
           smooth: Optional[Union[float, Sequence[float]]] = None,
-          axis: Optional[int] = None) -> Union[MultivariateDataType, ISmoothingSpline, AutoSmoothingResult]:
+          axis: Optional[int] = None,
+          normalizedsmooth: bool = False) -> Union[MultivariateDataType, ISmoothingSpline, AutoSmoothingResult]:
     """Smooths the univariate/multivariate/gridded data or computes the corresponding splines
 
     This function might be used as the main API for smoothing any data.
@@ -139,6 +146,10 @@ def csaps(xdata: Union[UnivariateDataType, NdGridDataType],
         .. note::
             Currently, `axis` will be ignored for nd-gridded ``ydata`` case.
 
+    normalizedsmooth : [*Optional*] bool
+        If True, the smooth parameter is normalized such that results are invariant to xdata range
+        and less sensitive to nonuniformity of weights and xdata clumping
+
     Returns
     -------
 
@@ -200,9 +211,10 @@ def csaps(xdata: Union[UnivariateDataType, NdGridDataType],
 
     if umv:
         axis = -1 if axis is None else axis
-        sp = CubicSmoothingSpline(xdata, ydata, weights=weights, smooth=smooth, axis=axis)
+        sp = CubicSmoothingSpline(xdata, ydata, weights=weights, smooth=smooth, axis=axis,
+                                  normalizedsmooth=normalizedsmooth)
     else:
-        sp = NdGridCubicSmoothingSpline(xdata, ydata, weights, smooth)
+        sp = NdGridCubicSmoothingSpline(xdata, ydata, weights, smooth, normalizedsmooth=normalizedsmooth)
 
     if xidata is None:
         return sp

csaps/_sspndg.py

@@ -193,6 +193,10 @@ class NdGridCubicSmoothingSpline(ISmoothingSpline[
             - 0: The smoothing spline is the least-squares straight line fit
             - 1: The cubic spline interpolant with natural condition
 
+    normalizedsmooth : [*Optional*] bool
+        If True, the smooth parameter is normalized such that results are invariant to xdata range
+        and less sensitive to nonuniformity of weights and xdata clumping
+
     """
 
     __module__ = 'csaps'
@@ -201,10 +205,11 @@ def __init__(self,
                  xdata: NdGridDataType,
                  ydata: np.ndarray,
                  weights: Optional[Union[UnivariateDataType, NdGridDataType]] = None,
-                 smooth: Optional[Union[float, Sequence[Optional[float]]]] = None) -> None:
+                 smooth: Optional[Union[float, Sequence[Optional[float]]]] = None,
+                 normalizedsmooth: bool = False) -> None:
 
         x, y, w, s = self._prepare_data(xdata, ydata, weights, smooth)
-        coeffs, smooth = self._make_spline(x, y, w, s)
+        coeffs, smooth = self._make_spline(x, y, w, s, normalizedsmooth)
 
         self._spline = NdGridSplinePPForm.construct_fast(coeffs, x)
         self._smooth = smooth
@@ -306,12 +311,12 @@ def _prepare_data(cls, xdata, ydata, weights, smooth):
         return xdata, ydata, weights, smooth
 
     @staticmethod
-    def _make_spline(xdata, ydata, weights, smooth):
+    def _make_spline(xdata, ydata, weights, smooth, normalizedsmooth):
         ndim = len(xdata)
 
         if ndim == 1:
             s = CubicSmoothingSpline(
-                xdata[0], ydata, weights=weights[0], smooth=smooth[0])
+                xdata[0], ydata, weights=weights[0], smooth=smooth[0], normalizedsmooth=normalizedsmooth)
             return s.spline.coeffs, (s.smooth,)
 
         shape = ydata.shape
@@ -327,7 +332,7 @@ def _make_spline(xdata, ydata, weights, smooth):
                 coeffs = coeffs.reshape(prod(coeffs.shape[:-1]), coeffs.shape[-1])
 
             s = CubicSmoothingSpline(
-                xdata[i], coeffs, weights=weights[i], smooth=smooth[i])
+                xdata[i], coeffs, weights=weights[i], smooth=smooth[i], normalizedsmooth=normalizedsmooth)
 
             smooths.append(s.smooth)
             coeffs = umv_coeffs_to_flatten(s.spline.coeffs)

csaps/_sspumv.py

@@ -114,6 +114,10 @@ class CubicSmoothingSpline(ISmoothingSpline[
         Axis along which ``ydata`` is assumed to be varying.
         Meaning that for x[i] the corresponding values are np.take(ydata, i, axis=axis).
         By default is -1 (the last axis).
+
+    normalizedsmooth : [*Optional*] bool
+        If True, the smooth parameter is normalized such that results are invariant to xdata range
+        and less sensitive to nonuniformity of weights and xdata clumping
     """
 
     __module__ = 'csaps'
@@ -123,10 +127,11 @@ def __init__(self,
                  ydata: MultivariateDataType,
                  weights: Optional[UnivariateDataType] = None,
                  smooth: Optional[float] = None,
-                 axis: int = -1):
+                 axis: int = -1,
+                 normalizedsmooth: bool = False):
 
         x, y, w, shape, axis = self._prepare_data(xdata, ydata, weights, axis)
-        coeffs, smooth = self._make_spline(x, y, w, smooth, shape)
+        coeffs, smooth = self._make_spline(x, y, w, smooth, shape, normalizedsmooth)
         spline = SplinePPForm.construct_fast(coeffs, x, axis=axis)
 
         self._smooth = smooth
@@ -236,7 +241,7 @@ def trace(m: sp.dia_matrix):
         return 1. / (1. + trace(a) / (6. * trace(b)))
 
     @staticmethod
-    def _make_spline(x, y, w, smooth, shape):
+    def _make_spline(x, y, w, smooth, shape, normalizedsmooth):
         pcount = x.size
         dx = np.diff(x)
 
@@ -270,10 +275,17 @@ def _make_spline(x, y, w, smooth, shape):
                sp.diags(diags_sqrw_recip, 0, (pcount, pcount)))
         qtw = qtw @ qtw.T
 
-        if smooth is None:
+        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)
+        elif smooth is None:
             p = CubicSmoothingSpline._compute_smooth(r, qtw)
-        else:
-            p = smooth
 
         pp = (6. * (1. - p))
 

tests/test_csaps.py

@@ -64,6 +64,29 @@ def test_shortcut_output(data, tolist):
     sp = csaps(x, y)
     assert isinstance(sp, sp_cls)
 
+    yi = csaps(x, y, xi, smooth=smooth, normalizedsmooth=True)
+    assert isinstance(yi, np.ndarray)
+
+    smoothed_data = csaps(x, y, xi, normalizedsmooth=True)
+    assert isinstance(smoothed_data, AutoSmoothingResult)
+
+    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),