[fix] Handle missing replia (#232)

* [fix] First version of a fix to cope with missing replica.

* [fix] added test for missing replica

* [fix] refactored fix for missing replica, modified tests

* [fix] refinement of tests

pyerrors/obs.py

@@ -1138,7 +1138,7 @@ def _intersection_idx(idl):
     return idinter
 
 
-def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
+def _expand_deltas_for_merge(deltas, idx, shape, new_idx, scalefactor):
     """Expand deltas defined on idx to the list of configs that is defined by new_idx.
        New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest
        common divisor of the step sizes is used as new step size.
@@ -1154,15 +1154,20 @@ def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
         Number of configs in idx.
     new_idx : list
         List of configs that defines the new range, has to be sorted in ascending order.
+    scalefactor : float
+        An additional scaling factor that can be applied to scale the fluctuations,
+        e.g., when Obs with differing numbers of replica are merged.
     """
-
     if type(idx) is range and type(new_idx) is range:
         if idx == new_idx:
-            return deltas
+            if scalefactor == 1:
+                return deltas
+            else:
+                return deltas * scalefactor
     ret = np.zeros(new_idx[-1] - new_idx[0] + 1)
     for i in range(shape):
         ret[idx[i] - new_idx[0]] = deltas[i]
-    return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx)
+    return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx) * scalefactor
 
 
 def derived_observable(func, data, array_mode=False, **kwargs):
@@ -1243,6 +1248,25 @@ def derived_observable(func, data, array_mode=False, **kwargs):
         new_r_values[name] = func(tmp_values, **kwargs)
         new_idl_d[name] = _merge_idx(idl)
 
+    def _compute_scalefactor_missing_rep(obs):
+        """
+        Computes the scale factor that is to be multiplied with the deltas
+        in the case where Obs with different subsets of replica are merged.
+        Returns a dictionary with the scale factor for each Monte Carlo name.
+
+        Parameters
+        ----------
+        obs : Obs
+            The observable corresponding to the deltas that are to be scaled
+        """
+        scalef_d = {}
+        for mc_name in obs.mc_names:
+            mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')]
+            new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')]
+            if len(mc_idl_d) > 0 and len(mc_idl_d) < len(new_mc_idl_d):
+                scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d])
+        return scalef_d
+
     if 'man_grad' in kwargs:
         deriv = np.asarray(kwargs.get('man_grad'))
         if new_values.shape + data.shape != deriv.shape:
@@ -1280,7 +1304,7 @@ def derived_observable(func, data, array_mode=False, **kwargs):
             d_extracted[name] = []
             ens_length = len(new_idl_d[name])
             for i_dat, dat in enumerate(data):
-                d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name]) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
+                d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name], _compute_scalefactor_missing_rep(o).get(name.split('|')[0], 1)) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
         for name in new_cov_names:
             g_extracted[name] = []
             zero_grad = _Zero_grad(new_covobs_lengths[name])
@@ -1302,11 +1326,12 @@ def derived_observable(func, data, array_mode=False, **kwargs):
                     new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
         else:
             for j_obs, obs in np.ndenumerate(data):
+                scalef_d = _compute_scalefactor_missing_rep(obs)
                 for name in obs.names:
                     if name in obs.cov_names:
                         new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
                     else:
-                        new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name])
+                        new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name], scalef_d.get(name.split('|')[0], 1))
 
         new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
 

tests/obs_test.py

@@ -5,6 +5,7 @@ import copy
 import matplotlib.pyplot as plt
 import pyerrors as pe
 import pytest
+import pyerrors.linalg
 from hypothesis import given, strategies as st
 
 np.random.seed(0)
@@ -1338,9 +1339,101 @@ def test_vec_gm():
     pe.gm(cc, S=4.12)
     assert np.all(np.vectorize(lambda x: x.S["qq"])(cc.content) == 4.12)
 
+
 def test_complex_addition():
     o = pe.pseudo_Obs(34.12, 1e-4, "testens")
     r = o + 2j
     assert r.real == o
     r = r * 1j
     assert r.imag == o
