diff --git a/CHANGELOG.md b/CHANGELOG.md index 5ee48b4..dd56828 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/csaps/_sspumv.py b/csaps/_sspumv.py index 585705c..e59205e 100644 --- a/csaps/_sspumv.py +++ b/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) diff --git a/csaps/_version.py b/csaps/_version.py index 39a1e61..a4a639a 100644 --- a/csaps/_version.py +++ b/csaps/_version.py @@ -1,3 +1,3 @@ # -*- coding: utf-8 -*- -__version__ = '1.0.4' +__version__ = '1.1.0.dev0' diff --git a/tests/test_csaps.py b/tests/test_csaps.py index bdb9415..3fbb4df 100644 --- a/tests/test_csaps.py +++ b/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)