Merge pull request #878 from Dmao233/cursor/fix-866-threshold-zero-nan-6b9b

BUG: return 0 not NaN when thresholding zeros at value 0

Author: rgommers

pywt/_thresholding.py

@@ -17,11 +17,14 @@ def soft(data, value, substitute=0):
     data = np.asarray(data)
     magnitude = np.absolute(data)
 
-    with np.errstate(divide='ignore'):
+    with np.errstate(divide='ignore', invalid='ignore'):
         # divide by zero okay as np.inf values get clipped, so ignore warning.
+        # 0/0 (value==0 and data==0) is invalid and yields NaN; zero those below.
         thresholded = (1 - value/magnitude)
         thresholded.clip(min=0, max=None, out=thresholded)
         thresholded = data * thresholded
+    # Keep dtype; sign(0) is undefined, but soft-threshold(0) is 0 for any λ≥0.
+    thresholded[magnitude == 0] = 0
 
     if substitute == 0:
         return thresholded
@@ -35,11 +38,13 @@ def nn_garrote(data, value, substitute=0):
     data = np.asarray(data)
     magnitude = np.absolute(data)
 
-    with np.errstate(divide='ignore'):
+    with np.errstate(divide='ignore', invalid='ignore'):
         # divide by zero okay as np.inf values get clipped, so ignore warning.
+        # 0/0 (value==0 and data==0) is invalid and yields NaN; zero those below.
         thresholded = (1 - value**2/magnitude**2)
         thresholded.clip(min=0, max=None, out=thresholded)
         thresholded = data * thresholded
+    thresholded[magnitude == 0] = 0
 
     if substitute == 0:
         return thresholded
@@ -235,12 +240,14 @@ def threshold_firm(data, value_low, value_high):
 
     data = np.asarray(data)
     magnitude = np.absolute(data)
-    with np.errstate(divide='ignore'):
+    with np.errstate(divide='ignore', invalid='ignore'):
         # divide by zero okay as np.inf values get clipped, so ignore warning.
+        # 0/0 (value_low==0 and data==0) is invalid and yields NaN; zero those below.
         vdiff = value_high - value_low
         thresholded = value_high * (1 - value_low/magnitude) / vdiff
         thresholded.clip(min=0, max=None, out=thresholded)
         thresholded = data * thresholded
+    thresholded[magnitude == 0] = 0
 
     # restore hard-thresholding behavior for values > value_high
     large_vals = np.where(magnitude > value_high)

pywt/tests/test_thresholding.py

@@ -165,3 +165,38 @@ def test_threshold_firm():
         mt_abs_firm = np.abs(d_firm[mt])
         assert_(np.all(mt_abs_firm < np.abs(d_hard[mt])))
         assert_(np.all(mt_abs_firm > np.abs(d_soft[mt])))
+
+
+def test_threshold_zero_value_with_zeros():
+    # Issue 866: value==0 and exact-zero data used to yield NaN (0/0).
+    data = np.array([0.0, 1.0, -2.0])
+    expected = np.array([0.0, 1.0, -2.0])
+
+    assert_allclose(pywt.threshold(data, 0.0, 'soft'), expected, rtol=1e-12)
+    assert_allclose(pywt.threshold(data, 0.0, 'garrote'), expected, rtol=1e-12)
+    assert_allclose(pywt.threshold_firm(data, 0.0, 0.0), expected, rtol=1e-12)
+
+    # all zeros remain zeros (and must not warn: pytest treats warnings as errors)
+    zeros = np.zeros(8)
+    assert_allclose(pywt.threshold(zeros, 0.0, 'soft'), zeros, rtol=1e-12)
+    assert_allclose(pywt.threshold(zeros, 0.0, 'garrote'), zeros, rtol=1e-12)
+    assert_allclose(pywt.threshold_firm(zeros, 0.0, 0.0), zeros, rtol=1e-12)
+
+    # complex: same 0/0 path, result must be 0+0j not nan+nanj
+    cdata = np.array([0.0, 1.0 + 1.0j])
+    cexpected = np.array([0.0 + 0.0j, 1.0 + 1.0j])
+    assert_allclose(pywt.threshold(cdata, 0.0, 'soft'), cexpected, rtol=1e-12)
+    assert_allclose(pywt.threshold(cdata, 0.0, 'garrote'), cexpected, rtol=1e-12)
+    assert_allclose(pywt.threshold_firm(cdata, 0.0, 0.0), cexpected, rtol=1e-12)
+
+    for dtype in float_dtypes:
+        typed = np.asarray(data if dtype in real_dtypes else cdata, dtype=dtype)
+        out_soft = pywt.threshold(typed, 0.0, 'soft')
+        out_garrote = pywt.threshold(typed, 0.0, 'garrote')
+        out_firm = pywt.threshold_firm(typed, 0.0, 0.0)
+        assert_equal(out_soft.dtype, typed.dtype)
+        assert_equal(out_garrote.dtype, typed.dtype)
+        assert_equal(out_firm.dtype, typed.dtype)
+        assert_(not np.isnan(out_soft).any())
+        assert_(not np.isnan(out_garrote).any())
+        assert_(not np.isnan(out_firm).any())