Merge branch 'develop' into documentation

pyerrors/__init__.py

@@ -137,6 +137,7 @@ from .obs import *
 from .correlators import *
 from .fits import *
 from . import dirac
+from . import input
 from . import linalg
 from . import misc
 from . import mpm

pyerrors/linalg.py

@@ -1,23 +1,25 @@
 import numpy as np
 from autograd import jacobian
 import autograd.numpy as anp  # Thinly-wrapped numpy
-from .obs import derived_observable, CObs, Obs
+from .obs import derived_observable, CObs, Obs, _merge_idx, _expand_deltas_for_merge, _filter_zeroes
 
 from functools import partial
 from autograd.extend import defvjp
 
 
 def derived_array(func, data, **kwargs):
-    """Construct a derived Obs according to func(data, **kwargs) of matrix value data
+    """Construct a derived Obs for a matrix valued function according to func(data, **kwargs) using automatic differentiation.
-    using automatic differentiation.
 
     Parameters
     ----------
-    func -- arbitrary function of the form func(data, **kwargs). For the
+    func : object
+        arbitrary function of the form func(data, **kwargs). For the
         automatic differentiation to work, all numpy functions have to have
         the autograd wrapper (use 'import autograd.numpy as anp').
-    data -- list of Obs, e.g. [obs1, obs2, obs3].
+    data : list
-    man_grad -- manually supply a list or an array which contains the jacobian
+        list of Obs, e.g. [obs1, obs2, obs3].
+    man_grad : list
+        manually supply a list or an array which contains the jacobian
         of func. Use cautiously, supplying the wrong derivative will
         not be intercepted.
     """
@@ -30,25 +32,29 @@ def derived_array(func, data, **kwargs):
         if isinstance(i_data, Obs):
             first_name = i_data.names[0]
             first_shape = i_data.shape[first_name]
+            first_idl = i_data.idl[first_name]
             break
 
     for i in range(len(raveled_data)):
         if isinstance(raveled_data[i], (int, float)):
-            raveled_data[i] = Obs([raveled_data[i] + np.zeros(first_shape)], [first_name])
+            raveled_data[i] = Obs([raveled_data[i] + np.zeros(first_shape)], [first_name], idl=[first_idl])
 
     n_obs = len(raveled_data)
     new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
 
-    new_shape = {}
+    is_merged = len(list(filter(lambda o: o.is_merged is True, raveled_data))) > 0
-    for i_data in raveled_data:
+    reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
+    new_idl_d = {}
     for name in new_names:
-            tmp = i_data.shape.get(name)
+        idl = []
+        for i_data in raveled_data:
+            tmp = i_data.idl.get(name)
             if tmp is not None:
-                if new_shape.get(name) is None:
+                idl.append(tmp)
-                    new_shape[name] = tmp
+        new_idl_d[name] = _merge_idx(idl)
-                else:
+        if not is_merged:
-                    if new_shape[name] != tmp:
+            is_merged = (1 != len(set([len(idx) for idx in [*idl, new_idl_d[name]]])))
-                        raise Exception('Shapes of ensemble', name, 'do not match.')
+
     if data.ndim == 1:
         values = np.array([o.value for o in data])
     else:
@@ -71,8 +77,6 @@ def derived_array(func, data, **kwargs):
         deriv = np.asarray(kwargs.get('man_grad'))
         if new_values.shape + data.shape != deriv.shape:
             raise Exception('Manual derivative does not have correct shape.')
-    elif kwargs.get('num_grad') is True:
-        raise Exception('Multi mode currently not supported for numerical derivative')
     else:
         deriv = jacobian(func)(values, **kwargs)
 
@@ -83,7 +87,7 @@ def derived_array(func, data, **kwargs):
         d_extracted[name] = []
         for i_dat, dat in enumerate(data):
             ens_length = dat.ravel()[0].shape[name]
-            d_extracted[name].append(np.array([o.deltas[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[name], o.idl[name], o.shape[name], new_idl_d[name]) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
 
     for i_val, new_val in np.ndenumerate(new_values):
         new_deltas = {}
@@ -95,12 +99,21 @@ def derived_array(func, data, **kwargs):
 
         new_samples = []
         new_means = []
-        for name in new_names:
+        new_idl = []
-            new_samples.append(new_deltas[name])
+        if is_merged:
+            filtered_names, filtered_deltas, filtered_idl_d = _filter_zeroes(new_names, new_deltas, new_idl_d)
+        else:
+            filtered_names = new_names
+            filtered_deltas = new_deltas
+            filtered_idl_d = new_idl_d
+        for name in filtered_names:
+            new_samples.append(filtered_deltas[name])
             new_means.append(new_r_values[name][i_val])
-
+            new_idl.append(filtered_idl_d[name])
-        final_result[i_val] = Obs(new_samples, new_names, means=new_means)
+        final_result[i_val] = Obs(new_samples, filtered_names, means=new_means, idl=new_idl)
         final_result[i_val]._value = new_val
+        final_result[i_val].is_merged = is_merged
+        final_result[i_val].reweighted = reweighted
 
     return final_result
 

pyerrors/obs.py

@@ -829,7 +829,7 @@ class CObs:
         return 'CObs[' + str(self) + ']'
 
 
-def merge_idx(idl):
+def _merge_idx(idl):
     """Returns the union of all lists in idl
 
     Parameters
