Weight threshold part duex (#319)

changes/319.bugfix.rst

@@ -0,0 +1 @@
+Update weight threshold calculation in outlier detection to work around numpy bug that introduces small numerical differences for a mean of a masked array.

src/stcal/outlier_detection/utils.py

@@ -68,20 +68,15 @@ def compute_weight_threshold(weight, maskpt):
     float
         The weight threshold for this integration.
     '''
-    # necessary in order to assure that mask gets applied correctly
-    if hasattr(weight, '_mask'):
-        del weight._mask
-    mask_zero_weight = np.equal(weight, 0.)
-    mask_nans = np.isnan(weight)
-    # Combine the masks
-    weight_masked = np.ma.array(weight, mask=np.logical_or(
-        mask_zero_weight, mask_nans))
-    # Sigma-clip the unmasked data
-    weight_masked = sigma_clip(weight_masked, sigma=3, maxiters=5)
-    mean_weight = np.mean(weight_masked)
-    # Mask pixels where weight falls below maskpt percent
-    weight_threshold = mean_weight * maskpt
-    return weight_threshold
+    return np.mean(
+        sigma_clip(
+            weight[np.isfinite(weight) & (weight != 0)],
+            sigma=3,
+            maxiters=5,
+            masked=False,
+            copy=False,
+        ),
+    dtype='f8') * maskpt
 
 
 def _abs_deriv(array):

tests/outlier_detection/test_utils.py

@@ -16,6 +16,7 @@
     reproject,
     medfilt,
 )
+from stcal.testing_helpers import MemoryThreshold
 
 
 @pytest.mark.parametrize("shape,diff", [
@@ -81,6 +82,23 @@ def test_compute_weight_threshold_zeros():
     np.testing.assert_allclose(result, 21)
 
 
+def test_compute_weight_threshold_memory():
+    """Test that weight threshold function modifies
+    the weight array in place"""
+    arr = np.zeros([500, 500], dtype=np.float32)
+    arr[:250, :250] = 42
+    arr[10,10] = 0
+    arr[-10,-10] = np.nan
+
+    # buffer to account for memory overhead needs to be small enough
+    # to ensure that the array was not copied
+    fractional_memory_buffer = 0.9
+    expected_mem = int(arr.nbytes*fractional_memory_buffer)
+    with MemoryThreshold(str(expected_mem) + " B"):
+        result = compute_weight_threshold(arr, 0.5)
+    np.testing.assert_allclose(result, 21)
+
+
 def test_flag_crs():
     sci = np.zeros((10, 10), dtype=np.float32)
     err = np.ones_like(sci)