Repository: wright-group/WrightTools

Issue:
BUG: Trim writes channel to all NaNs when no outliers are present
Example code:
```
# import
import numpy as np
from matplotlib import pyplot as plt
import WrightTools as wt
# create arrays
x = np.linspace(-3, 3, 31)[:, None]
y = np.linspace(-3, 3, 31)[None, :]
arr = np.exp(-1 * (x ** 2 + y ** 2))
# create damaged array
arr2 = arr.copy()
arr2[15,15] = 20
# create data object
d = wt.data.Data()
d.create_variable("x", values=x)
d.create_variable("y", values=y)
d.create_channel("original", arr)
d.create_channel("damaged", arr2)
d.create_channel("trimmed", arr2)
d.create_channel("broken", arr)
d.transform("x", "y")
# now trim
d.trimmed.trim([2,2]) # this one works as expected
d.broken.trim([2,2]) # this one will write the channel to all NaNs
# create figure
fig, gs = wt.artists.create_figure(cols=[1, 1, 1, 1])
for i in range(4):
    ax = plt.subplot(gs[i])
    ax.pcolor(d, channel=i)
    ax.set_title(d.channel_names[i])
# pretty up
ticks = [-2, 0, 2]
wt.artists.set_fig_labels(
    xlabel=d.axes[0].label, ylabel=d.axes[1].label, xticks=ticks, yticks=ticks
)
```
stdout from execution is
```
1 outliers removed
0 outliers removed
/home/darien/source/WrightTools/WrightTools/_dataset.py:373: RuntimeWarning: All-NaN slice encountered
  return np.nanmax(dataset[s])
/home/darien/source/WrightTools/WrightTools/_dataset.py:375: RuntimeWarning: All-NaN axis encountered
  self.attrs["max"] = np.nanmax(list(self.chunkwise(f).values()))
```
![image](https://user-images.githubusercontent.com/15792657/54384984-8c3c9d00-4663-11e9-971b-e911f95489d3.png)


Discussion:

Patch:
diff --git a/WrightTools/data/_channel.py b/WrightTools/data/_channel.py
index 91a04ed..03d4973 100644
--- a/WrightTools/data/_channel.py
+++ b/WrightTools/data/_channel.py
@@ -152,14 +152,17 @@ class Channel(Dataset):
 
         factor : number (optional)
             Tolerance factor.  Default is 3.
-        replace : {'nan', 'mean', number} (optional)
+        replace : {'nan', 'mean', 'exclusive_mean', number} (optional)
             Behavior of outlier replacement. Default is nan.
 
             nan
                 Outliers are replaced by numpy nans.
 
             mean
-                Outliers are replaced by the mean of its neighborhood.
+                Outliers are replaced by the mean of its neighborhood, including itself.
+
+            exclusive_mean
+                Outilers are replaced by the mean of its neighborhood, not including itself.
 
             number
                 Array becomes given number.
@@ -177,6 +180,7 @@ class Channel(Dataset):
         warnings.warn("trim", category=wt_exceptions.EntireDatasetInMemoryWarning)
         outliers = []
         means = []
+        ex_means = []
         # find outliers
         for idx in np.ndindex(self.shape):
             slices = []
@@ -186,26 +190,33 @@ class Channel(Dataset):
                 slices.append(slice(start, stop, 1))
             neighbors = self[slices]
             mean = np.nanmean(neighbors)
+            sum_ = np.nansum(neighbors)
             limit = np.nanstd(neighbors) * factor
             if np.abs(self[idx] - mean) > limit:
                 outliers.append(idx)
                 means.append(mean)
+                # Note, "- 1" is to exclude the point itself, which is not nan, in order
+                # to enter this if block, as `np.abs(nan - mean)` is nan, which would
+                # evaluate to False
+                ex_means.append((sum_ - self[idx]) / (np.sum(~np.isnan(neighbors)) - 1))
+
         # replace outliers
         i = tuple(zip(*outliers))
-        if replace == "nan":
-            arr = self[:]
-            arr[i] = np.nan
-            self[:] = arr
-        elif replace == "mean":
-            arr = self[:]
-            arr[i] = means
-            self[:] = arr
-        elif isinstance(replace, numbers.Number):
-            arr = self[:]
-            arr[i] = replace
-            self[:] = arr
-        else:
-            raise KeyError("replace must be one of {nan, mean} or some number")
+
+        if len(i) == 0:
+            if verbose:
+                print("No outliers found")
+            return []
+
+        replace = {"nan": np.nan, "mean": means, "exclusive_mean": ex_means}.get(replace, replace)
+
+        # This may someday be available in h5py directly, but seems that day is not yet.
+        # This is annoying because it is the only reason we hold the whole set in memory.
+        # KFS 2019-03-21
+        arr = self[:]
+        arr[i] = replace
+        self[:] = arr
+
         # finish
         if verbose:
             print("%i outliers removed" % len(outliers))


Test patch:
diff --git a/tests/data/trim.py b/tests/data/trim.py
index 2a94167..f096904 100644
--- a/tests/data/trim.py
+++ b/tests/data/trim.py
@@ -32,18 +32,23 @@ def test_trim_2Dgauss():
     d.create_channel("damaged1", arr2)
     d.create_channel("damaged2", arr2)
     d.create_channel("damaged3", arr2)
+    d.create_channel("damaged4", arr2)
     d.transform("x", "y")
     # trim
+    d.original.trim([2, 2], factor=2)
     d.damaged1.trim([2, 2], factor=2)
     d.damaged2.trim([2, 2], factor=2, replace="mean")
     d.damaged3.trim([2, 2], factor=2, replace=0.5)
+    d.damaged4.trim([2, 2], factor=2, replace="exclusive_mean")
     # now heal
     d.create_channel("healed_linear", d.damaged1[:])
     d.heal(channel="healed_linear", fill_value=0, method="linear")
     # check
-    assert np.allclose(d.original[:], d.healed_linear[:], rtol=1e-1, atol=1e-1)
-    assert np.allclose(d.original[:], d.damaged2[:], rtol=1e-1, atol=9e-1)
-    assert np.allclose(d.original[:], d.damaged3[:], rtol=1e-1, atol=5e-1)
+    np.testing.assert_allclose(d.original[:], d.original[:], rtol=1e-1, atol=1e-1)
+    np.testing.assert_allclose(d.original[:], d.healed_linear[:], rtol=1e-1, atol=1e-1)
+    np.testing.assert_allclose(d.original[:], d.damaged2[:], rtol=1e-1, atol=9e-1)
+    np.testing.assert_allclose(d.original[:], d.damaged3[:], rtol=1e-1, atol=5e-1)
+    np.testing.assert_allclose(d.original[:], d.damaged4[:], rtol=1e-1, atol=3e-1)
 
 
 def test_trim_3Dgauss():
@@ -67,7 +72,7 @@ def test_trim_3Dgauss():
     # trim
     d.damaged.trim([2, 2, 2], factor=2, replace="mean")
     # check
-    assert np.allclose(d.original[:], d.damaged[:], rtol=1e-1, atol=9e-1)
+    np.testing.assert_allclose(d.original[:], d.damaged[:], rtol=1e-1, atol=9e-1)
 
 
 if __name__ == "__main__":