diff --git a/pywt/_thresholding.py b/pywt/_thresholding.py index 33af65b8..10a0127d 100644 --- a/pywt/_thresholding.py +++ b/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) diff --git a/pywt/tests/test_thresholding.py b/pywt/tests/test_thresholding.py index 8e1dc940..a4fdc945 100644 --- a/pywt/tests/test_thresholding.py +++ b/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())