+
+
+def test_missing_replica():
+    N1 = 3000
+    N2 = 2000
+    O1 = np.random.normal(1.0, .1, N1 + N2)
+    O2 = .5 * O1[:N1]
+
+    w1 = N1 / (N1 + N2)
+    w2 = N2 / (N1 + N2)
+    m12 = np.mean(O1[N1:])
+    m2 = np.mean(O2)
+    d12 = np.std(O1[N1:]) / np.sqrt(N2)  # error of <O1> from second rep
+    d2 = np.std(O2) / np.sqrt(N1)  # error of <O2> from first rep
+    dval = np.sqrt((w2 * d12 / m2)**2 + (w2 * m12 * d2 / m2**2)**2)  # complete error of <O1>/<O2>
+
+    # pyerrors version that should give the same result
+    O1dobs = pe.Obs([O1[:N1], O1[N1:]], names=['E|1', 'E|2'])
+    O2dobs = pe.Obs([O2], names=['E|1'])
+    O1O2 = O1dobs / O2dobs
+    O1O2.gm(S=0)
+
+    # explicit construction with different ensembles
+    O1a = pe.Obs([O1[:N1]], names=['E|1'])
+    O1b = pe.Obs([O1[N1:]], names=['F|2'])
+    O1O2b = (w1 * O1a + w2 * O1b) / O2dobs
+    O1O2b.gm(S=0)
+
+    # pyerrors version without replica (missing configs)
+    O1c = pe.Obs([O1], names=['E|1'])
+    O1O2c = O1c / O2dobs
+    O1O2c.gm(S=0)
+
+    for o in [O1O2, O1O2b, O1O2c]:
+        assert(np.isclose(dval, o.dvalue, atol=0, rtol=5e-2))
+
+    o = O1O2 * O2dobs - O1dobs
+    o.gm()
+    assert(o.is_zero())
+
+    o = O1dobs / O1O2 - O2dobs
+    o.gm()
+    assert(o.is_zero())
+
+    # bring more randomness and complexity into the game
+    Nl = [int(np.random.uniform(low=500, high=5000)) for i in range(4)]
+    wl = np.array(Nl) / sum(Nl)
+    O1 = np.random.normal(1.0, .1, sum(Nl))
+
+    # pyerrors replica version
+    datl = [O1[:Nl[0]], O1[Nl[0]:sum(Nl[:2])], O1[sum(Nl[:2]):sum(Nl[:3])], O1[sum(Nl[:3]):sum(Nl[:4])]]
+    O1dobs = pe.Obs(datl, names=['E|%d' % (d) for d in range(len(Nl))])
+    O2dobs = .5 * pe.Obs([datl[0]], names=['E|0'])
+    O3dobs = 2. / pe.Obs([datl[1]], names=['E|1'])
+    O1O2 = O1dobs / O2dobs
+    O1O2.gm(S=0)
+    O1O2O3 = O1O2 * np.sinh(O3dobs)
+    O1O2O3.gm(S=0)
+
+    # explicit construction with different ensembles
+    charl = ['E', 'F', 'G', 'H']
+    Ol = [pe.Obs([datl[i]], names=['%s|%d' % (charl[i], i)]) for i in range(len(Nl))]
+    O1O2b = sum(np.array(Ol) * wl) / O2dobs
+    O1O2b.gm(S=0)
+    i = 1
+    O3dobsb = 2. / pe.Obs([datl[i]], names=['%s|%d' % (charl[i], i)])
+    O1O2O3b = O1O2b * np.sinh(O3dobsb)
+    O1O2O3b.gm(S=0)
+
+    for op in [[O1O2, O1O2b], [O1O2O3, O1O2O3b]]:
+        assert np.isclose(op[0].value, op[1].value)
+        assert np.isclose(op[0].dvalue, op[1].dvalue, atol=0, rtol=5e-2)
+
+    # perform the same test using the array_mode of derived_observable
+    O1O2 = pyerrors.linalg.matmul(np.diag(np.diag(np.reshape(4 * [O1dobs], (2, 2)))), np.diag(np.diag(np.reshape(4 * [1. / O2dobs], (2, 2)))))
+    O1O2O3 = pyerrors.linalg.matmul(O1O2, np.diag(np.diag(np.sinh(np.reshape(4 * [O3dobs], (2, 2))))))
+    O1O2 = O1O2[0][0]
+    O1O2.gm(S=0)
+    O1O2O3 = O1O2O3[0][0]
+    O1O2O3.gm(S=0)
+
+    O1O2b = pyerrors.linalg.matmul(np.diag(np.diag(np.reshape(4 * [sum(np.array(Ol) * wl)], (2, 2)))), np.diag(np.diag(np.reshape(4 * [1. / O2dobs], (2, 2)))))
+    O1O2O3b = pyerrors.linalg.matmul(O1O2b, np.diag(np.diag(np.sinh(np.reshape(4 * [O3dobsb], (2, 2))))))
+    O1O2b = O1O2b[0][0]
+    O1O2b.gm(S=0)
+    O1O2O3b = O1O2O3b[0][0]
+    O1O2O3b.gm(S=0)
+
+    for op in [[O1O2, O1O2b], [O1O2O3, O1O2O3b]]:
+        assert np.isclose(op[1].value, op[0].value)
+        assert np.isclose(op[1].dvalue, op[0].dvalue, atol=0, rtol=5e-2)