diff --git a/pyerrors/__init__.py b/pyerrors/__init__.py index 4ef64352..c639c887 100644 --- a/pyerrors/__init__.py +++ b/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 diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py index 210f79fa..d152b05b 100644 --- a/pyerrors/linalg.py +++ b/pyerrors/linalg.py @@ -1,25 +1,27 @@ 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 - using automatic differentiation. + """Construct a derived Obs for a matrix valued function according to func(data, **kwargs) using automatic differentiation. Parameters ---------- - func -- 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]. - man_grad -- manually supply a list or an array which contains the jacobian - of func. Use cautiously, supplying the wrong derivative will - not be intercepted. + 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 + 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. """ data = np.asarray(data) @@ -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 = {} - for i_data in raveled_data: - for name in new_names: - tmp = i_data.shape.get(name) + is_merged = len(list(filter(lambda o: o.is_merged is True, raveled_data))) > 0 + reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0 + new_idl_d = {} + for name in new_names: + 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: - new_shape[name] = tmp - else: - if new_shape[name] != tmp: - raise Exception('Shapes of ensemble', name, 'do not match.') + idl.append(tmp) + 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]]]))) + 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_samples.append(new_deltas[name]) + new_idl = [] + 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]) - - final_result[i_val] = Obs(new_samples, new_names, means=new_means) + new_idl.append(filtered_idl_d[name]) + 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 diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 87f40b58..6a2b414a 100644 --- a/pyerrors/obs.py +++ b/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'):