Repository: fjosw/pyerrors

Issue:
Issues with _filter_zeroes and Corr
I have noticed issues with `_filter_zeroes`. These appear when two observables with vastly different variances, defined on different subsets of an ensemble, are combined (e.g. by renormalizing a noisy correlator). In this case it can happen that configurations that belong only to the observable with the small variance are filtered out by chance. This is something that I did not have in mind when designing the method.
When this does happen to elements of correlators, the export to json fails because the Obs of the Corr are not defined on the same set of configurations anymore.

One can construct a case artificially:
```python 
li = [-pe.pseudo_Obs(.2, .1, 'a', samples=10) for i in range(6)]
for i in range(len(li)):
    li[i].idl['a'] = range(1, 21, 2)
b = pe.pseudo_Obs(1, 1e-11, 'a', samples=30)
c *= b
print(c[2].idl)
print(c[3].idl)
{'a': [1, 3, 5, 7, 9, 10, 11, 13, 15, 17, 19]}
{'a': [1, 3, 4, 5, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17, 19, 20, 26, 30]}
```
This is certainly a rare event but I have observed it a number of times in the past and this is quite annoying. 

A possible fix would be to remove the automatic filtering altogether and instead to put the feature in a function that may be called by the user. In this case, exact cancellations would manifest themselves as vanishing deltas and would not be filtered out anymore. I tend to think that this is better to live with than silent removals of configurations from the obs.

If we like to keep the feature, I see a number of possible improvements (use mean instead of the median in `_filter_zeroes`, lower filter_eps, ...). 

Discussion:

Patch:
diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 0799842..7d11da9 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -57,7 +57,6 @@ class Obs:
     tau_exp_dict = {}
     N_sigma_global = 1.0
     N_sigma_dict = {}
-    filter_eps = 1e-10
 
     def __init__(self, samples, names, idl=None, **kwargs):
         """ Initialize Obs object.
@@ -1102,35 +1101,6 @@ def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
     return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))])
 
 
-def _filter_zeroes(deltas, idx, eps=Obs.filter_eps):
-    """Filter out all configurations with vanishing fluctuation such that they do not
-       contribute to the error estimate anymore. Returns the new deltas and
-       idx according to the filtering.
-       A fluctuation is considered to be vanishing, if it is smaller than eps times
-       the mean of the absolute values of all deltas in one list.
-
-    Parameters
-    ----------
-    deltas : list
-        List of fluctuations
-    idx : list
-        List or ranges of configs on which the deltas are defined.
-    eps : float
-        Prefactor that enters the filter criterion.
-    """
-    new_deltas = []
-    new_idx = []
-    maxd = np.mean(np.fabs(deltas))
-    for i in range(len(deltas)):
-        if abs(deltas[i]) > eps * maxd:
-            new_deltas.append(deltas[i])
-            new_idx.append(idx[i])
-    if new_idx:
-        return np.array(new_deltas), new_idx
-    else:
-        return deltas, idx
-
-
 def derived_observable(func, data, array_mode=False, **kwargs):
     """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
 
@@ -1287,14 +1257,8 @@ def derived_observable(func, data, array_mode=False, **kwargs):
         new_names_obs = []
         for name in new_names:
             if name not in new_covobs:
-                if is_merged[name]:
-                    filtered_deltas, filtered_idl_d = _filter_zeroes(new_deltas[name], new_idl_d[name])
-                else:
-                    filtered_deltas = new_deltas[name]
-                    filtered_idl_d = new_idl_d[name]
-
-                new_samples.append(filtered_deltas)
-                new_idl.append(filtered_idl_d)
+                new_samples.append(new_deltas[name])
+                new_idl.append(new_idl_d[name])
                 new_means.append(new_r_values[name][i_val])
                 new_names_obs.append(name)
         final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)


Test patch:
diff --git a/tests/correlators_test.py b/tests/correlators_test.py
index 1ab31f1..0b408ba 100644
--- a/tests/correlators_test.py
+++ b/tests/correlators_test.py
@@ -542,3 +542,13 @@ def test_complex_Corr():
     print(ccorr)
     mcorr = pe.Corr(np.array([[ccorr, ccorr], [ccorr, ccorr]]))
     assert np.all([mcorr.imag[i] == -mcorr.real[i] for i in range(mcorr.T)])
+
+
+def test_corr_no_filtering():
+    li = [-pe.pseudo_Obs(.2, .1, 'a', samples=10) for i in range(96)]
+    for i in range(len(li)):
+        li[i].idl['a'] = range(1, 21, 2)
+    c= pe.Corr(li)
+    b = pe.pseudo_Obs(1, 1e-11, 'a', samples=30)
+    c *= b
+    assert np.all([c[0].idl == o.idl for o in c])