@@ -856,7 +856,7 @@ def merge_idx(idl):
     return list(set().union(*idl))
 
 
-def expand_deltas_for_merge(deltas, idx, shape, new_idx):
+def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
     """Expand deltas defined on idx to the list of configs that is defined by new_idx.
        New, empy 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.
@@ -883,7 +883,7 @@ 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(names, deltas, idl, eps=Obs.filter_eps):
+def _filter_zeroes(names, deltas, idl, 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 names, deltas and
        idl according to the filtering.
@@ -976,7 +976,7 @@ def derived_observable(func, data, **kwargs):
             tmp = i_data.idl.get(name)
             if tmp is not None:
                 idl.append(tmp)
-        new_idl_d[name] = merge_idx(idl)
+        new_idl_d[name] = _merge_idx(idl)
         if not is_merged:
             is_merged = (1 != len(set([len(idx) for idx in [*idl, new_idl_d[name]]])))
 
@@ -1038,13 +1038,13 @@ def derived_observable(func, data, **kwargs):
         new_deltas = {}
         for j_obs, obs in np.ndenumerate(data):
             for name in obs.names:
-                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])
 
         new_samples = []
         new_means = []
         new_idl = []
         if is_merged:
-            filtered_names, filtered_deltas, filtered_idl_d = filter_zeroes(new_names, new_deltas, new_idl_d)
+            filtered_names, filtered_deltas, filtered_idl_d = _filter_zeroes(new_names, new_deltas, new_idl_d)
         else:
             filtered_names = new_names
             filtered_deltas = new_deltas
@@ -1064,7 +1064,7 @@ def derived_observable(func, data, **kwargs):
     return final_result
 
 
-def reduce_deltas(deltas, idx_old, idx_new):
+def _reduce_deltas(deltas, idx_old, idx_new):
     """Extract deltas defined on idx_old on all configs of idx_new.
 
     Parameters
@@ -1094,7 +1094,7 @@ def reduce_deltas(deltas, idx_old, idx_new):
                 pos = j
                 break
         if pos < 0:
-            raise Exception('Error in reduce_deltas: Config %d not in idx_old' % (idx_new[i]))
+            raise Exception('Error in _reduce_deltas: Config %d not in idx_old' % (idx_new[i]))
         ret[i] = deltas[j]
     return np.array(ret)
 
@@ -1124,7 +1124,7 @@ def reweight(weight, obs, **kwargs):
         new_samples = []
         w_deltas = {}
         for name in sorted(weight.names):
-            w_deltas[name] = reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
+            w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
             new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
         tmp_obs = Obs(new_samples, sorted(weight.names), idl=[obs[i].idl[name] for name in sorted(weight.names)])
 
@@ -1200,7 +1200,7 @@ def covariance(obs1, obs2, correlation=False, **kwargs):
     for name in sorted(set(obs1.names + obs2.names)):
         if (obs1.shape.get(name) != obs2.shape.get(name)) and (obs1.shape.get(name) is not None) and (obs2.shape.get(name) is not None):
             raise Exception('Shapes of ensemble', name, 'do not fit')
-        if (1 != len(set([len(idx) for idx in [obs1.idl[name], obs2.idl[name], merge_idx([obs1.idl[name], obs2.idl[name]])]]))):
+        if (1 != len(set([len(idx) for idx in [obs1.idl[name], obs2.idl[name], _merge_idx([obs1.idl[name], obs2.idl[name]])]]))):
             raise Exception('Shapes of ensemble', name, 'do not fit')
 
     if not hasattr(obs1, 'e_names') or not hasattr(obs2, 'e_names'):
@@ -1306,7 +1306,7 @@ def covariance2(obs1, obs2, correlation=False, **kwargs):
         for r_name in obs1.e_content[e_name]:
             if r_name not in obs2.e_content[e_name]:
                 continue
-            idl_d[r_name] = merge_idx([obs1.idl[r_name], obs2.idl[r_name]])
+            idl_d[r_name] = _merge_idx([obs1.idl[r_name], obs2.idl[r_name]])
             if idl_d[r_name] is range:
                 r_length.append(len(idl_d[r_name]))
             else:
@@ -1380,7 +1380,7 @@ def covariance3(obs1, obs2, correlation=False, **kwargs):
     for name in sorted(set(obs1.names + obs2.names)):
         if (obs1.shape.get(name) != obs2.shape.get(name)) and (obs1.shape.get(name) is not None) and (obs2.shape.get(name) is not None):
             raise Exception('Shapes of ensemble', name, 'do not fit')
-        if (1 != len(set([len(idx) for idx in [obs1.idl[name], obs2.idl[name], merge_idx([obs1.idl[name], obs2.idl[name]])]]))):
+        if (1 != len(set([len(idx) for idx in [obs1.idl[name], obs2.idl[name], _merge_idx([obs1.idl[name], obs2.idl[name]])]]))):
             raise Exception('Shapes of ensemble', name, 'do not fit')
 
     if not hasattr(obs1, 'e_names') or not hasattr(obs2, 'e_names'):