diff --git a/docs/pyerrors/correlators.html b/docs/pyerrors/correlators.html index de97c0d9..c9c25c2e 100644 --- a/docs/pyerrors/correlators.html +++ b/docs/pyerrors/correlators.html @@ -1238,349 +1238,347 @@ 998 content_string += "Description: " + self.tag + "\n" 999 if self.N != 1: 1000 return content_string -1001 if isinstance(self[0], CObs): -1002 return content_string -1003 -1004 if print_range[1]: -1005 print_range[1] += 1 -1006 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' -1007 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): -1008 if sub_corr is None: -1009 content_string += str(i + print_range[0]) + '\n' -1010 else: -1011 content_string += str(i + print_range[0]) -1012 for element in sub_corr: -1013 content_string += '\t' + ' ' * int(element >= 0) + str(element) -1014 content_string += '\n' -1015 return content_string -1016 -1017 def __str__(self): -1018 return self.__repr__() -1019 -1020 # We define the basic operations, that can be performed with correlators. -1021 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. -1022 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. -1023 # One could try and tell Obs to check if the y in __mul__ is a Corr and -1024 -1025 def __add__(self, y): -1026 if isinstance(y, Corr): -1027 if ((self.N != y.N) or (self.T != y.T)): -1028 raise Exception("Addition of Corrs with different shape") -1029 newcontent = [] -1030 for t in range(self.T): -1031 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1032 newcontent.append(None) -1033 else: -1034 newcontent.append(self.content[t] + y.content[t]) -1035 return Corr(newcontent) -1036 -1037 elif isinstance(y, (Obs, int, float, CObs)): -1038 newcontent = [] -1039 for t in range(self.T): -1040 if _check_for_none(self, self.content[t]): -1041 newcontent.append(None) -1042 else: -1043 newcontent.append(self.content[t] + y) -1044 return Corr(newcontent, prange=self.prange) -1045 elif isinstance(y, np.ndarray): -1046 if y.shape == (self.T,): -1047 return Corr(list((np.array(self.content).T + y).T)) -1048 else: -1049 raise ValueError("operands could not be broadcast together") -1050 else: -1051 raise TypeError("Corr + wrong type") -1052 -1053 def __mul__(self, y): -1054 if isinstance(y, Corr): -1055 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1056 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1057 newcontent = [] -1058 for t in range(self.T): -1059 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1060 newcontent.append(None) -1061 else: -1062 newcontent.append(self.content[t] * y.content[t]) -1063 return Corr(newcontent) -1064 -1065 elif isinstance(y, (Obs, int, float, CObs)): -1066 newcontent = [] -1067 for t in range(self.T): -1068 if _check_for_none(self, self.content[t]): -1069 newcontent.append(None) -1070 else: -1071 newcontent.append(self.content[t] * y) -1072 return Corr(newcontent, prange=self.prange) -1073 elif isinstance(y, np.ndarray): -1074 if y.shape == (self.T,): -1075 return Corr(list((np.array(self.content).T * y).T)) -1076 else: -1077 raise ValueError("operands could not be broadcast together") -1078 else: -1079 raise TypeError("Corr * wrong type") -1080 -1081 def __truediv__(self, y): -1082 if isinstance(y, Corr): -1083 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1084 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1085 newcontent = [] -1086 for t in range(self.T): -1087 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1088 newcontent.append(None) -1089 else: -1090 newcontent.append(self.content[t] / y.content[t]) -1091 for t in range(self.T): -1092 if _check_for_none(self, newcontent[t]): -1093 continue -1094 if np.isnan(np.sum(newcontent[t]).value): -1095 newcontent[t] = None -1096 -1097 if all([item is None for item in newcontent]): -1098 raise Exception("Division returns completely undefined correlator") -1099 return Corr(newcontent) -1100 -1101 elif isinstance(y, (Obs, CObs)): -1102 if isinstance(y, Obs): -1103 if y.value == 0: -1104 raise Exception('Division by zero will return undefined correlator') -1105 if isinstance(y, CObs): -1106 if y.is_zero(): -1107 raise Exception('Division by zero will return undefined correlator') -1108 -1109 newcontent = [] -1110 for t in range(self.T): -1111 if _check_for_none(self, self.content[t]): -1112 newcontent.append(None) -1113 else: -1114 newcontent.append(self.content[t] / y) -1115 return Corr(newcontent, prange=self.prange) -1116 -1117 elif isinstance(y, (int, float)): -1118 if y == 0: -1119 raise Exception('Division by zero will return undefined correlator') -1120 newcontent = [] -1121 for t in range(self.T): -1122 if _check_for_none(self, self.content[t]): -1123 newcontent.append(None) -1124 else: -1125 newcontent.append(self.content[t] / y) -1126 return Corr(newcontent, prange=self.prange) -1127 elif isinstance(y, np.ndarray): -1128 if y.shape == (self.T,): -1129 return Corr(list((np.array(self.content).T / y).T)) -1130 else: -1131 raise ValueError("operands could not be broadcast together") -1132 else: -1133 raise TypeError('Corr / wrong type') -1134 -1135 def __neg__(self): -1136 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] -1137 return Corr(newcontent, prange=self.prange) -1138 -1139 def __sub__(self, y): -1140 return self + (-y) -1141 -1142 def __pow__(self, y): -1143 if isinstance(y, (Obs, int, float, CObs)): -1144 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] -1145 return Corr(newcontent, prange=self.prange) -1146 else: -1147 raise TypeError('Type of exponent not supported') -1148 -1149 def __abs__(self): -1150 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] -1151 return Corr(newcontent, prange=self.prange) -1152 -1153 # The numpy functions: -1154 def sqrt(self): -1155 return self ** 0.5 -1156 -1157 def log(self): -1158 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] -1159 return Corr(newcontent, prange=self.prange) -1160 -1161 def exp(self): -1162 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] -1163 return Corr(newcontent, prange=self.prange) -1164 -1165 def _apply_func_to_corr(self, func): -1166 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] -1167 for t in range(self.T): -1168 if _check_for_none(self, newcontent[t]): -1169 continue -1170 tmp_sum = np.sum(newcontent[t]) -1171 if hasattr(tmp_sum, "value"): -1172 if np.isnan(tmp_sum.value): -1173 newcontent[t] = None -1174 if all([item is None for item in newcontent]): -1175 raise Exception('Operation returns undefined correlator') -1176 return Corr(newcontent) -1177 -1178 def sin(self): -1179 return self._apply_func_to_corr(np.sin) -1180 -1181 def cos(self): -1182 return self._apply_func_to_corr(np.cos) -1183 -1184 def tan(self): -1185 return self._apply_func_to_corr(np.tan) -1186 -1187 def sinh(self): -1188 return self._apply_func_to_corr(np.sinh) -1189 -1190 def cosh(self): -1191 return self._apply_func_to_corr(np.cosh) -1192 -1193 def tanh(self): -1194 return self._apply_func_to_corr(np.tanh) -1195 -1196 def arcsin(self): -1197 return self._apply_func_to_corr(np.arcsin) -1198 -1199 def arccos(self): -1200 return self._apply_func_to_corr(np.arccos) -1201 -1202 def arctan(self): -1203 return self._apply_func_to_corr(np.arctan) -1204 -1205 def arcsinh(self): -1206 return self._apply_func_to_corr(np.arcsinh) -1207 -1208 def arccosh(self): -1209 return self._apply_func_to_corr(np.arccosh) -1210 -1211 def arctanh(self): -1212 return self._apply_func_to_corr(np.arctanh) -1213 -1214 # Right hand side operations (require tweak in main module to work) -1215 def __radd__(self, y): -1216 return self + y -1217 -1218 def __rsub__(self, y): -1219 return -self + y -1220 -1221 def __rmul__(self, y): -1222 return self * y -1223 -1224 def __rtruediv__(self, y): -1225 return (self / y) ** (-1) -1226 -1227 @property -1228 def real(self): -1229 def return_real(obs_OR_cobs): -1230 if isinstance(obs_OR_cobs.flatten()[0], CObs): -1231 return np.vectorize(lambda x: x.real)(obs_OR_cobs) -1232 else: -1233 return obs_OR_cobs +1001 +1002 if print_range[1]: +1003 print_range[1] += 1 +1004 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' +1005 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): +1006 if sub_corr is None: +1007 content_string += str(i + print_range[0]) + '\n' +1008 else: +1009 content_string += str(i + print_range[0]) +1010 for element in sub_corr: +1011 content_string += f"\t{element:+2}" +1012 content_string += '\n' +1013 return content_string +1014 +1015 def __str__(self): +1016 return self.__repr__() +1017 +1018 # We define the basic operations, that can be performed with correlators. +1019 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. +1020 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. +1021 # One could try and tell Obs to check if the y in __mul__ is a Corr and +1022 +1023 def __add__(self, y): +1024 if isinstance(y, Corr): +1025 if ((self.N != y.N) or (self.T != y.T)): +1026 raise Exception("Addition of Corrs with different shape") +1027 newcontent = [] +1028 for t in range(self.T): +1029 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1030 newcontent.append(None) +1031 else: +1032 newcontent.append(self.content[t] + y.content[t]) +1033 return Corr(newcontent) +1034 +1035 elif isinstance(y, (Obs, int, float, CObs)): +1036 newcontent = [] +1037 for t in range(self.T): +1038 if _check_for_none(self, self.content[t]): +1039 newcontent.append(None) +1040 else: +1041 newcontent.append(self.content[t] + y) +1042 return Corr(newcontent, prange=self.prange) +1043 elif isinstance(y, np.ndarray): +1044 if y.shape == (self.T,): +1045 return Corr(list((np.array(self.content).T + y).T)) +1046 else: +1047 raise ValueError("operands could not be broadcast together") +1048 else: +1049 raise TypeError("Corr + wrong type") +1050 +1051 def __mul__(self, y): +1052 if isinstance(y, Corr): +1053 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1054 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1055 newcontent = [] +1056 for t in range(self.T): +1057 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1058 newcontent.append(None) +1059 else: +1060 newcontent.append(self.content[t] * y.content[t]) +1061 return Corr(newcontent) +1062 +1063 elif isinstance(y, (Obs, int, float, CObs)): +1064 newcontent = [] +1065 for t in range(self.T): +1066 if _check_for_none(self, self.content[t]): +1067 newcontent.append(None) +1068 else: +1069 newcontent.append(self.content[t] * y) +1070 return Corr(newcontent, prange=self.prange) +1071 elif isinstance(y, np.ndarray): +1072 if y.shape == (self.T,): +1073 return Corr(list((np.array(self.content).T * y).T)) +1074 else: +1075 raise ValueError("operands could not be broadcast together") +1076 else: +1077 raise TypeError("Corr * wrong type") +1078 +1079 def __truediv__(self, y): +1080 if isinstance(y, Corr): +1081 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1082 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1083 newcontent = [] +1084 for t in range(self.T): +1085 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1086 newcontent.append(None) +1087 else: +1088 newcontent.append(self.content[t] / y.content[t]) +1089 for t in range(self.T): +1090 if _check_for_none(self, newcontent[t]): +1091 continue +1092 if np.isnan(np.sum(newcontent[t]).value): +1093 newcontent[t] = None +1094 +1095 if all([item is None for item in newcontent]): +1096 raise Exception("Division returns completely undefined correlator") +1097 return Corr(newcontent) +1098 +1099 elif isinstance(y, (Obs, CObs)): +1100 if isinstance(y, Obs): +1101 if y.value == 0: +1102 raise Exception('Division by zero will return undefined correlator') +1103 if isinstance(y, CObs): +1104 if y.is_zero(): +1105 raise Exception('Division by zero will return undefined correlator') +1106 +1107 newcontent = [] +1108 for t in range(self.T): +1109 if _check_for_none(self, self.content[t]): +1110 newcontent.append(None) +1111 else: +1112 newcontent.append(self.content[t] / y) +1113 return Corr(newcontent, prange=self.prange) +1114 +1115 elif isinstance(y, (int, float)): +1116 if y == 0: +1117 raise Exception('Division by zero will return undefined correlator') +1118 newcontent = [] +1119 for t in range(self.T): +1120 if _check_for_none(self, self.content[t]): +1121 newcontent.append(None) +1122 else: +1123 newcontent.append(self.content[t] / y) +1124 return Corr(newcontent, prange=self.prange) +1125 elif isinstance(y, np.ndarray): +1126 if y.shape == (self.T,): +1127 return Corr(list((np.array(self.content).T / y).T)) +1128 else: +1129 raise ValueError("operands could not be broadcast together") +1130 else: +1131 raise TypeError('Corr / wrong type') +1132 +1133 def __neg__(self): +1134 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] +1135 return Corr(newcontent, prange=self.prange) +1136 +1137 def __sub__(self, y): +1138 return self + (-y) +1139 +1140 def __pow__(self, y): +1141 if isinstance(y, (Obs, int, float, CObs)): +1142 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] +1143 return Corr(newcontent, prange=self.prange) +1144 else: +1145 raise TypeError('Type of exponent not supported') +1146 +1147 def __abs__(self): +1148 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] +1149 return Corr(newcontent, prange=self.prange) +1150 +1151 # The numpy functions: +1152 def sqrt(self): +1153 return self ** 0.5 +1154 +1155 def log(self): +1156 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] +1157 return Corr(newcontent, prange=self.prange) +1158 +1159 def exp(self): +1160 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] +1161 return Corr(newcontent, prange=self.prange) +1162 +1163 def _apply_func_to_corr(self, func): +1164 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] +1165 for t in range(self.T): +1166 if _check_for_none(self, newcontent[t]): +1167 continue +1168 tmp_sum = np.sum(newcontent[t]) +1169 if hasattr(tmp_sum, "value"): +1170 if np.isnan(tmp_sum.value): +1171 newcontent[t] = None +1172 if all([item is None for item in newcontent]): +1173 raise Exception('Operation returns undefined correlator') +1174 return Corr(newcontent) +1175 +1176 def sin(self): +1177 return self._apply_func_to_corr(np.sin) +1178 +1179 def cos(self): +1180 return self._apply_func_to_corr(np.cos) +1181 +1182 def tan(self): +1183 return self._apply_func_to_corr(np.tan) +1184 +1185 def sinh(self): +1186 return self._apply_func_to_corr(np.sinh) +1187 +1188 def cosh(self): +1189 return self._apply_func_to_corr(np.cosh) +1190 +1191 def tanh(self): +1192 return self._apply_func_to_corr(np.tanh) +1193 +1194 def arcsin(self): +1195 return self._apply_func_to_corr(np.arcsin) +1196 +1197 def arccos(self): +1198 return self._apply_func_to_corr(np.arccos) +1199 +1200 def arctan(self): +1201 return self._apply_func_to_corr(np.arctan) +1202 +1203 def arcsinh(self): +1204 return self._apply_func_to_corr(np.arcsinh) +1205 +1206 def arccosh(self): +1207 return self._apply_func_to_corr(np.arccosh) +1208 +1209 def arctanh(self): +1210 return self._apply_func_to_corr(np.arctanh) +1211 +1212 # Right hand side operations (require tweak in main module to work) +1213 def __radd__(self, y): +1214 return self + y +1215 +1216 def __rsub__(self, y): +1217 return -self + y +1218 +1219 def __rmul__(self, y): +1220 return self * y +1221 +1222 def __rtruediv__(self, y): +1223 return (self / y) ** (-1) +1224 +1225 @property +1226 def real(self): +1227 def return_real(obs_OR_cobs): +1228 if isinstance(obs_OR_cobs.flatten()[0], CObs): +1229 return np.vectorize(lambda x: x.real)(obs_OR_cobs) +1230 else: +1231 return obs_OR_cobs +1232 +1233 return self._apply_func_to_corr(return_real) 1234 -1235 return self._apply_func_to_corr(return_real) -1236 -1237 @property -1238 def imag(self): -1239 def return_imag(obs_OR_cobs): -1240 if isinstance(obs_OR_cobs.flatten()[0], CObs): -1241 return np.vectorize(lambda x: x.imag)(obs_OR_cobs) -1242 else: -1243 return obs_OR_cobs * 0 # So it stays the right type +1235 @property +1236 def imag(self): +1237 def return_imag(obs_OR_cobs): +1238 if isinstance(obs_OR_cobs.flatten()[0], CObs): +1239 return np.vectorize(lambda x: x.imag)(obs_OR_cobs) +1240 else: +1241 return obs_OR_cobs * 0 # So it stays the right type +1242 +1243 return self._apply_func_to_corr(return_imag) 1244 -1245 return self._apply_func_to_corr(return_imag) -1246 -1247 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): -1248 r''' Project large correlation matrix to lowest states -1249 -1250 This method can be used to reduce the size of an (N x N) correlation matrix -1251 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise -1252 is still small. -1253 -1254 Parameters -1255 ---------- -1256 Ntrunc: int -1257 Rank of the target matrix. -1258 tproj: int -1259 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. -1260 The default value is 3. -1261 t0proj: int -1262 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly -1263 discouraged for O(a) improved theories, since the correctness of the procedure -1264 cannot be granted in this case. The default value is 2. -1265 basematrix : Corr -1266 Correlation matrix that is used to determine the eigenvectors of the -1267 lowest states based on a GEVP. basematrix is taken to be the Corr itself if -1268 is is not specified. -1269 -1270 Notes -1271 ----- -1272 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving -1273 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$ -1274 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the -1275 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via -1276 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large -1277 correlation matrix and to remove some noise that is added by irrelevant operators. -1278 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated -1279 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. -1280 ''' -1281 -1282 if self.N == 1: -1283 raise Exception('Method cannot be applied to one-dimensional correlators.') -1284 if basematrix is None: -1285 basematrix = self -1286 if Ntrunc >= basematrix.N: -1287 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) -1288 if basematrix.N != self.N: -1289 raise Exception('basematrix and targetmatrix have to be of the same size.') +1245 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): +1246 r''' Project large correlation matrix to lowest states +1247 +1248 This method can be used to reduce the size of an (N x N) correlation matrix +1249 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise +1250 is still small. +1251 +1252 Parameters +1253 ---------- +1254 Ntrunc: int +1255 Rank of the target matrix. +1256 tproj: int +1257 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. +1258 The default value is 3. +1259 t0proj: int +1260 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly +1261 discouraged for O(a) improved theories, since the correctness of the procedure +1262 cannot be granted in this case. The default value is 2. +1263 basematrix : Corr +1264 Correlation matrix that is used to determine the eigenvectors of the +1265 lowest states based on a GEVP. basematrix is taken to be the Corr itself if +1266 is is not specified. +1267 +1268 Notes +1269 ----- +1270 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving +1271 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$ +1272 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the +1273 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via +1274 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large +1275 correlation matrix and to remove some noise that is added by irrelevant operators. +1276 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated +1277 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. +1278 ''' +1279 +1280 if self.N == 1: +1281 raise Exception('Method cannot be applied to one-dimensional correlators.') +1282 if basematrix is None: +1283 basematrix = self +1284 if Ntrunc >= basematrix.N: +1285 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) +1286 if basematrix.N != self.N: +1287 raise Exception('basematrix and targetmatrix have to be of the same size.') +1288 +1289 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] 1290 -1291 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] -1292 -1293 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) -1294 rmat = [] -1295 for t in range(basematrix.T): -1296 for i in range(Ntrunc): -1297 for j in range(Ntrunc): -1298 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] -1299 rmat.append(np.copy(tmpmat)) -1300 -1301 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] -1302 return Corr(newcontent) -1303 -1304 -1305def _sort_vectors(vec_set, ts): -1306 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" -1307 reference_sorting = np.array(vec_set[ts]) -1308 N = reference_sorting.shape[0] -1309 sorted_vec_set = [] -1310 for t in range(len(vec_set)): -1311 if vec_set[t] is None: -1312 sorted_vec_set.append(None) -1313 elif not t == ts: -1314 perms = [list(o) for o in permutations([i for i in range(N)], N)] -1315 best_score = 0 -1316 for perm in perms: -1317 current_score = 1 -1318 for k in range(N): -1319 new_sorting = reference_sorting.copy() -1320 new_sorting[perm[k], :] = vec_set[t][k] -1321 current_score *= abs(np.linalg.det(new_sorting)) -1322 if current_score > best_score: -1323 best_score = current_score -1324 best_perm = perm -1325 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) -1326 else: -1327 sorted_vec_set.append(vec_set[t]) +1291 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) +1292 rmat = [] +1293 for t in range(basematrix.T): +1294 for i in range(Ntrunc): +1295 for j in range(Ntrunc): +1296 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] +1297 rmat.append(np.copy(tmpmat)) +1298 +1299 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] +1300 return Corr(newcontent) +1301 +1302 +1303def _sort_vectors(vec_set, ts): +1304 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" +1305 reference_sorting = np.array(vec_set[ts]) +1306 N = reference_sorting.shape[0] +1307 sorted_vec_set = [] +1308 for t in range(len(vec_set)): +1309 if vec_set[t] is None: +1310 sorted_vec_set.append(None) +1311 elif not t == ts: +1312 perms = [list(o) for o in permutations([i for i in range(N)], N)] +1313 best_score = 0 +1314 for perm in perms: +1315 current_score = 1 +1316 for k in range(N): +1317 new_sorting = reference_sorting.copy() +1318 new_sorting[perm[k], :] = vec_set[t][k] +1319 current_score *= abs(np.linalg.det(new_sorting)) +1320 if current_score > best_score: +1321 best_score = current_score +1322 best_perm = perm +1323 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) +1324 else: +1325 sorted_vec_set.append(vec_set[t]) +1326 +1327 return sorted_vec_set 1328 -1329 return sorted_vec_set -1330 -1331 -1332def _check_for_none(corr, entry): -1333 """Checks if entry for correlator corr is None""" -1334 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2 -1335 -1336 -1337def _GEVP_solver(Gt, G0): -1338 """Helper function for solving the GEVP and sorting the eigenvectors. -1339 -1340 The helper function assumes that both provided matrices are symmetric and -1341 only processes the lower triangular part of both matrices. In case the matrices -1342 are not symmetric the upper triangular parts are effectively discarded.""" -1343 return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1] +1329 +1330def _check_for_none(corr, entry): +1331 """Checks if entry for correlator corr is None""" +1332 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2 +1333 +1334 +1335def _GEVP_solver(Gt, G0): +1336 """Helper function for solving the GEVP and sorting the eigenvectors. +1337 +1338 The helper function assumes that both provided matrices are symmetric and +1339 only processes the lower triangular part of both matrices. In case the matrices +1340 are not symmetric the upper triangular parts are effectively discarded.""" +1341 return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1] @@ -2584,308 +2582,306 @@ 999 content_string += "Description: " + self.tag + "\n" 1000 if self.N != 1: 1001 return content_string -1002 if isinstance(self[0], CObs): -1003 return content_string -1004 -1005 if print_range[1]: -1006 print_range[1] += 1 -1007 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' -1008 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): -1009 if sub_corr is None: -1010 content_string += str(i + print_range[0]) + '\n' -1011 else: -1012 content_string += str(i + print_range[0]) -1013 for element in sub_corr: -1014 content_string += '\t' + ' ' * int(element >= 0) + str(element) -1015 content_string += '\n' -1016 return content_string -1017 -1018 def __str__(self): -1019 return self.__repr__() -1020 -1021 # We define the basic operations, that can be performed with correlators. -1022 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. -1023 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. -1024 # One could try and tell Obs to check if the y in __mul__ is a Corr and -1025 -1026 def __add__(self, y): -1027 if isinstance(y, Corr): -1028 if ((self.N != y.N) or (self.T != y.T)): -1029 raise Exception("Addition of Corrs with different shape") -1030 newcontent = [] -1031 for t in range(self.T): -1032 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1033 newcontent.append(None) -1034 else: -1035 newcontent.append(self.content[t] + y.content[t]) -1036 return Corr(newcontent) -1037 -1038 elif isinstance(y, (Obs, int, float, CObs)): -1039 newcontent = [] -1040 for t in range(self.T): -1041 if _check_for_none(self, self.content[t]): -1042 newcontent.append(None) -1043 else: -1044 newcontent.append(self.content[t] + y) -1045 return Corr(newcontent, prange=self.prange) -1046 elif isinstance(y, np.ndarray): -1047 if y.shape == (self.T,): -1048 return Corr(list((np.array(self.content).T + y).T)) -1049 else: -1050 raise ValueError("operands could not be broadcast together") -1051 else: -1052 raise TypeError("Corr + wrong type") -1053 -1054 def __mul__(self, y): -1055 if isinstance(y, Corr): -1056 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1057 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1058 newcontent = [] -1059 for t in range(self.T): -1060 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1061 newcontent.append(None) -1062 else: -1063 newcontent.append(self.content[t] * y.content[t]) -1064 return Corr(newcontent) -1065 -1066 elif isinstance(y, (Obs, int, float, CObs)): -1067 newcontent = [] -1068 for t in range(self.T): -1069 if _check_for_none(self, self.content[t]): -1070 newcontent.append(None) -1071 else: -1072 newcontent.append(self.content[t] * y) -1073 return Corr(newcontent, prange=self.prange) -1074 elif isinstance(y, np.ndarray): -1075 if y.shape == (self.T,): -1076 return Corr(list((np.array(self.content).T * y).T)) -1077 else: -1078 raise ValueError("operands could not be broadcast together") -1079 else: -1080 raise TypeError("Corr * wrong type") -1081 -1082 def __truediv__(self, y): -1083 if isinstance(y, Corr): -1084 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1085 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1086 newcontent = [] -1087 for t in range(self.T): -1088 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1089 newcontent.append(None) -1090 else: -1091 newcontent.append(self.content[t] / y.content[t]) -1092 for t in range(self.T): -1093 if _check_for_none(self, newcontent[t]): -1094 continue -1095 if np.isnan(np.sum(newcontent[t]).value): -1096 newcontent[t] = None -1097 -1098 if all([item is None for item in newcontent]): -1099 raise Exception("Division returns completely undefined correlator") -1100 return Corr(newcontent) -1101 -1102 elif isinstance(y, (Obs, CObs)): -1103 if isinstance(y, Obs): -1104 if y.value == 0: -1105 raise Exception('Division by zero will return undefined correlator') -1106 if isinstance(y, CObs): -1107 if y.is_zero(): -1108 raise Exception('Division by zero will return undefined correlator') -1109 -1110 newcontent = [] -1111 for t in range(self.T): -1112 if _check_for_none(self, self.content[t]): -1113 newcontent.append(None) -1114 else: -1115 newcontent.append(self.content[t] / y) -1116 return Corr(newcontent, prange=self.prange) -1117 -1118 elif isinstance(y, (int, float)): -1119 if y == 0: -1120 raise Exception('Division by zero will return undefined correlator') -1121 newcontent = [] -1122 for t in range(self.T): -1123 if _check_for_none(self, self.content[t]): -1124 newcontent.append(None) -1125 else: -1126 newcontent.append(self.content[t] / y) -1127 return Corr(newcontent, prange=self.prange) -1128 elif isinstance(y, np.ndarray): -1129 if y.shape == (self.T,): -1130 return Corr(list((np.array(self.content).T / y).T)) -1131 else: -1132 raise ValueError("operands could not be broadcast together") -1133 else: -1134 raise TypeError('Corr / wrong type') -1135 -1136 def __neg__(self): -1137 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] -1138 return Corr(newcontent, prange=self.prange) -1139 -1140 def __sub__(self, y): -1141 return self + (-y) -1142 -1143 def __pow__(self, y): -1144 if isinstance(y, (Obs, int, float, CObs)): -1145 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] -1146 return Corr(newcontent, prange=self.prange) -1147 else: -1148 raise TypeError('Type of exponent not supported') -1149 -1150 def __abs__(self): -1151 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] -1152 return Corr(newcontent, prange=self.prange) -1153 -1154 # The numpy functions: -1155 def sqrt(self): -1156 return self ** 0.5 -1157 -1158 def log(self): -1159 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] -1160 return Corr(newcontent, prange=self.prange) -1161 -1162 def exp(self): -1163 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] -1164 return Corr(newcontent, prange=self.prange) -1165 -1166 def _apply_func_to_corr(self, func): -1167 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] -1168 for t in range(self.T): -1169 if _check_for_none(self, newcontent[t]): -1170 continue -1171 tmp_sum = np.sum(newcontent[t]) -1172 if hasattr(tmp_sum, "value"): -1173 if np.isnan(tmp_sum.value): -1174 newcontent[t] = None -1175 if all([item is None for item in newcontent]): -1176 raise Exception('Operation returns undefined correlator') -1177 return Corr(newcontent) -1178 -1179 def sin(self): -1180 return self._apply_func_to_corr(np.sin) -1181 -1182 def cos(self): -1183 return self._apply_func_to_corr(np.cos) -1184 -1185 def tan(self): -1186 return self._apply_func_to_corr(np.tan) -1187 -1188 def sinh(self): -1189 return self._apply_func_to_corr(np.sinh) -1190 -1191 def cosh(self): -1192 return self._apply_func_to_corr(np.cosh) -1193 -1194 def tanh(self): -1195 return self._apply_func_to_corr(np.tanh) -1196 -1197 def arcsin(self): -1198 return self._apply_func_to_corr(np.arcsin) -1199 -1200 def arccos(self): -1201 return self._apply_func_to_corr(np.arccos) -1202 -1203 def arctan(self): -1204 return self._apply_func_to_corr(np.arctan) -1205 -1206 def arcsinh(self): -1207 return self._apply_func_to_corr(np.arcsinh) -1208 -1209 def arccosh(self): -1210 return self._apply_func_to_corr(np.arccosh) -1211 -1212 def arctanh(self): -1213 return self._apply_func_to_corr(np.arctanh) -1214 -1215 # Right hand side operations (require tweak in main module to work) -1216 def __radd__(self, y): -1217 return self + y -1218 -1219 def __rsub__(self, y): -1220 return -self + y -1221 -1222 def __rmul__(self, y): -1223 return self * y -1224 -1225 def __rtruediv__(self, y): -1226 return (self / y) ** (-1) -1227 -1228 @property -1229 def real(self): -1230 def return_real(obs_OR_cobs): -1231 if isinstance(obs_OR_cobs.flatten()[0], CObs): -1232 return np.vectorize(lambda x: x.real)(obs_OR_cobs) -1233 else: -1234 return obs_OR_cobs +1002 +1003 if print_range[1]: +1004 print_range[1] += 1 +1005 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' +1006 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): +1007 if sub_corr is None: +1008 content_string += str(i + print_range[0]) + '\n' +1009 else: +1010 content_string += str(i + print_range[0]) +1011 for element in sub_corr: +1012 content_string += f"\t{element:+2}" +1013 content_string += '\n' +1014 return content_string +1015 +1016 def __str__(self): +1017 return self.__repr__() +1018 +1019 # We define the basic operations, that can be performed with correlators. +1020 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. +1021 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. +1022 # One could try and tell Obs to check if the y in __mul__ is a Corr and +1023 +1024 def __add__(self, y): +1025 if isinstance(y, Corr): +1026 if ((self.N != y.N) or (self.T != y.T)): +1027 raise Exception("Addition of Corrs with different shape") +1028 newcontent = [] +1029 for t in range(self.T): +1030 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1031 newcontent.append(None) +1032 else: +1033 newcontent.append(self.content[t] + y.content[t]) +1034 return Corr(newcontent) +1035 +1036 elif isinstance(y, (Obs, int, float, CObs)): +1037 newcontent = [] +1038 for t in range(self.T): +1039 if _check_for_none(self, self.content[t]): +1040 newcontent.append(None) +1041 else: +1042 newcontent.append(self.content[t] + y) +1043 return Corr(newcontent, prange=self.prange) +1044 elif isinstance(y, np.ndarray): +1045 if y.shape == (self.T,): +1046 return Corr(list((np.array(self.content).T + y).T)) +1047 else: +1048 raise ValueError("operands could not be broadcast together") +1049 else: +1050 raise TypeError("Corr + wrong type") +1051 +1052 def __mul__(self, y): +1053 if isinstance(y, Corr): +1054 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1055 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1056 newcontent = [] +1057 for t in range(self.T): +1058 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1059 newcontent.append(None) +1060 else: +1061 newcontent.append(self.content[t] * y.content[t]) +1062 return Corr(newcontent) +1063 +1064 elif isinstance(y, (Obs, int, float, CObs)): +1065 newcontent = [] +1066 for t in range(self.T): +1067 if _check_for_none(self, self.content[t]): +1068 newcontent.append(None) +1069 else: +1070 newcontent.append(self.content[t] * y) +1071 return Corr(newcontent, prange=self.prange) +1072 elif isinstance(y, np.ndarray): +1073 if y.shape == (self.T,): +1074 return Corr(list((np.array(self.content).T * y).T)) +1075 else: +1076 raise ValueError("operands could not be broadcast together") +1077 else: +1078 raise TypeError("Corr * wrong type") +1079 +1080 def __truediv__(self, y): +1081 if isinstance(y, Corr): +1082 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1083 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1084 newcontent = [] +1085 for t in range(self.T): +1086 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1087 newcontent.append(None) +1088 else: +1089 newcontent.append(self.content[t] / y.content[t]) +1090 for t in range(self.T): +1091 if _check_for_none(self, newcontent[t]): +1092 continue +1093 if np.isnan(np.sum(newcontent[t]).value): +1094 newcontent[t] = None +1095 +1096 if all([item is None for item in newcontent]): +1097 raise Exception("Division returns completely undefined correlator") +1098 return Corr(newcontent) +1099 +1100 elif isinstance(y, (Obs, CObs)): +1101 if isinstance(y, Obs): +1102 if y.value == 0: +1103 raise Exception('Division by zero will return undefined correlator') +1104 if isinstance(y, CObs): +1105 if y.is_zero(): +1106 raise Exception('Division by zero will return undefined correlator') +1107 +1108 newcontent = [] +1109 for t in range(self.T): +1110 if _check_for_none(self, self.content[t]): +1111 newcontent.append(None) +1112 else: +1113 newcontent.append(self.content[t] / y) +1114 return Corr(newcontent, prange=self.prange) +1115 +1116 elif isinstance(y, (int, float)): +1117 if y == 0: +1118 raise Exception('Division by zero will return undefined correlator') +1119 newcontent = [] +1120 for t in range(self.T): +1121 if _check_for_none(self, self.content[t]): +1122 newcontent.append(None) +1123 else: +1124 newcontent.append(self.content[t] / y) +1125 return Corr(newcontent, prange=self.prange) +1126 elif isinstance(y, np.ndarray): +1127 if y.shape == (self.T,): +1128 return Corr(list((np.array(self.content).T / y).T)) +1129 else: +1130 raise ValueError("operands could not be broadcast together") +1131 else: +1132 raise TypeError('Corr / wrong type') +1133 +1134 def __neg__(self): +1135 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] +1136 return Corr(newcontent, prange=self.prange) +1137 +1138 def __sub__(self, y): +1139 return self + (-y) +1140 +1141 def __pow__(self, y): +1142 if isinstance(y, (Obs, int, float, CObs)): +1143 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] +1144 return Corr(newcontent, prange=self.prange) +1145 else: +1146 raise TypeError('Type of exponent not supported') +1147 +1148 def __abs__(self): +1149 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] +1150 return Corr(newcontent, prange=self.prange) +1151 +1152 # The numpy functions: +1153 def sqrt(self): +1154 return self ** 0.5 +1155 +1156 def log(self): +1157 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] +1158 return Corr(newcontent, prange=self.prange) +1159 +1160 def exp(self): +1161 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] +1162 return Corr(newcontent, prange=self.prange) +1163 +1164 def _apply_func_to_corr(self, func): +1165 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] +1166 for t in range(self.T): +1167 if _check_for_none(self, newcontent[t]): +1168 continue +1169 tmp_sum = np.sum(newcontent[t]) +1170 if hasattr(tmp_sum, "value"): +1171 if np.isnan(tmp_sum.value): +1172 newcontent[t] = None +1173 if all([item is None for item in newcontent]): +1174 raise Exception('Operation returns undefined correlator') +1175 return Corr(newcontent) +1176 +1177 def sin(self): +1178 return self._apply_func_to_corr(np.sin) +1179 +1180 def cos(self): +1181 return self._apply_func_to_corr(np.cos) +1182 +1183 def tan(self): +1184 return self._apply_func_to_corr(np.tan) +1185 +1186 def sinh(self): +1187 return self._apply_func_to_corr(np.sinh) +1188 +1189 def cosh(self): +1190 return self._apply_func_to_corr(np.cosh) +1191 +1192 def tanh(self): +1193 return self._apply_func_to_corr(np.tanh) +1194 +1195 def arcsin(self): +1196 return self._apply_func_to_corr(np.arcsin) +1197 +1198 def arccos(self): +1199 return self._apply_func_to_corr(np.arccos) +1200 +1201 def arctan(self): +1202 return self._apply_func_to_corr(np.arctan) +1203 +1204 def arcsinh(self): +1205 return self._apply_func_to_corr(np.arcsinh) +1206 +1207 def arccosh(self): +1208 return self._apply_func_to_corr(np.arccosh) +1209 +1210 def arctanh(self): +1211 return self._apply_func_to_corr(np.arctanh) +1212 +1213 # Right hand side operations (require tweak in main module to work) +1214 def __radd__(self, y): +1215 return self + y +1216 +1217 def __rsub__(self, y): +1218 return -self + y +1219 +1220 def __rmul__(self, y): +1221 return self * y +1222 +1223 def __rtruediv__(self, y): +1224 return (self / y) ** (-1) +1225 +1226 @property +1227 def real(self): +1228 def return_real(obs_OR_cobs): +1229 if isinstance(obs_OR_cobs.flatten()[0], CObs): +1230 return np.vectorize(lambda x: x.real)(obs_OR_cobs) +1231 else: +1232 return obs_OR_cobs +1233 +1234 return self._apply_func_to_corr(return_real) 1235 -1236 return self._apply_func_to_corr(return_real) -1237 -1238 @property -1239 def imag(self): -1240 def return_imag(obs_OR_cobs): -1241 if isinstance(obs_OR_cobs.flatten()[0], CObs): -1242 return np.vectorize(lambda x: x.imag)(obs_OR_cobs) -1243 else: -1244 return obs_OR_cobs * 0 # So it stays the right type +1236 @property +1237 def imag(self): +1238 def return_imag(obs_OR_cobs): +1239 if isinstance(obs_OR_cobs.flatten()[0], CObs): +1240 return np.vectorize(lambda x: x.imag)(obs_OR_cobs) +1241 else: +1242 return obs_OR_cobs * 0 # So it stays the right type +1243 +1244 return self._apply_func_to_corr(return_imag) 1245 -1246 return self._apply_func_to_corr(return_imag) -1247 -1248 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): -1249 r''' Project large correlation matrix to lowest states -1250 -1251 This method can be used to reduce the size of an (N x N) correlation matrix -1252 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise -1253 is still small. -1254 -1255 Parameters -1256 ---------- -1257 Ntrunc: int -1258 Rank of the target matrix. -1259 tproj: int -1260 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. -1261 The default value is 3. -1262 t0proj: int -1263 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly -1264 discouraged for O(a) improved theories, since the correctness of the procedure -1265 cannot be granted in this case. The default value is 2. -1266 basematrix : Corr -1267 Correlation matrix that is used to determine the eigenvectors of the -1268 lowest states based on a GEVP. basematrix is taken to be the Corr itself if -1269 is is not specified. -1270 -1271 Notes -1272 ----- -1273 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving -1274 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$ -1275 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the -1276 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via -1277 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large -1278 correlation matrix and to remove some noise that is added by irrelevant operators. -1279 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated -1280 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. -1281 ''' -1282 -1283 if self.N == 1: -1284 raise Exception('Method cannot be applied to one-dimensional correlators.') -1285 if basematrix is None: -1286 basematrix = self -1287 if Ntrunc >= basematrix.N: -1288 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) -1289 if basematrix.N != self.N: -1290 raise Exception('basematrix and targetmatrix have to be of the same size.') +1246 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): +1247 r''' Project large correlation matrix to lowest states +1248 +1249 This method can be used to reduce the size of an (N x N) correlation matrix +1250 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise +1251 is still small. +1252 +1253 Parameters +1254 ---------- +1255 Ntrunc: int +1256 Rank of the target matrix. +1257 tproj: int +1258 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. +1259 The default value is 3. +1260 t0proj: int +1261 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly +1262 discouraged for O(a) improved theories, since the correctness of the procedure +1263 cannot be granted in this case. The default value is 2. +1264 basematrix : Corr +1265 Correlation matrix that is used to determine the eigenvectors of the +1266 lowest states based on a GEVP. basematrix is taken to be the Corr itself if +1267 is is not specified. +1268 +1269 Notes +1270 ----- +1271 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving +1272 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$ +1273 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the +1274 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via +1275 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large +1276 correlation matrix and to remove some noise that is added by irrelevant operators. +1277 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated +1278 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. +1279 ''' +1280 +1281 if self.N == 1: +1282 raise Exception('Method cannot be applied to one-dimensional correlators.') +1283 if basematrix is None: +1284 basematrix = self +1285 if Ntrunc >= basematrix.N: +1286 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) +1287 if basematrix.N != self.N: +1288 raise Exception('basematrix and targetmatrix have to be of the same size.') +1289 +1290 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] 1291 -1292 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] -1293 -1294 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) -1295 rmat = [] -1296 for t in range(basematrix.T): -1297 for i in range(Ntrunc): -1298 for j in range(Ntrunc): -1299 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] -1300 rmat.append(np.copy(tmpmat)) -1301 -1302 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] -1303 return Corr(newcontent) +1292 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) +1293 rmat = [] +1294 for t in range(basematrix.T): +1295 for i in range(Ntrunc): +1296 for j in range(Ntrunc): +1297 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] +1298 rmat.append(np.copy(tmpmat)) +1299 +1300 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] +1301 return Corr(newcontent) @@ -4684,8 +4680,8 @@ specifies a custom path for the file (default '.') - 
1155    def sqrt(self):
-1156        return self ** 0.5
+            
1153    def sqrt(self):
+1154        return self ** 0.5
 

 
 
@@ -4703,9 +4699,9 @@ specifies a custom path for the file (default '.')
- 
1158    def log(self):
-1159        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
-1160        return Corr(newcontent, prange=self.prange)
+            
1156    def log(self):
+1157        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
+1158        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4723,9 +4719,9 @@ specifies a custom path for the file (default '.')
- 
1162    def exp(self):
-1163        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
-1164        return Corr(newcontent, prange=self.prange)
+            
1160    def exp(self):
+1161        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
+1162        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4743,8 +4739,8 @@ specifies a custom path for the file (default '.')
- 
1179    def sin(self):
-1180        return self._apply_func_to_corr(np.sin)
+            
1177    def sin(self):
+1178        return self._apply_func_to_corr(np.sin)
 

 
 
@@ -4762,8 +4758,8 @@ specifies a custom path for the file (default '.')
- 
1182    def cos(self):
-1183        return self._apply_func_to_corr(np.cos)
+            
1180    def cos(self):
+1181        return self._apply_func_to_corr(np.cos)
 

 
 
@@ -4781,8 +4777,8 @@ specifies a custom path for the file (default '.')
- 
1185    def tan(self):
-1186        return self._apply_func_to_corr(np.tan)
+            
1183    def tan(self):
+1184        return self._apply_func_to_corr(np.tan)
 

 
 
@@ -4800,8 +4796,8 @@ specifies a custom path for the file (default '.')
- 
1188    def sinh(self):
-1189        return self._apply_func_to_corr(np.sinh)
+            
1186    def sinh(self):
+1187        return self._apply_func_to_corr(np.sinh)
 

 
 
@@ -4819,8 +4815,8 @@ specifies a custom path for the file (default '.')
- 
1191    def cosh(self):
-1192        return self._apply_func_to_corr(np.cosh)
+            
1189    def cosh(self):
+1190        return self._apply_func_to_corr(np.cosh)
 

 
 
@@ -4838,8 +4834,8 @@ specifies a custom path for the file (default '.')
- 
1194    def tanh(self):
-1195        return self._apply_func_to_corr(np.tanh)
+            
1192    def tanh(self):
+1193        return self._apply_func_to_corr(np.tanh)
 

 
 
@@ -4857,8 +4853,8 @@ specifies a custom path for the file (default '.')
- 
1197    def arcsin(self):
-1198        return self._apply_func_to_corr(np.arcsin)
+            
1195    def arcsin(self):
+1196        return self._apply_func_to_corr(np.arcsin)
 

 
 
@@ -4876,8 +4872,8 @@ specifies a custom path for the file (default '.')
- 
1200    def arccos(self):
-1201        return self._apply_func_to_corr(np.arccos)
+            
1198    def arccos(self):
+1199        return self._apply_func_to_corr(np.arccos)
 

 
 
@@ -4895,8 +4891,8 @@ specifies a custom path for the file (default '.')
- 
1203    def arctan(self):
-1204        return self._apply_func_to_corr(np.arctan)
+            
1201    def arctan(self):
+1202        return self._apply_func_to_corr(np.arctan)
 

 
 
@@ -4914,8 +4910,8 @@ specifies a custom path for the file (default '.')
- 
1206    def arcsinh(self):
-1207        return self._apply_func_to_corr(np.arcsinh)
+            
1204    def arcsinh(self):
+1205        return self._apply_func_to_corr(np.arcsinh)
 

 
 
@@ -4933,8 +4929,8 @@ specifies a custom path for the file (default '.')
- 
1209    def arccosh(self):
-1210        return self._apply_func_to_corr(np.arccosh)
+            
1207    def arccosh(self):
+1208        return self._apply_func_to_corr(np.arccosh)
 

 
 
@@ -4952,8 +4948,8 @@ specifies a custom path for the file (default '.')
- 
1212    def arctanh(self):
-1213        return self._apply_func_to_corr(np.arctanh)
+            
1210    def arctanh(self):
+1211        return self._apply_func_to_corr(np.arctanh)
 

 
 
@@ -4993,62 +4989,62 @@ specifies a custom path for the file (default '.')
- 
1248    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1249        r''' Project large correlation matrix to lowest states
-1250
-1251        This method can be used to reduce the size of an (N x N) correlation matrix
-1252        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1253        is still small.
-1254
-1255        Parameters
-1256        ----------
-1257        Ntrunc: int
-1258            Rank of the target matrix.
-1259        tproj: int
-1260            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1261            The default value is 3.
-1262        t0proj: int
-1263            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1264            discouraged for O(a) improved theories, since the correctness of the procedure
-1265            cannot be granted in this case. The default value is 2.
-1266        basematrix : Corr
-1267            Correlation matrix that is used to determine the eigenvectors of the
-1268            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1269            is is not specified.
-1270
-1271        Notes
-1272        -----
-1273        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1274        the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
-1275        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1276        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1277        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1278        correlation matrix and to remove some noise that is added by irrelevant operators.
-1279        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1280        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1281        '''
-1282
-1283        if self.N == 1:
-1284            raise Exception('Method cannot be applied to one-dimensional correlators.')
-1285        if basematrix is None:
-1286            basematrix = self
-1287        if Ntrunc >= basematrix.N:
-1288            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
-1289        if basematrix.N != self.N:
-1290            raise Exception('basematrix and targetmatrix have to be of the same size.')
+            
1246    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1247        r''' Project large correlation matrix to lowest states
+1248
+1249        This method can be used to reduce the size of an (N x N) correlation matrix
+1250        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1251        is still small.
+1252
+1253        Parameters
+1254        ----------
+1255        Ntrunc: int
+1256            Rank of the target matrix.
+1257        tproj: int
+1258            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1259            The default value is 3.
+1260        t0proj: int
+1261            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1262            discouraged for O(a) improved theories, since the correctness of the procedure
+1263            cannot be granted in this case. The default value is 2.
+1264        basematrix : Corr
+1265            Correlation matrix that is used to determine the eigenvectors of the
+1266            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1267            is is not specified.
+1268
+1269        Notes
+1270        -----
+1271        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1272        the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
+1273        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1274        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1275        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1276        correlation matrix and to remove some noise that is added by irrelevant operators.
+1277        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1278        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1279        '''
+1280
+1281        if self.N == 1:
+1282            raise Exception('Method cannot be applied to one-dimensional correlators.')
+1283        if basematrix is None:
+1284            basematrix = self
+1285        if Ntrunc >= basematrix.N:
+1286            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
+1287        if basematrix.N != self.N:
+1288            raise Exception('basematrix and targetmatrix have to be of the same size.')
+1289
+1290        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
 1291
-1292        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
-1293
-1294        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1295        rmat = []
-1296        for t in range(basematrix.T):
-1297            for i in range(Ntrunc):
-1298                for j in range(Ntrunc):
-1299                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1300            rmat.append(np.copy(tmpmat))
-1301
-1302        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
-1303        return Corr(newcontent)
+1292        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1293        rmat = []
+1294        for t in range(basematrix.T):
+1295            for i in range(Ntrunc):
+1296                for j in range(Ntrunc):
+1297                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1298            rmat.append(np.copy(tmpmat))
+1299
+1300        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
+1301        return Corr(newcontent)
 

 
 
diff --git a/docs/pyerrors/obs.html b/docs/pyerrors/obs.html
index 313720b4..4bf2fc9f 100644
--- a/docs/pyerrors/obs.html
+++ b/docs/pyerrors/obs.html
@@ -1347,708 +1347,716 @@
 1023    def __repr__(self):
 1024        return 'CObs[' + str(self) + ']'
 1025
-1026
-1027def _format_uncertainty(value, dvalue, significance=2):
-1028    """Creates a string of a value and its error in paranthesis notation, e.g., 13.02(45)"""
-1029    if dvalue == 0.0 or (not np.isfinite(dvalue)):
-1030        return str(value)
-1031    if not isinstance(significance, int):
-1032        raise TypeError("significance needs to be an integer.")
-1033    if significance < 1:
-1034        raise ValueError("significance needs to be larger than zero.")
-1035    fexp = np.floor(np.log10(dvalue))
-1036    if fexp < 0.0:
-1037        return '{:{form}}({:1.0f})'.format(value, dvalue * 10 ** (-fexp + significance - 1), form='.' + str(-int(fexp) + significance - 1) + 'f')
-1038    elif fexp == 0.0:
-1039        return f"{value:.{significance - 1}f}({dvalue:1.{significance - 1}f})"
-1040    else:
-1041        return f"{value:.{max(0, int(significance - fexp - 1))}f}({dvalue:2.{max(0, int(significance - fexp - 1))}f})"
-1042
-1043
-1044def _expand_deltas(deltas, idx, shape, gapsize):
-1045    """Expand deltas defined on idx to a regular range with spacing gapsize between two
-1046       configurations and where holes are filled by 0.
-1047       If idx is of type range, the deltas are not changed if the idx.step == gapsize.
-1048
-1049    Parameters
-1050    ----------
-1051    deltas : list
-1052        List of fluctuations
-1053    idx : list
-1054        List or range of configs on which the deltas are defined, has to be sorted in ascending order.
-1055    shape : int
-1056        Number of configs in idx.
-1057    gapsize : int
-1058        The target distance between two configurations. If longer distances
-1059        are found in idx, the data is expanded.
-1060    """
-1061    if isinstance(idx, range):
-1062        if (idx.step == gapsize):
-1063            return deltas
-1064    ret = np.zeros((idx[-1] - idx[0] + gapsize) // gapsize)
-1065    for i in range(shape):
-1066        ret[(idx[i] - idx[0]) // gapsize] = deltas[i]
-1067    return ret
-1068
-1069
-1070def _merge_idx(idl):
-1071    """Returns the union of all lists in idl as range or sorted list
-1072
-1073    Parameters
-1074    ----------
-1075    idl : list
-1076        List of lists or ranges.
-1077    """
-1078
-1079    if _check_lists_equal(idl):
-1080        return idl[0]
-1081
-1082    idunion = sorted(set().union(*idl))
-1083
-1084    # Check whether idunion can be expressed as range
-1085    idrange = range(idunion[0], idunion[-1] + 1, idunion[1] - idunion[0])
-1086    idtest = [list(idrange), idunion]
-1087    if _check_lists_equal(idtest):
-1088        return idrange
+1026    def __format__(self, format_type):
+1027        if format_type == "":
+1028            significance = 2
+1029            format_type = "2"
+1030        else:
+1031            significance = int(float(format_type.replace("+", "").replace("-", "")))
+1032        return f"({self.real:{format_type}}{self.imag:+{significance}}j)"
+1033
+1034
+1035def _format_uncertainty(value, dvalue, significance=2):
+1036    """Creates a string of a value and its error in paranthesis notation, e.g., 13.02(45)"""
+1037    if dvalue == 0.0 or (not np.isfinite(dvalue)):
+1038        return str(value)
+1039    if not isinstance(significance, int):
+1040        raise TypeError("significance needs to be an integer.")
+1041    if significance < 1:
+1042        raise ValueError("significance needs to be larger than zero.")
+1043    fexp = np.floor(np.log10(dvalue))
+1044    if fexp < 0.0:
+1045        return '{:{form}}({:1.0f})'.format(value, dvalue * 10 ** (-fexp + significance - 1), form='.' + str(-int(fexp) + significance - 1) + 'f')
+1046    elif fexp == 0.0:
+1047        return f"{value:.{significance - 1}f}({dvalue:1.{significance - 1}f})"
+1048    else:
+1049        return f"{value:.{max(0, int(significance - fexp - 1))}f}({dvalue:2.{max(0, int(significance - fexp - 1))}f})"
+1050
+1051
+1052def _expand_deltas(deltas, idx, shape, gapsize):
+1053    """Expand deltas defined on idx to a regular range with spacing gapsize between two
+1054       configurations and where holes are filled by 0.
+1055       If idx is of type range, the deltas are not changed if the idx.step == gapsize.
+1056
+1057    Parameters
+1058    ----------
+1059    deltas : list
+1060        List of fluctuations
+1061    idx : list
+1062        List or range of configs on which the deltas are defined, has to be sorted in ascending order.
+1063    shape : int
+1064        Number of configs in idx.
+1065    gapsize : int
+1066        The target distance between two configurations. If longer distances
+1067        are found in idx, the data is expanded.
+1068    """
+1069    if isinstance(idx, range):
+1070        if (idx.step == gapsize):
+1071            return deltas
+1072    ret = np.zeros((idx[-1] - idx[0] + gapsize) // gapsize)
+1073    for i in range(shape):
+1074        ret[(idx[i] - idx[0]) // gapsize] = deltas[i]
+1075    return ret
+1076
+1077
+1078def _merge_idx(idl):
+1079    """Returns the union of all lists in idl as range or sorted list
+1080
+1081    Parameters
+1082    ----------
+1083    idl : list
+1084        List of lists or ranges.
+1085    """
+1086
+1087    if _check_lists_equal(idl):
+1088        return idl[0]
 1089
-1090    return idunion
+1090    idunion = sorted(set().union(*idl))
 1091
-1092
-1093def _intersection_idx(idl):
-1094    """Returns the intersection of all lists in idl as range or sorted list
-1095
-1096    Parameters
-1097    ----------
-1098    idl : list
-1099        List of lists or ranges.
-1100    """
-1101
-1102    if _check_lists_equal(idl):
-1103        return idl[0]
-1104
-1105    idinter = sorted(set.intersection(*[set(o) for o in idl]))
-1106
-1107    # Check whether idinter can be expressed as range
-1108    try:
-1109        idrange = range(idinter[0], idinter[-1] + 1, idinter[1] - idinter[0])
-1110        idtest = [list(idrange), idinter]
-1111        if _check_lists_equal(idtest):
-1112            return idrange
-1113    except IndexError:
-1114        pass
-1115
-1116    return idinter
-1117
-1118
-1119def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
-1120    """Expand deltas defined on idx to the list of configs that is defined by new_idx.
-1121       New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest
-1122       common divisor of the step sizes is used as new step size.
+1092    # Check whether idunion can be expressed as range
+1093    idrange = range(idunion[0], idunion[-1] + 1, idunion[1] - idunion[0])
+1094    idtest = [list(idrange), idunion]
+1095    if _check_lists_equal(idtest):
+1096        return idrange
+1097
+1098    return idunion
+1099
+1100
+1101def _intersection_idx(idl):
+1102    """Returns the intersection of all lists in idl as range or sorted list
+1103
+1104    Parameters
+1105    ----------
+1106    idl : list
+1107        List of lists or ranges.
+1108    """
+1109
+1110    if _check_lists_equal(idl):
+1111        return idl[0]
+1112
+1113    idinter = sorted(set.intersection(*[set(o) for o in idl]))
+1114
+1115    # Check whether idinter can be expressed as range
+1116    try:
+1117        idrange = range(idinter[0], idinter[-1] + 1, idinter[1] - idinter[0])
+1118        idtest = [list(idrange), idinter]
+1119        if _check_lists_equal(idtest):
+1120            return idrange
+1121    except IndexError:
+1122        pass
 1123
-1124    Parameters
-1125    ----------
-1126    deltas : list
-1127        List of fluctuations
-1128    idx : list
-1129        List or range of configs on which the deltas are defined.
-1130        Has to be a subset of new_idx and has to be sorted in ascending order.
-1131    shape : list
-1132        Number of configs in idx.
-1133    new_idx : list
-1134        List of configs that defines the new range, has to be sorted in ascending order.
-1135    """
-1136
-1137    if type(idx) is range and type(new_idx) is range:
-1138        if idx == new_idx:
-1139            return deltas
-1140    ret = np.zeros(new_idx[-1] - new_idx[0] + 1)
-1141    for i in range(shape):
-1142        ret[idx[i] - new_idx[0]] = deltas[i]
-1143    return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx)
+1124    return idinter
+1125
+1126
+1127def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
+1128    """Expand deltas defined on idx to the list of configs that is defined by new_idx.
+1129       New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest
+1130       common divisor of the step sizes is used as new step size.
+1131
+1132    Parameters
+1133    ----------
+1134    deltas : list
+1135        List of fluctuations
+1136    idx : list
+1137        List or range of configs on which the deltas are defined.
+1138        Has to be a subset of new_idx and has to be sorted in ascending order.
+1139    shape : list
+1140        Number of configs in idx.
+1141    new_idx : list
+1142        List of configs that defines the new range, has to be sorted in ascending order.
+1143    """
 1144
-1145
-1146def derived_observable(func, data, array_mode=False, **kwargs):
-1147    """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
-1148
-1149    Parameters
-1150    ----------
-1151    func : object
-1152        arbitrary function of the form func(data, **kwargs). For the
-1153        automatic differentiation to work, all numpy functions have to have
-1154        the autograd wrapper (use 'import autograd.numpy as anp').
-1155    data : list
-1156        list of Obs, e.g. [obs1, obs2, obs3].
-1157    num_grad : bool
-1158        if True, numerical derivatives are used instead of autograd
-1159        (default False). To control the numerical differentiation the
-1160        kwargs of numdifftools.step_generators.MaxStepGenerator
-1161        can be used.
-1162    man_grad : list
-1163        manually supply a list or an array which contains the jacobian
-1164        of func. Use cautiously, supplying the wrong derivative will
-1165        not be intercepted.
-1166
-1167    Notes
-1168    -----
-1169    For simple mathematical operations it can be practical to use anonymous
-1170    functions. For the ratio of two observables one can e.g. use
-1171
-1172    new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
-1173    """
+1145    if type(idx) is range and type(new_idx) is range:
+1146        if idx == new_idx:
+1147            return deltas
+1148    ret = np.zeros(new_idx[-1] - new_idx[0] + 1)
+1149    for i in range(shape):
+1150        ret[idx[i] - new_idx[0]] = deltas[i]
+1151    return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx)
+1152
+1153
+1154def derived_observable(func, data, array_mode=False, **kwargs):
+1155    """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
+1156
+1157    Parameters
+1158    ----------
+1159    func : object
+1160        arbitrary function of the form func(data, **kwargs). For the
+1161        automatic differentiation to work, all numpy functions have to have
+1162        the autograd wrapper (use 'import autograd.numpy as anp').
+1163    data : list
+1164        list of Obs, e.g. [obs1, obs2, obs3].
+1165    num_grad : bool
+1166        if True, numerical derivatives are used instead of autograd
+1167        (default False). To control the numerical differentiation the
+1168        kwargs of numdifftools.step_generators.MaxStepGenerator
+1169        can be used.
+1170    man_grad : list
+1171        manually supply a list or an array which contains the jacobian
+1172        of func. Use cautiously, supplying the wrong derivative will
+1173        not be intercepted.
 1174
-1175    data = np.asarray(data)
-1176    raveled_data = data.ravel()
-1177
-1178    # Workaround for matrix operations containing non Obs data
-1179    if not all(isinstance(x, Obs) for x in raveled_data):
-1180        for i in range(len(raveled_data)):
-1181            if isinstance(raveled_data[i], (int, float)):
-1182                raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
-1183
-1184    allcov = {}
-1185    for o in raveled_data:
-1186        for name in o.cov_names:
-1187            if name in allcov:
-1188                if not np.allclose(allcov[name], o.covobs[name].cov):
-1189                    raise Exception('Inconsistent covariance matrices for %s!' % (name))
-1190            else:
-1191                allcov[name] = o.covobs[name].cov
-1192
-1193    n_obs = len(raveled_data)
-1194    new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
-1195    new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
-1196    new_sample_names = sorted(set(new_names) - set(new_cov_names))
-1197
-1198    reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
-1199
-1200    if data.ndim == 1:
-1201        values = np.array([o.value for o in data])
-1202    else:
-1203        values = np.vectorize(lambda x: x.value)(data)
-1204
-1205    new_values = func(values, **kwargs)
-1206
-1207    multi = int(isinstance(new_values, np.ndarray))
-1208
-1209    new_r_values = {}
-1210    new_idl_d = {}
-1211    for name in new_sample_names:
-1212        idl = []
-1213        tmp_values = np.zeros(n_obs)
-1214        for i, item in enumerate(raveled_data):
-1215            tmp_values[i] = item.r_values.get(name, item.value)
-1216            tmp_idl = item.idl.get(name)
-1217            if tmp_idl is not None:
-1218                idl.append(tmp_idl)
-1219        if multi > 0:
-1220            tmp_values = np.array(tmp_values).reshape(data.shape)
-1221        new_r_values[name] = func(tmp_values, **kwargs)
-1222        new_idl_d[name] = _merge_idx(idl)
-1223
-1224    if 'man_grad' in kwargs:
-1225        deriv = np.asarray(kwargs.get('man_grad'))
-1226        if new_values.shape + data.shape != deriv.shape:
-1227            raise Exception('Manual derivative does not have correct shape.')
-1228    elif kwargs.get('num_grad') is True:
-1229        if multi > 0:
-1230            raise Exception('Multi mode currently not supported for numerical derivative')
-1231        options = {
-1232            'base_step': 0.1,
-1233            'step_ratio': 2.5}
-1234        for key in options.keys():
-1235            kwarg = kwargs.get(key)
-1236            if kwarg is not None:
-1237                options[key] = kwarg
-1238        tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
-1239        if tmp_df.size == 1:
-1240            deriv = np.array([tmp_df.real])
-1241        else:
-1242            deriv = tmp_df.real
-1243    else:
-1244        deriv = jacobian(func)(values, **kwargs)
-1245
-1246    final_result = np.zeros(new_values.shape, dtype=object)
-1247
-1248    if array_mode is True:
-1249
-1250        class _Zero_grad():
-1251            def __init__(self, N):
-1252                self.grad = np.zeros((N, 1))
+1175    Notes
+1176    -----
+1177    For simple mathematical operations it can be practical to use anonymous
+1178    functions. For the ratio of two observables one can e.g. use
+1179
+1180    new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
+1181    """
+1182
+1183    data = np.asarray(data)
+1184    raveled_data = data.ravel()
+1185
+1186    # Workaround for matrix operations containing non Obs data
+1187    if not all(isinstance(x, Obs) for x in raveled_data):
+1188        for i in range(len(raveled_data)):
+1189            if isinstance(raveled_data[i], (int, float)):
+1190                raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
+1191
+1192    allcov = {}
+1193    for o in raveled_data:
+1194        for name in o.cov_names:
+1195            if name in allcov:
+1196                if not np.allclose(allcov[name], o.covobs[name].cov):
+1197                    raise Exception('Inconsistent covariance matrices for %s!' % (name))
+1198            else:
+1199                allcov[name] = o.covobs[name].cov
+1200
+1201    n_obs = len(raveled_data)
+1202    new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
+1203    new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
+1204    new_sample_names = sorted(set(new_names) - set(new_cov_names))
+1205
+1206    reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
+1207
+1208    if data.ndim == 1:
+1209        values = np.array([o.value for o in data])
+1210    else:
+1211        values = np.vectorize(lambda x: x.value)(data)
+1212
+1213    new_values = func(values, **kwargs)
+1214
+1215    multi = int(isinstance(new_values, np.ndarray))
+1216
+1217    new_r_values = {}
+1218    new_idl_d = {}
+1219    for name in new_sample_names:
+1220        idl = []
+1221        tmp_values = np.zeros(n_obs)
+1222        for i, item in enumerate(raveled_data):
+1223            tmp_values[i] = item.r_values.get(name, item.value)
+1224            tmp_idl = item.idl.get(name)
+1225            if tmp_idl is not None:
+1226                idl.append(tmp_idl)
+1227        if multi > 0:
+1228            tmp_values = np.array(tmp_values).reshape(data.shape)
+1229        new_r_values[name] = func(tmp_values, **kwargs)
+1230        new_idl_d[name] = _merge_idx(idl)
+1231
+1232    if 'man_grad' in kwargs:
+1233        deriv = np.asarray(kwargs.get('man_grad'))
+1234        if new_values.shape + data.shape != deriv.shape:
+1235            raise Exception('Manual derivative does not have correct shape.')
+1236    elif kwargs.get('num_grad') is True:
+1237        if multi > 0:
+1238            raise Exception('Multi mode currently not supported for numerical derivative')
+1239        options = {
+1240            'base_step': 0.1,
+1241            'step_ratio': 2.5}
+1242        for key in options.keys():
+1243            kwarg = kwargs.get(key)
+1244            if kwarg is not None:
+1245                options[key] = kwarg
+1246        tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
+1247        if tmp_df.size == 1:
+1248            deriv = np.array([tmp_df.real])
+1249        else:
+1250            deriv = tmp_df.real
+1251    else:
+1252        deriv = jacobian(func)(values, **kwargs)
 1253
-1254        new_covobs_lengths = dict(set([y for x in [[(n, o.covobs[n].N) for n in o.cov_names] for o in raveled_data] for y in x]))
-1255        d_extracted = {}
-1256        g_extracted = {}
-1257        for name in new_sample_names:
-1258            d_extracted[name] = []
-1259            ens_length = len(new_idl_d[name])
-1260            for i_dat, dat in enumerate(data):
-1261                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, )))
-1262        for name in new_cov_names:
-1263            g_extracted[name] = []
-1264            zero_grad = _Zero_grad(new_covobs_lengths[name])
-1265            for i_dat, dat in enumerate(data):
-1266                g_extracted[name].append(np.array([o.covobs.get(name, zero_grad).grad for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (new_covobs_lengths[name], 1)))
-1267
-1268    for i_val, new_val in np.ndenumerate(new_values):
-1269        new_deltas = {}
-1270        new_grad = {}
-1271        if array_mode is True:
-1272            for name in new_sample_names:
-1273                ens_length = d_extracted[name][0].shape[-1]
-1274                new_deltas[name] = np.zeros(ens_length)
-1275                for i_dat, dat in enumerate(d_extracted[name]):
-1276                    new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1277            for name in new_cov_names:
-1278                new_grad[name] = 0
-1279                for i_dat, dat in enumerate(g_extracted[name]):
-1280                    new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1281        else:
-1282            for j_obs, obs in np.ndenumerate(data):
-1283                for name in obs.names:
-1284                    if name in obs.cov_names:
-1285                        new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
-1286                    else:
-1287                        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])
-1288
-1289        new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
-1290
-1291        if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
-1292            raise Exception('The same name has been used for deltas and covobs!')
-1293        new_samples = []
-1294        new_means = []
-1295        new_idl = []
-1296        new_names_obs = []
-1297        for name in new_names:
-1298            if name not in new_covobs:
-1299                new_samples.append(new_deltas[name])
-1300                new_idl.append(new_idl_d[name])
-1301                new_means.append(new_r_values[name][i_val])
-1302                new_names_obs.append(name)
-1303        final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
-1304        for name in new_covobs:
-1305            final_result[i_val].names.append(name)
-1306        final_result[i_val]._covobs = new_covobs
-1307        final_result[i_val]._value = new_val
-1308        final_result[i_val].reweighted = reweighted
-1309
-1310    if multi == 0:
-1311        final_result = final_result.item()
-1312
-1313    return final_result
-1314
-1315
-1316def _reduce_deltas(deltas, idx_old, idx_new):
-1317    """Extract deltas defined on idx_old on all configs of idx_new.
-1318
-1319    Assumes, that idx_old and idx_new are correctly defined idl, i.e., they
-1320    are ordered in an ascending order.
-1321
-1322    Parameters
-1323    ----------
-1324    deltas : list
-1325        List of fluctuations
-1326    idx_old : list
-1327        List or range of configs on which the deltas are defined
-1328    idx_new : list
-1329        List of configs for which we want to extract the deltas.
-1330        Has to be a subset of idx_old.
-1331    """
-1332    if not len(deltas) == len(idx_old):
-1333        raise Exception('Length of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old)))
-1334    if type(idx_old) is range and type(idx_new) is range:
-1335        if idx_old == idx_new:
-1336            return deltas
-1337    if _check_lists_equal([idx_old, idx_new]):
-1338        return deltas
-1339    indices = np.intersect1d(idx_old, idx_new, assume_unique=True, return_indices=True)[1]
-1340    if len(indices) < len(idx_new):
-1341        raise Exception('Error in _reduce_deltas: Config of idx_new not in idx_old')
-1342    return np.array(deltas)[indices]
-1343
-1344
-1345def reweight(weight, obs, **kwargs):
-1346    """Reweight a list of observables.
-1347
-1348    Parameters
-1349    ----------
-1350    weight : Obs
-1351        Reweighting factor. An Observable that has to be defined on a superset of the
-1352        configurations in obs[i].idl for all i.
-1353    obs : list
-1354        list of Obs, e.g. [obs1, obs2, obs3].
-1355    all_configs : bool
-1356        if True, the reweighted observables are normalized by the average of
-1357        the reweighting factor on all configurations in weight.idl and not
-1358        on the configurations in obs[i].idl. Default False.
-1359    """
-1360    result = []
-1361    for i in range(len(obs)):
-1362        if len(obs[i].cov_names):
-1363            raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
-1364        if not set(obs[i].names).issubset(weight.names):
-1365            raise Exception('Error: Ensembles do not fit')
-1366        for name in obs[i].names:
-1367            if not set(obs[i].idl[name]).issubset(weight.idl[name]):
-1368                raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
-1369        new_samples = []
-1370        w_deltas = {}
-1371        for name in sorted(obs[i].names):
-1372            w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
-1373            new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
-1374        tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
-1375
-1376        if kwargs.get('all_configs'):
-1377            new_weight = weight
-1378        else:
-1379            new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
-1380
-1381        result.append(tmp_obs / new_weight)
-1382        result[-1].reweighted = True
+1254    final_result = np.zeros(new_values.shape, dtype=object)
+1255
+1256    if array_mode is True:
+1257
+1258        class _Zero_grad():
+1259            def __init__(self, N):
+1260                self.grad = np.zeros((N, 1))
+1261
+1262        new_covobs_lengths = dict(set([y for x in [[(n, o.covobs[n].N) for n in o.cov_names] for o in raveled_data] for y in x]))
+1263        d_extracted = {}
+1264        g_extracted = {}
+1265        for name in new_sample_names:
+1266            d_extracted[name] = []
+1267            ens_length = len(new_idl_d[name])
+1268            for i_dat, dat in enumerate(data):
+1269                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, )))
+1270        for name in new_cov_names:
+1271            g_extracted[name] = []
+1272            zero_grad = _Zero_grad(new_covobs_lengths[name])
+1273            for i_dat, dat in enumerate(data):
+1274                g_extracted[name].append(np.array([o.covobs.get(name, zero_grad).grad for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (new_covobs_lengths[name], 1)))
+1275
+1276    for i_val, new_val in np.ndenumerate(new_values):
+1277        new_deltas = {}
+1278        new_grad = {}
+1279        if array_mode is True:
+1280            for name in new_sample_names:
+1281                ens_length = d_extracted[name][0].shape[-1]
+1282                new_deltas[name] = np.zeros(ens_length)
+1283                for i_dat, dat in enumerate(d_extracted[name]):
+1284                    new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1285            for name in new_cov_names:
+1286                new_grad[name] = 0
+1287                for i_dat, dat in enumerate(g_extracted[name]):
+1288                    new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1289        else:
+1290            for j_obs, obs in np.ndenumerate(data):
+1291                for name in obs.names:
+1292                    if name in obs.cov_names:
+1293                        new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
+1294                    else:
+1295                        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])
+1296
+1297        new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
+1298
+1299        if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
+1300            raise Exception('The same name has been used for deltas and covobs!')
+1301        new_samples = []
+1302        new_means = []
+1303        new_idl = []
+1304        new_names_obs = []
+1305        for name in new_names:
+1306            if name not in new_covobs:
+1307                new_samples.append(new_deltas[name])
+1308                new_idl.append(new_idl_d[name])
+1309                new_means.append(new_r_values[name][i_val])
+1310                new_names_obs.append(name)
+1311        final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
+1312        for name in new_covobs:
+1313            final_result[i_val].names.append(name)
+1314        final_result[i_val]._covobs = new_covobs
+1315        final_result[i_val]._value = new_val
+1316        final_result[i_val].reweighted = reweighted
+1317
+1318    if multi == 0:
+1319        final_result = final_result.item()
+1320
+1321    return final_result
+1322
+1323
+1324def _reduce_deltas(deltas, idx_old, idx_new):
+1325    """Extract deltas defined on idx_old on all configs of idx_new.
+1326
+1327    Assumes, that idx_old and idx_new are correctly defined idl, i.e., they
+1328    are ordered in an ascending order.
+1329
+1330    Parameters
+1331    ----------
+1332    deltas : list
+1333        List of fluctuations
+1334    idx_old : list
+1335        List or range of configs on which the deltas are defined
+1336    idx_new : list
+1337        List of configs for which we want to extract the deltas.
+1338        Has to be a subset of idx_old.
+1339    """
+1340    if not len(deltas) == len(idx_old):
+1341        raise Exception('Length of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old)))
+1342    if type(idx_old) is range and type(idx_new) is range:
+1343        if idx_old == idx_new:
+1344            return deltas
+1345    if _check_lists_equal([idx_old, idx_new]):
+1346        return deltas
+1347    indices = np.intersect1d(idx_old, idx_new, assume_unique=True, return_indices=True)[1]
+1348    if len(indices) < len(idx_new):
+1349        raise Exception('Error in _reduce_deltas: Config of idx_new not in idx_old')
+1350    return np.array(deltas)[indices]
+1351
+1352
+1353def reweight(weight, obs, **kwargs):
+1354    """Reweight a list of observables.
+1355
+1356    Parameters
+1357    ----------
+1358    weight : Obs
+1359        Reweighting factor. An Observable that has to be defined on a superset of the
+1360        configurations in obs[i].idl for all i.
+1361    obs : list
+1362        list of Obs, e.g. [obs1, obs2, obs3].
+1363    all_configs : bool
+1364        if True, the reweighted observables are normalized by the average of
+1365        the reweighting factor on all configurations in weight.idl and not
+1366        on the configurations in obs[i].idl. Default False.
+1367    """
+1368    result = []
+1369    for i in range(len(obs)):
+1370        if len(obs[i].cov_names):
+1371            raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
+1372        if not set(obs[i].names).issubset(weight.names):
+1373            raise Exception('Error: Ensembles do not fit')
+1374        for name in obs[i].names:
+1375            if not set(obs[i].idl[name]).issubset(weight.idl[name]):
+1376                raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
+1377        new_samples = []
+1378        w_deltas = {}
+1379        for name in sorted(obs[i].names):
+1380            w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
+1381            new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
+1382        tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
 1383
-1384    return result
-1385
-1386
-1387def correlate(obs_a, obs_b):
-1388    """Correlate two observables.
-1389
-1390    Parameters
-1391    ----------
-1392    obs_a : Obs
-1393        First observable
-1394    obs_b : Obs
-1395        Second observable
-1396
-1397    Notes
-1398    -----
-1399    Keep in mind to only correlate primary observables which have not been reweighted
-1400    yet. The reweighting has to be applied after correlating the observables.
-1401    Currently only works if ensembles are identical (this is not strictly necessary).
-1402    """
-1403
-1404    if sorted(obs_a.names) != sorted(obs_b.names):
-1405        raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
-1406    if len(obs_a.cov_names) or len(obs_b.cov_names):
-1407        raise Exception('Error: Not possible to correlate Obs that contain covobs!')
-1408    for name in obs_a.names:
-1409        if obs_a.shape[name] != obs_b.shape[name]:
-1410            raise Exception('Shapes of ensemble', name, 'do not fit')
-1411        if obs_a.idl[name] != obs_b.idl[name]:
-1412            raise Exception('idl of ensemble', name, 'do not fit')
-1413
-1414    if obs_a.reweighted is True:
-1415        warnings.warn("The first observable is already reweighted.", RuntimeWarning)
-1416    if obs_b.reweighted is True:
-1417        warnings.warn("The second observable is already reweighted.", RuntimeWarning)
-1418
-1419    new_samples = []
-1420    new_idl = []
-1421    for name in sorted(obs_a.names):
-1422        new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
-1423        new_idl.append(obs_a.idl[name])
-1424
-1425    o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
-1426    o.reweighted = obs_a.reweighted or obs_b.reweighted
-1427    return o
-1428
-1429
-1430def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
-1431    r'''Calculates the error covariance matrix of a set of observables.
+1384        if kwargs.get('all_configs'):
+1385            new_weight = weight
+1386        else:
+1387            new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
+1388
+1389        result.append(tmp_obs / new_weight)
+1390        result[-1].reweighted = True
+1391
+1392    return result
+1393
+1394
+1395def correlate(obs_a, obs_b):
+1396    """Correlate two observables.
+1397
+1398    Parameters
+1399    ----------
+1400    obs_a : Obs
+1401        First observable
+1402    obs_b : Obs
+1403        Second observable
+1404
+1405    Notes
+1406    -----
+1407    Keep in mind to only correlate primary observables which have not been reweighted
+1408    yet. The reweighting has to be applied after correlating the observables.
+1409    Currently only works if ensembles are identical (this is not strictly necessary).
+1410    """
+1411
+1412    if sorted(obs_a.names) != sorted(obs_b.names):
+1413        raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
+1414    if len(obs_a.cov_names) or len(obs_b.cov_names):
+1415        raise Exception('Error: Not possible to correlate Obs that contain covobs!')
+1416    for name in obs_a.names:
+1417        if obs_a.shape[name] != obs_b.shape[name]:
+1418            raise Exception('Shapes of ensemble', name, 'do not fit')
+1419        if obs_a.idl[name] != obs_b.idl[name]:
+1420            raise Exception('idl of ensemble', name, 'do not fit')
+1421
+1422    if obs_a.reweighted is True:
+1423        warnings.warn("The first observable is already reweighted.", RuntimeWarning)
+1424    if obs_b.reweighted is True:
+1425        warnings.warn("The second observable is already reweighted.", RuntimeWarning)
+1426
+1427    new_samples = []
+1428    new_idl = []
+1429    for name in sorted(obs_a.names):
+1430        new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
+1431        new_idl.append(obs_a.idl[name])
 1432
-1433    WARNING: This function should be used with care, especially for observables with support on multiple
-1434             ensembles with differing autocorrelations. See the notes below for details.
-1435
-1436    The gamma method has to be applied first to all observables.
+1433    o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
+1434    o.reweighted = obs_a.reweighted or obs_b.reweighted
+1435    return o
+1436
 1437
-1438    Parameters
-1439    ----------
-1440    obs : list or numpy.ndarray
-1441        List or one dimensional array of Obs
-1442    visualize : bool
-1443        If True plots the corresponding normalized correlation matrix (default False).
-1444    correlation : bool
-1445        If True the correlation matrix instead of the error covariance matrix is returned (default False).
-1446    smooth : None or int
-1447        If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
-1448        smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
-1449        largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
-1450        small ones.
-1451
-1452    Notes
-1453    -----
-1454    The error covariance is defined such that it agrees with the squared standard error for two identical observables
-1455    $$\operatorname{cov}(a,a)=\sum_{s=1}^N\delta_a^s\delta_a^s/N^2=\Gamma_{aa}(0)/N=\operatorname{var}(a)/N=\sigma_a^2$$
-1456    in the absence of autocorrelation.
-1457    The error covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. The covariance at windowsize 0 is guaranteed to be positive semi-definite
-1458    $$\sum_{i,j}v_i\Gamma_{ij}(0)v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i,j}v_i\delta_i^s\delta_j^s v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i}|v_i\delta_i^s|^2\geq 0\,,$$ for every $v\in\mathbb{R}^M$, while such an identity does not hold for larger windows/lags.
-1459    For observables defined on a single ensemble our approximation is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
-1460    $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
-1461    This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
-1462    '''
-1463
-1464    length = len(obs)
-1465
-1466    max_samples = np.max([o.N for o in obs])
-1467    if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
-1468        warnings.warn(f"The dimension of the covariance matrix ({length}) is larger or equal to the number of samples ({max_samples}). This will result in a rank deficient matrix.", RuntimeWarning)
-1469
-1470    cov = np.zeros((length, length))
-1471    for i in range(length):
-1472        for j in range(i, length):
-1473            cov[i, j] = _covariance_element(obs[i], obs[j])
-1474    cov = cov + cov.T - np.diag(np.diag(cov))
-1475
-1476    corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+1438def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
+1439    r'''Calculates the error covariance matrix of a set of observables.
+1440
+1441    WARNING: This function should be used with care, especially for observables with support on multiple
+1442             ensembles with differing autocorrelations. See the notes below for details.
+1443
+1444    The gamma method has to be applied first to all observables.
+1445
+1446    Parameters
+1447    ----------
+1448    obs : list or numpy.ndarray
+1449        List or one dimensional array of Obs
+1450    visualize : bool
+1451        If True plots the corresponding normalized correlation matrix (default False).
+1452    correlation : bool
+1453        If True the correlation matrix instead of the error covariance matrix is returned (default False).
+1454    smooth : None or int
+1455        If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
+1456        smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
+1457        largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
+1458        small ones.
+1459
+1460    Notes
+1461    -----
+1462    The error covariance is defined such that it agrees with the squared standard error for two identical observables
+1463    $$\operatorname{cov}(a,a)=\sum_{s=1}^N\delta_a^s\delta_a^s/N^2=\Gamma_{aa}(0)/N=\operatorname{var}(a)/N=\sigma_a^2$$
+1464    in the absence of autocorrelation.
+1465    The error covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. The covariance at windowsize 0 is guaranteed to be positive semi-definite
+1466    $$\sum_{i,j}v_i\Gamma_{ij}(0)v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i,j}v_i\delta_i^s\delta_j^s v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i}|v_i\delta_i^s|^2\geq 0\,,$$ for every $v\in\mathbb{R}^M$, while such an identity does not hold for larger windows/lags.
+1467    For observables defined on a single ensemble our approximation is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
+1468    $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
+1469    This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
+1470    '''
+1471
+1472    length = len(obs)
+1473
+1474    max_samples = np.max([o.N for o in obs])
+1475    if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
+1476        warnings.warn(f"The dimension of the covariance matrix ({length}) is larger or equal to the number of samples ({max_samples}). This will result in a rank deficient matrix.", RuntimeWarning)
 1477
-1478    if isinstance(smooth, int):
-1479        corr = _smooth_eigenvalues(corr, smooth)
-1480
-1481    if visualize:
-1482        plt.matshow(corr, vmin=-1, vmax=1)
-1483        plt.set_cmap('RdBu')
-1484        plt.colorbar()
-1485        plt.draw()
-1486
-1487    if correlation is True:
-1488        return corr
-1489
-1490    errors = [o.dvalue for o in obs]
-1491    cov = np.diag(errors) @ corr @ np.diag(errors)
-1492
-1493    eigenvalues = np.linalg.eigh(cov)[0]
-1494    if not np.all(eigenvalues >= 0):
-1495        warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
-1496
-1497    return cov
-1498
-1499
-1500def _smooth_eigenvalues(corr, E):
-1501    """Eigenvalue smoothing as described in hep-lat/9412087
-1502
-1503    corr : np.ndarray
-1504        correlation matrix
-1505    E : integer
-1506        Number of eigenvalues to be left substantially unchanged
-1507    """
-1508    if not (2 < E < corr.shape[0] - 1):
-1509        raise Exception(f"'E' has to be between 2 and the dimension of the correlation matrix minus 1 ({corr.shape[0] - 1}).")
-1510    vals, vec = np.linalg.eigh(corr)
-1511    lambda_min = np.mean(vals[:-E])
-1512    vals[vals < lambda_min] = lambda_min
-1513    vals /= np.mean(vals)
-1514    return vec @ np.diag(vals) @ vec.T
-1515
-1516
-1517def _covariance_element(obs1, obs2):
-1518    """Estimates the covariance of two Obs objects, neglecting autocorrelations."""
-1519
-1520    def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx):
-1521        deltas1 = _reduce_deltas(deltas1, idx1, new_idx)
-1522        deltas2 = _reduce_deltas(deltas2, idx2, new_idx)
-1523        return np.sum(deltas1 * deltas2)
+1478    cov = np.zeros((length, length))
+1479    for i in range(length):
+1480        for j in range(i, length):
+1481            cov[i, j] = _covariance_element(obs[i], obs[j])
+1482    cov = cov + cov.T - np.diag(np.diag(cov))
+1483
+1484    corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+1485
+1486    if isinstance(smooth, int):
+1487        corr = _smooth_eigenvalues(corr, smooth)
+1488
+1489    if visualize:
+1490        plt.matshow(corr, vmin=-1, vmax=1)
+1491        plt.set_cmap('RdBu')
+1492        plt.colorbar()
+1493        plt.draw()
+1494
+1495    if correlation is True:
+1496        return corr
+1497
+1498    errors = [o.dvalue for o in obs]
+1499    cov = np.diag(errors) @ corr @ np.diag(errors)
+1500
+1501    eigenvalues = np.linalg.eigh(cov)[0]
+1502    if not np.all(eigenvalues >= 0):
+1503        warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
+1504
+1505    return cov
+1506
+1507
+1508def _smooth_eigenvalues(corr, E):
+1509    """Eigenvalue smoothing as described in hep-lat/9412087
+1510
+1511    corr : np.ndarray
+1512        correlation matrix
+1513    E : integer
+1514        Number of eigenvalues to be left substantially unchanged
+1515    """
+1516    if not (2 < E < corr.shape[0] - 1):
+1517        raise Exception(f"'E' has to be between 2 and the dimension of the correlation matrix minus 1 ({corr.shape[0] - 1}).")
+1518    vals, vec = np.linalg.eigh(corr)
+1519    lambda_min = np.mean(vals[:-E])
+1520    vals[vals < lambda_min] = lambda_min
+1521    vals /= np.mean(vals)
+1522    return vec @ np.diag(vals) @ vec.T
+1523
 1524
-1525    if set(obs1.names).isdisjoint(set(obs2.names)):
-1526        return 0.0
+1525def _covariance_element(obs1, obs2):
+1526    """Estimates the covariance of two Obs objects, neglecting autocorrelations."""
 1527
-1528    if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'):
-1529        raise Exception('The gamma method has to be applied to both Obs first.')
-1530
-1531    dvalue = 0.0
+1528    def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx):
+1529        deltas1 = _reduce_deltas(deltas1, idx1, new_idx)
+1530        deltas2 = _reduce_deltas(deltas2, idx2, new_idx)
+1531        return np.sum(deltas1 * deltas2)
 1532
-1533    for e_name in obs1.mc_names:
-1534
-1535        if e_name not in obs2.mc_names:
-1536            continue
-1537
-1538        idl_d = {}
-1539        for r_name in obs1.e_content[e_name]:
-1540            if r_name not in obs2.e_content[e_name]:
-1541                continue
-1542            idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]])
-1543
-1544        gamma = 0.0
+1533    if set(obs1.names).isdisjoint(set(obs2.names)):
+1534        return 0.0
+1535
+1536    if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'):
+1537        raise Exception('The gamma method has to be applied to both Obs first.')
+1538
+1539    dvalue = 0.0
+1540
+1541    for e_name in obs1.mc_names:
+1542
+1543        if e_name not in obs2.mc_names:
+1544            continue
 1545
-1546        for r_name in obs1.e_content[e_name]:
-1547            if r_name not in obs2.e_content[e_name]:
-1548                continue
-1549            if len(idl_d[r_name]) == 0:
-1550                continue
-1551            gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name])
-1552
-1553        if gamma == 0.0:
-1554            continue
-1555
-1556        gamma_div = 0.0
-1557        for r_name in obs1.e_content[e_name]:
-1558            if r_name not in obs2.e_content[e_name]:
-1559                continue
-1560            if len(idl_d[r_name]) == 0:
-1561                continue
-1562            gamma_div += np.sqrt(calc_gamma(obs1.deltas[r_name], obs1.deltas[r_name], obs1.idl[r_name], obs1.idl[r_name], idl_d[r_name]) * calc_gamma(obs2.deltas[r_name], obs2.deltas[r_name], obs2.idl[r_name], obs2.idl[r_name], idl_d[r_name]))
-1563        gamma /= gamma_div
-1564
-1565        dvalue += gamma
-1566
-1567    for e_name in obs1.cov_names:
-1568
-1569        if e_name not in obs2.cov_names:
-1570            continue
-1571
-1572        dvalue += np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad)).item()
-1573
-1574    return dvalue
-1575
+1546        idl_d = {}
+1547        for r_name in obs1.e_content[e_name]:
+1548            if r_name not in obs2.e_content[e_name]:
+1549                continue
+1550            idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]])
+1551
+1552        gamma = 0.0
+1553
+1554        for r_name in obs1.e_content[e_name]:
+1555            if r_name not in obs2.e_content[e_name]:
+1556                continue
+1557            if len(idl_d[r_name]) == 0:
+1558                continue
+1559            gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name])
+1560
+1561        if gamma == 0.0:
+1562            continue
+1563
+1564        gamma_div = 0.0
+1565        for r_name in obs1.e_content[e_name]:
+1566            if r_name not in obs2.e_content[e_name]:
+1567                continue
+1568            if len(idl_d[r_name]) == 0:
+1569                continue
+1570            gamma_div += np.sqrt(calc_gamma(obs1.deltas[r_name], obs1.deltas[r_name], obs1.idl[r_name], obs1.idl[r_name], idl_d[r_name]) * calc_gamma(obs2.deltas[r_name], obs2.deltas[r_name], obs2.idl[r_name], obs2.idl[r_name], idl_d[r_name]))
+1571        gamma /= gamma_div
+1572
+1573        dvalue += gamma
+1574
+1575    for e_name in obs1.cov_names:
 1576
-1577def import_jackknife(jacks, name, idl=None):
-1578    """Imports jackknife samples and returns an Obs
+1577        if e_name not in obs2.cov_names:
+1578            continue
 1579
-1580    Parameters
-1581    ----------
-1582    jacks : numpy.ndarray
-1583        numpy array containing the mean value as zeroth entry and
-1584        the N jackknife samples as first to Nth entry.
-1585    name : str
-1586        name of the ensemble the samples are defined on.
-1587    """
-1588    length = len(jacks) - 1
-1589    prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
-1590    samples = jacks[1:] @ prj
-1591    mean = np.mean(samples)
-1592    new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
-1593    new_obs._value = jacks[0]
-1594    return new_obs
-1595
-1596
-1597def import_bootstrap(boots, name, random_numbers):
-1598    """Imports bootstrap samples and returns an Obs
-1599
-1600    Parameters
-1601    ----------
-1602    boots : numpy.ndarray
-1603        numpy array containing the mean value as zeroth entry and
-1604        the N bootstrap samples as first to Nth entry.
-1605    name : str
-1606        name of the ensemble the samples are defined on.
-1607    random_numbers : np.ndarray
-1608        Array of shape (samples, length) containing the random numbers to generate the bootstrap samples,
-1609        where samples is the number of bootstrap samples and length is the length of the original Monte Carlo
-1610        chain to be reconstructed.
-1611    """
-1612    samples, length = random_numbers.shape
-1613    if samples != len(boots) - 1:
-1614        raise ValueError("Random numbers do not have the correct shape.")
-1615
-1616    if samples < length:
-1617        raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.")
-1618
-1619    proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
-1620
-1621    samples = scipy.linalg.lstsq(proj, boots[1:])[0]
-1622    ret = Obs([samples], [name])
-1623    ret._value = boots[0]
-1624    return ret
-1625
+1580        dvalue += np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad)).item()
+1581
+1582    return dvalue
+1583
+1584
+1585def import_jackknife(jacks, name, idl=None):
+1586    """Imports jackknife samples and returns an Obs
+1587
+1588    Parameters
+1589    ----------
+1590    jacks : numpy.ndarray
+1591        numpy array containing the mean value as zeroth entry and
+1592        the N jackknife samples as first to Nth entry.
+1593    name : str
+1594        name of the ensemble the samples are defined on.
+1595    """
+1596    length = len(jacks) - 1
+1597    prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
+1598    samples = jacks[1:] @ prj
+1599    mean = np.mean(samples)
+1600    new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
+1601    new_obs._value = jacks[0]
+1602    return new_obs
+1603
+1604
+1605def import_bootstrap(boots, name, random_numbers):
+1606    """Imports bootstrap samples and returns an Obs
+1607
+1608    Parameters
+1609    ----------
+1610    boots : numpy.ndarray
+1611        numpy array containing the mean value as zeroth entry and
+1612        the N bootstrap samples as first to Nth entry.
+1613    name : str
+1614        name of the ensemble the samples are defined on.
+1615    random_numbers : np.ndarray
+1616        Array of shape (samples, length) containing the random numbers to generate the bootstrap samples,
+1617        where samples is the number of bootstrap samples and length is the length of the original Monte Carlo
+1618        chain to be reconstructed.
+1619    """
+1620    samples, length = random_numbers.shape
+1621    if samples != len(boots) - 1:
+1622        raise ValueError("Random numbers do not have the correct shape.")
+1623
+1624    if samples < length:
+1625        raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.")
 1626
-1627def merge_obs(list_of_obs):
-1628    """Combine all observables in list_of_obs into one new observable
-1629
-1630    Parameters
-1631    ----------
-1632    list_of_obs : list
-1633        list of the Obs object to be combined
+1627    proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
+1628
+1629    samples = scipy.linalg.lstsq(proj, boots[1:])[0]
+1630    ret = Obs([samples], [name])
+1631    ret._value = boots[0]
+1632    return ret
+1633
 1634
-1635    Notes
-1636    -----
-1637    It is not possible to combine obs which are based on the same replicum
-1638    """
-1639    replist = [item for obs in list_of_obs for item in obs.names]
-1640    if (len(replist) == len(set(replist))) is False:
-1641        raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
-1642    if any([len(o.cov_names) for o in list_of_obs]):
-1643        raise Exception('Not possible to merge data that contains covobs!')
-1644    new_dict = {}
-1645    idl_dict = {}
-1646    for o in list_of_obs:
-1647        new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
-1648                        for key in set(o.deltas) | set(o.r_values)})
-1649        idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
-1650
-1651    names = sorted(new_dict.keys())
-1652    o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
-1653    o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
-1654    return o
-1655
-1656
-1657def cov_Obs(means, cov, name, grad=None):
-1658    """Create an Obs based on mean(s) and a covariance matrix
-1659
-1660    Parameters
-1661    ----------
-1662    mean : list of floats or float
-1663        N mean value(s) of the new Obs
-1664    cov : list or array
-1665        2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
-1666    name : str
-1667        identifier for the covariance matrix
-1668    grad : list or array
-1669        Gradient of the Covobs wrt. the means belonging to cov.
-1670    """
-1671
-1672    def covobs_to_obs(co):
-1673        """Make an Obs out of a Covobs
-1674
-1675        Parameters
-1676        ----------
-1677        co : Covobs
-1678            Covobs to be embedded into the Obs
-1679        """
-1680        o = Obs([], [], means=[])
-1681        o._value = co.value
-1682        o.names.append(co.name)
-1683        o._covobs[co.name] = co
-1684        o._dvalue = np.sqrt(co.errsq())
-1685        return o
-1686
-1687    ol = []
-1688    if isinstance(means, (float, int)):
-1689        means = [means]
-1690
-1691    for i in range(len(means)):
-1692        ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
-1693    if ol[0].covobs[name].N != len(means):
-1694        raise Exception('You have to provide %d mean values!' % (ol[0].N))
-1695    if len(ol) == 1:
-1696        return ol[0]
-1697    return ol
+1635def merge_obs(list_of_obs):
+1636    """Combine all observables in list_of_obs into one new observable
+1637
+1638    Parameters
+1639    ----------
+1640    list_of_obs : list
+1641        list of the Obs object to be combined
+1642
+1643    Notes
+1644    -----
+1645    It is not possible to combine obs which are based on the same replicum
+1646    """
+1647    replist = [item for obs in list_of_obs for item in obs.names]
+1648    if (len(replist) == len(set(replist))) is False:
+1649        raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
+1650    if any([len(o.cov_names) for o in list_of_obs]):
+1651        raise Exception('Not possible to merge data that contains covobs!')
+1652    new_dict = {}
+1653    idl_dict = {}
+1654    for o in list_of_obs:
+1655        new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
+1656                        for key in set(o.deltas) | set(o.r_values)})
+1657        idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
+1658
+1659    names = sorted(new_dict.keys())
+1660    o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
+1661    o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
+1662    return o
+1663
+1664
+1665def cov_Obs(means, cov, name, grad=None):
+1666    """Create an Obs based on mean(s) and a covariance matrix
+1667
+1668    Parameters
+1669    ----------
+1670    mean : list of floats or float
+1671        N mean value(s) of the new Obs
+1672    cov : list or array
+1673        2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
+1674    name : str
+1675        identifier for the covariance matrix
+1676    grad : list or array
+1677        Gradient of the Covobs wrt. the means belonging to cov.
+1678    """
+1679
+1680    def covobs_to_obs(co):
+1681        """Make an Obs out of a Covobs
+1682
+1683        Parameters
+1684        ----------
+1685        co : Covobs
+1686            Covobs to be embedded into the Obs
+1687        """
+1688        o = Obs([], [], means=[])
+1689        o._value = co.value
+1690        o.names.append(co.name)
+1691        o._covobs[co.name] = co
+1692        o._dvalue = np.sqrt(co.errsq())
+1693        return o
+1694
+1695    ol = []
+1696    if isinstance(means, (float, int)):
+1697        means = [means]
 1698
-1699
-1700def _determine_gap(o, e_content, e_name):
-1701    gaps = []
-1702    for r_name in e_content[e_name]:
-1703        if isinstance(o.idl[r_name], range):
-1704            gaps.append(o.idl[r_name].step)
-1705        else:
-1706            gaps.append(np.min(np.diff(o.idl[r_name])))
+1699    for i in range(len(means)):
+1700        ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
+1701    if ol[0].covobs[name].N != len(means):
+1702        raise Exception('You have to provide %d mean values!' % (ol[0].N))
+1703    if len(ol) == 1:
+1704        return ol[0]
+1705    return ol
+1706
 1707
-1708    gap = min(gaps)
-1709    if not np.all([gi % gap == 0 for gi in gaps]):
-1710        raise Exception(f"Replica for ensemble {e_name} do not have a common spacing.", gaps)
-1711
-1712    return gap
-1713
-1714
-1715def _check_lists_equal(idl):
-1716    '''
-1717    Use groupby to efficiently check whether all elements of idl are identical.
-1718    Returns True if all elements are equal, otherwise False.
+1708def _determine_gap(o, e_content, e_name):
+1709    gaps = []
+1710    for r_name in e_content[e_name]:
+1711        if isinstance(o.idl[r_name], range):
+1712            gaps.append(o.idl[r_name].step)
+1713        else:
+1714            gaps.append(np.min(np.diff(o.idl[r_name])))
+1715
+1716    gap = min(gaps)
+1717    if not np.all([gi % gap == 0 for gi in gaps]):
+1718        raise Exception(f"Replica for ensemble {e_name} do not have a common spacing.", gaps)
 1719
-1720    Parameters
-1721    ----------
-1722    idl : list of lists, ranges or np.ndarrays
-1723    '''
-1724    g = groupby([np.nditer(el) if isinstance(el, np.ndarray) else el for el in idl])
-1725    if next(g, True) and not next(g, False):
-1726        return True
-1727    return False
+1720    return gap
+1721
+1722
+1723def _check_lists_equal(idl):
+1724    '''
+1725    Use groupby to efficiently check whether all elements of idl are identical.
+1726    Returns True if all elements are equal, otherwise False.
+1727
+1728    Parameters
+1729    ----------
+1730    idl : list of lists, ranges or np.ndarrays
+1731    '''
+1732    g = groupby([np.nditer(el) if isinstance(el, np.ndarray) else el for el in idl])
+1733    if next(g, True) and not next(g, False):
+1734        return True
+1735    return False
@@ -4976,6 +4984,14 @@ should agree with samples from a full bootstrap analysis up to O(1/N). 1023 1024 def __repr__(self): 1025 return 'CObs[' + str(self) + ']' +1026 +1027 def __format__(self, format_type): +1028 if format_type == "": +1029 significance = 2 +1030 format_type = "2" +1031 else: +1032 significance = int(float(format_type.replace("+", "").replace("-", ""))) +1033 return f"({self.real:{format_type}}{self.imag:+{significance}}j)" @@ -5114,174 +5130,174 @@ should agree with samples from a full bootstrap analysis up to O(1/N). - 
1147def derived_observable(func, data, array_mode=False, **kwargs):
-1148    """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
-1149
-1150    Parameters
-1151    ----------
-1152    func : object
-1153        arbitrary function of the form func(data, **kwargs). For the
-1154        automatic differentiation to work, all numpy functions have to have
-1155        the autograd wrapper (use 'import autograd.numpy as anp').
-1156    data : list
-1157        list of Obs, e.g. [obs1, obs2, obs3].
-1158    num_grad : bool
-1159        if True, numerical derivatives are used instead of autograd
-1160        (default False). To control the numerical differentiation the
-1161        kwargs of numdifftools.step_generators.MaxStepGenerator
-1162        can be used.
-1163    man_grad : list
-1164        manually supply a list or an array which contains the jacobian
-1165        of func. Use cautiously, supplying the wrong derivative will
-1166        not be intercepted.
-1167
-1168    Notes
-1169    -----
-1170    For simple mathematical operations it can be practical to use anonymous
-1171    functions. For the ratio of two observables one can e.g. use
-1172
-1173    new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
-1174    """
+            
1155def derived_observable(func, data, array_mode=False, **kwargs):
+1156    """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
+1157
+1158    Parameters
+1159    ----------
+1160    func : object
+1161        arbitrary function of the form func(data, **kwargs). For the
+1162        automatic differentiation to work, all numpy functions have to have
+1163        the autograd wrapper (use 'import autograd.numpy as anp').
+1164    data : list
+1165        list of Obs, e.g. [obs1, obs2, obs3].
+1166    num_grad : bool
+1167        if True, numerical derivatives are used instead of autograd
+1168        (default False). To control the numerical differentiation the
+1169        kwargs of numdifftools.step_generators.MaxStepGenerator
+1170        can be used.
+1171    man_grad : list
+1172        manually supply a list or an array which contains the jacobian
+1173        of func. Use cautiously, supplying the wrong derivative will
+1174        not be intercepted.
 1175
-1176    data = np.asarray(data)
-1177    raveled_data = data.ravel()
-1178
-1179    # Workaround for matrix operations containing non Obs data
-1180    if not all(isinstance(x, Obs) for x in raveled_data):
-1181        for i in range(len(raveled_data)):
-1182            if isinstance(raveled_data[i], (int, float)):
-1183                raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
-1184
-1185    allcov = {}
-1186    for o in raveled_data:
-1187        for name in o.cov_names:
-1188            if name in allcov:
-1189                if not np.allclose(allcov[name], o.covobs[name].cov):
-1190                    raise Exception('Inconsistent covariance matrices for %s!' % (name))
-1191            else:
-1192                allcov[name] = o.covobs[name].cov
-1193
-1194    n_obs = len(raveled_data)
-1195    new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
-1196    new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
-1197    new_sample_names = sorted(set(new_names) - set(new_cov_names))
-1198
-1199    reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
-1200
-1201    if data.ndim == 1:
-1202        values = np.array([o.value for o in data])
-1203    else:
-1204        values = np.vectorize(lambda x: x.value)(data)
-1205
-1206    new_values = func(values, **kwargs)
-1207
-1208    multi = int(isinstance(new_values, np.ndarray))
-1209
-1210    new_r_values = {}
-1211    new_idl_d = {}
-1212    for name in new_sample_names:
-1213        idl = []
-1214        tmp_values = np.zeros(n_obs)
-1215        for i, item in enumerate(raveled_data):
-1216            tmp_values[i] = item.r_values.get(name, item.value)
-1217            tmp_idl = item.idl.get(name)
-1218            if tmp_idl is not None:
-1219                idl.append(tmp_idl)
-1220        if multi > 0:
-1221            tmp_values = np.array(tmp_values).reshape(data.shape)
-1222        new_r_values[name] = func(tmp_values, **kwargs)
-1223        new_idl_d[name] = _merge_idx(idl)
-1224
-1225    if 'man_grad' in kwargs:
-1226        deriv = np.asarray(kwargs.get('man_grad'))
-1227        if new_values.shape + data.shape != deriv.shape:
-1228            raise Exception('Manual derivative does not have correct shape.')
-1229    elif kwargs.get('num_grad') is True:
-1230        if multi > 0:
-1231            raise Exception('Multi mode currently not supported for numerical derivative')
-1232        options = {
-1233            'base_step': 0.1,
-1234            'step_ratio': 2.5}
-1235        for key in options.keys():
-1236            kwarg = kwargs.get(key)
-1237            if kwarg is not None:
-1238                options[key] = kwarg
-1239        tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
-1240        if tmp_df.size == 1:
-1241            deriv = np.array([tmp_df.real])
-1242        else:
-1243            deriv = tmp_df.real
-1244    else:
-1245        deriv = jacobian(func)(values, **kwargs)
-1246
-1247    final_result = np.zeros(new_values.shape, dtype=object)
-1248
-1249    if array_mode is True:
-1250
-1251        class _Zero_grad():
-1252            def __init__(self, N):
-1253                self.grad = np.zeros((N, 1))
+1176    Notes
+1177    -----
+1178    For simple mathematical operations it can be practical to use anonymous
+1179    functions. For the ratio of two observables one can e.g. use
+1180
+1181    new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
+1182    """
+1183
+1184    data = np.asarray(data)
+1185    raveled_data = data.ravel()
+1186
+1187    # Workaround for matrix operations containing non Obs data
+1188    if not all(isinstance(x, Obs) for x in raveled_data):
+1189        for i in range(len(raveled_data)):
+1190            if isinstance(raveled_data[i], (int, float)):
+1191                raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
+1192
+1193    allcov = {}
+1194    for o in raveled_data:
+1195        for name in o.cov_names:
+1196            if name in allcov:
+1197                if not np.allclose(allcov[name], o.covobs[name].cov):
+1198                    raise Exception('Inconsistent covariance matrices for %s!' % (name))
+1199            else:
+1200                allcov[name] = o.covobs[name].cov
+1201
+1202    n_obs = len(raveled_data)
+1203    new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
+1204    new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
+1205    new_sample_names = sorted(set(new_names) - set(new_cov_names))
+1206
+1207    reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
+1208
+1209    if data.ndim == 1:
+1210        values = np.array([o.value for o in data])
+1211    else:
+1212        values = np.vectorize(lambda x: x.value)(data)
+1213
+1214    new_values = func(values, **kwargs)
+1215
+1216    multi = int(isinstance(new_values, np.ndarray))
+1217
+1218    new_r_values = {}
+1219    new_idl_d = {}
+1220    for name in new_sample_names:
+1221        idl = []
+1222        tmp_values = np.zeros(n_obs)
+1223        for i, item in enumerate(raveled_data):
+1224            tmp_values[i] = item.r_values.get(name, item.value)
+1225            tmp_idl = item.idl.get(name)
+1226            if tmp_idl is not None:
+1227                idl.append(tmp_idl)
+1228        if multi > 0:
+1229            tmp_values = np.array(tmp_values).reshape(data.shape)
+1230        new_r_values[name] = func(tmp_values, **kwargs)
+1231        new_idl_d[name] = _merge_idx(idl)
+1232
+1233    if 'man_grad' in kwargs:
+1234        deriv = np.asarray(kwargs.get('man_grad'))
+1235        if new_values.shape + data.shape != deriv.shape:
+1236            raise Exception('Manual derivative does not have correct shape.')
+1237    elif kwargs.get('num_grad') is True:
+1238        if multi > 0:
+1239            raise Exception('Multi mode currently not supported for numerical derivative')
+1240        options = {
+1241            'base_step': 0.1,
+1242            'step_ratio': 2.5}
+1243        for key in options.keys():
+1244            kwarg = kwargs.get(key)
+1245            if kwarg is not None:
+1246                options[key] = kwarg
+1247        tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
+1248        if tmp_df.size == 1:
+1249            deriv = np.array([tmp_df.real])
+1250        else:
+1251            deriv = tmp_df.real
+1252    else:
+1253        deriv = jacobian(func)(values, **kwargs)
 1254
-1255        new_covobs_lengths = dict(set([y for x in [[(n, o.covobs[n].N) for n in o.cov_names] for o in raveled_data] for y in x]))
-1256        d_extracted = {}
-1257        g_extracted = {}
-1258        for name in new_sample_names:
-1259            d_extracted[name] = []
-1260            ens_length = len(new_idl_d[name])
-1261            for i_dat, dat in enumerate(data):
-1262                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, )))
-1263        for name in new_cov_names:
-1264            g_extracted[name] = []
-1265            zero_grad = _Zero_grad(new_covobs_lengths[name])
-1266            for i_dat, dat in enumerate(data):
-1267                g_extracted[name].append(np.array([o.covobs.get(name, zero_grad).grad for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (new_covobs_lengths[name], 1)))
-1268
-1269    for i_val, new_val in np.ndenumerate(new_values):
-1270        new_deltas = {}
-1271        new_grad = {}
-1272        if array_mode is True:
-1273            for name in new_sample_names:
-1274                ens_length = d_extracted[name][0].shape[-1]
-1275                new_deltas[name] = np.zeros(ens_length)
-1276                for i_dat, dat in enumerate(d_extracted[name]):
-1277                    new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1278            for name in new_cov_names:
-1279                new_grad[name] = 0
-1280                for i_dat, dat in enumerate(g_extracted[name]):
-1281                    new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1282        else:
-1283            for j_obs, obs in np.ndenumerate(data):
-1284                for name in obs.names:
-1285                    if name in obs.cov_names:
-1286                        new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
-1287                    else:
-1288                        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])
-1289
-1290        new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
-1291
-1292        if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
-1293            raise Exception('The same name has been used for deltas and covobs!')
-1294        new_samples = []
-1295        new_means = []
-1296        new_idl = []
-1297        new_names_obs = []
-1298        for name in new_names:
-1299            if name not in new_covobs:
-1300                new_samples.append(new_deltas[name])
-1301                new_idl.append(new_idl_d[name])
-1302                new_means.append(new_r_values[name][i_val])
-1303                new_names_obs.append(name)
-1304        final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
-1305        for name in new_covobs:
-1306            final_result[i_val].names.append(name)
-1307        final_result[i_val]._covobs = new_covobs
-1308        final_result[i_val]._value = new_val
-1309        final_result[i_val].reweighted = reweighted
-1310
-1311    if multi == 0:
-1312        final_result = final_result.item()
-1313
-1314    return final_result
+1255    final_result = np.zeros(new_values.shape, dtype=object)
+1256
+1257    if array_mode is True:
+1258
+1259        class _Zero_grad():
+1260            def __init__(self, N):
+1261                self.grad = np.zeros((N, 1))
+1262
+1263        new_covobs_lengths = dict(set([y for x in [[(n, o.covobs[n].N) for n in o.cov_names] for o in raveled_data] for y in x]))
+1264        d_extracted = {}
+1265        g_extracted = {}
+1266        for name in new_sample_names:
+1267            d_extracted[name] = []
+1268            ens_length = len(new_idl_d[name])
+1269            for i_dat, dat in enumerate(data):
+1270                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, )))
+1271        for name in new_cov_names:
+1272            g_extracted[name] = []
+1273            zero_grad = _Zero_grad(new_covobs_lengths[name])
+1274            for i_dat, dat in enumerate(data):
+1275                g_extracted[name].append(np.array([o.covobs.get(name, zero_grad).grad for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (new_covobs_lengths[name], 1)))
+1276
+1277    for i_val, new_val in np.ndenumerate(new_values):
+1278        new_deltas = {}
+1279        new_grad = {}
+1280        if array_mode is True:
+1281            for name in new_sample_names:
+1282                ens_length = d_extracted[name][0].shape[-1]
+1283                new_deltas[name] = np.zeros(ens_length)
+1284                for i_dat, dat in enumerate(d_extracted[name]):
+1285                    new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1286            for name in new_cov_names:
+1287                new_grad[name] = 0
+1288                for i_dat, dat in enumerate(g_extracted[name]):
+1289                    new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1290        else:
+1291            for j_obs, obs in np.ndenumerate(data):
+1292                for name in obs.names:
+1293                    if name in obs.cov_names:
+1294                        new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
+1295                    else:
+1296                        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])
+1297
+1298        new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
+1299
+1300        if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
+1301            raise Exception('The same name has been used for deltas and covobs!')
+1302        new_samples = []
+1303        new_means = []
+1304        new_idl = []
+1305        new_names_obs = []
+1306        for name in new_names:
+1307            if name not in new_covobs:
+1308                new_samples.append(new_deltas[name])
+1309                new_idl.append(new_idl_d[name])
+1310                new_means.append(new_r_values[name][i_val])
+1311                new_names_obs.append(name)
+1312        final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
+1313        for name in new_covobs:
+1314            final_result[i_val].names.append(name)
+1315        final_result[i_val]._covobs = new_covobs
+1316        final_result[i_val]._value = new_val
+1317        final_result[i_val].reweighted = reweighted
+1318
+1319    if multi == 0:
+1320        final_result = final_result.item()
+1321
+1322    return final_result
 

 
 
@@ -5328,46 +5344,46 @@ functions. For the ratio of two observables one can e.g. use
- 
1346def reweight(weight, obs, **kwargs):
-1347    """Reweight a list of observables.
-1348
-1349    Parameters
-1350    ----------
-1351    weight : Obs
-1352        Reweighting factor. An Observable that has to be defined on a superset of the
-1353        configurations in obs[i].idl for all i.
-1354    obs : list
-1355        list of Obs, e.g. [obs1, obs2, obs3].
-1356    all_configs : bool
-1357        if True, the reweighted observables are normalized by the average of
-1358        the reweighting factor on all configurations in weight.idl and not
-1359        on the configurations in obs[i].idl. Default False.
-1360    """
-1361    result = []
-1362    for i in range(len(obs)):
-1363        if len(obs[i].cov_names):
-1364            raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
-1365        if not set(obs[i].names).issubset(weight.names):
-1366            raise Exception('Error: Ensembles do not fit')
-1367        for name in obs[i].names:
-1368            if not set(obs[i].idl[name]).issubset(weight.idl[name]):
-1369                raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
-1370        new_samples = []
-1371        w_deltas = {}
-1372        for name in sorted(obs[i].names):
-1373            w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
-1374            new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
-1375        tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
-1376
-1377        if kwargs.get('all_configs'):
-1378            new_weight = weight
-1379        else:
-1380            new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
-1381
-1382        result.append(tmp_obs / new_weight)
-1383        result[-1].reweighted = True
+            
1354def reweight(weight, obs, **kwargs):
+1355    """Reweight a list of observables.
+1356
+1357    Parameters
+1358    ----------
+1359    weight : Obs
+1360        Reweighting factor. An Observable that has to be defined on a superset of the
+1361        configurations in obs[i].idl for all i.
+1362    obs : list
+1363        list of Obs, e.g. [obs1, obs2, obs3].
+1364    all_configs : bool
+1365        if True, the reweighted observables are normalized by the average of
+1366        the reweighting factor on all configurations in weight.idl and not
+1367        on the configurations in obs[i].idl. Default False.
+1368    """
+1369    result = []
+1370    for i in range(len(obs)):
+1371        if len(obs[i].cov_names):
+1372            raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
+1373        if not set(obs[i].names).issubset(weight.names):
+1374            raise Exception('Error: Ensembles do not fit')
+1375        for name in obs[i].names:
+1376            if not set(obs[i].idl[name]).issubset(weight.idl[name]):
+1377                raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
+1378        new_samples = []
+1379        w_deltas = {}
+1380        for name in sorted(obs[i].names):
+1381            w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
+1382            new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
+1383        tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
 1384
-1385    return result
+1385        if kwargs.get('all_configs'):
+1386            new_weight = weight
+1387        else:
+1388            new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
+1389
+1390        result.append(tmp_obs / new_weight)
+1391        result[-1].reweighted = True
+1392
+1393    return result
 

 
 
@@ -5401,47 +5417,47 @@ on the configurations in obs[i].idl. Default False.
- 
1388def correlate(obs_a, obs_b):
-1389    """Correlate two observables.
-1390
-1391    Parameters
-1392    ----------
-1393    obs_a : Obs
-1394        First observable
-1395    obs_b : Obs
-1396        Second observable
-1397
-1398    Notes
-1399    -----
-1400    Keep in mind to only correlate primary observables which have not been reweighted
-1401    yet. The reweighting has to be applied after correlating the observables.
-1402    Currently only works if ensembles are identical (this is not strictly necessary).
-1403    """
-1404
-1405    if sorted(obs_a.names) != sorted(obs_b.names):
-1406        raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
-1407    if len(obs_a.cov_names) or len(obs_b.cov_names):
-1408        raise Exception('Error: Not possible to correlate Obs that contain covobs!')
-1409    for name in obs_a.names:
-1410        if obs_a.shape[name] != obs_b.shape[name]:
-1411            raise Exception('Shapes of ensemble', name, 'do not fit')
-1412        if obs_a.idl[name] != obs_b.idl[name]:
-1413            raise Exception('idl of ensemble', name, 'do not fit')
-1414
-1415    if obs_a.reweighted is True:
-1416        warnings.warn("The first observable is already reweighted.", RuntimeWarning)
-1417    if obs_b.reweighted is True:
-1418        warnings.warn("The second observable is already reweighted.", RuntimeWarning)
-1419
-1420    new_samples = []
-1421    new_idl = []
-1422    for name in sorted(obs_a.names):
-1423        new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
-1424        new_idl.append(obs_a.idl[name])
-1425
-1426    o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
-1427    o.reweighted = obs_a.reweighted or obs_b.reweighted
-1428    return o
+            
1396def correlate(obs_a, obs_b):
+1397    """Correlate two observables.
+1398
+1399    Parameters
+1400    ----------
+1401    obs_a : Obs
+1402        First observable
+1403    obs_b : Obs
+1404        Second observable
+1405
+1406    Notes
+1407    -----
+1408    Keep in mind to only correlate primary observables which have not been reweighted
+1409    yet. The reweighting has to be applied after correlating the observables.
+1410    Currently only works if ensembles are identical (this is not strictly necessary).
+1411    """
+1412
+1413    if sorted(obs_a.names) != sorted(obs_b.names):
+1414        raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
+1415    if len(obs_a.cov_names) or len(obs_b.cov_names):
+1416        raise Exception('Error: Not possible to correlate Obs that contain covobs!')
+1417    for name in obs_a.names:
+1418        if obs_a.shape[name] != obs_b.shape[name]:
+1419            raise Exception('Shapes of ensemble', name, 'do not fit')
+1420        if obs_a.idl[name] != obs_b.idl[name]:
+1421            raise Exception('idl of ensemble', name, 'do not fit')
+1422
+1423    if obs_a.reweighted is True:
+1424        warnings.warn("The first observable is already reweighted.", RuntimeWarning)
+1425    if obs_b.reweighted is True:
+1426        warnings.warn("The second observable is already reweighted.", RuntimeWarning)
+1427
+1428    new_samples = []
+1429    new_idl = []
+1430    for name in sorted(obs_a.names):
+1431        new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
+1432        new_idl.append(obs_a.idl[name])
+1433
+1434    o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
+1435    o.reweighted = obs_a.reweighted or obs_b.reweighted
+1436    return o
 

 
 
@@ -5476,74 +5492,74 @@ Currently only works if ensembles are identical (this is not strictly necessary)
- 
1431def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
-1432    r'''Calculates the error covariance matrix of a set of observables.
-1433
-1434    WARNING: This function should be used with care, especially for observables with support on multiple
-1435             ensembles with differing autocorrelations. See the notes below for details.
-1436
-1437    The gamma method has to be applied first to all observables.
-1438
-1439    Parameters
-1440    ----------
-1441    obs : list or numpy.ndarray
-1442        List or one dimensional array of Obs
-1443    visualize : bool
-1444        If True plots the corresponding normalized correlation matrix (default False).
-1445    correlation : bool
-1446        If True the correlation matrix instead of the error covariance matrix is returned (default False).
-1447    smooth : None or int
-1448        If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
-1449        smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
-1450        largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
-1451        small ones.
-1452
-1453    Notes
-1454    -----
-1455    The error covariance is defined such that it agrees with the squared standard error for two identical observables
-1456    $$\operatorname{cov}(a,a)=\sum_{s=1}^N\delta_a^s\delta_a^s/N^2=\Gamma_{aa}(0)/N=\operatorname{var}(a)/N=\sigma_a^2$$
-1457    in the absence of autocorrelation.
-1458    The error covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. The covariance at windowsize 0 is guaranteed to be positive semi-definite
-1459    $$\sum_{i,j}v_i\Gamma_{ij}(0)v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i,j}v_i\delta_i^s\delta_j^s v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i}|v_i\delta_i^s|^2\geq 0\,,$$ for every $v\in\mathbb{R}^M$, while such an identity does not hold for larger windows/lags.
-1460    For observables defined on a single ensemble our approximation is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
-1461    $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
-1462    This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
-1463    '''
-1464
-1465    length = len(obs)
-1466
-1467    max_samples = np.max([o.N for o in obs])
-1468    if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
-1469        warnings.warn(f"The dimension of the covariance matrix ({length}) is larger or equal to the number of samples ({max_samples}). This will result in a rank deficient matrix.", RuntimeWarning)
-1470
-1471    cov = np.zeros((length, length))
-1472    for i in range(length):
-1473        for j in range(i, length):
-1474            cov[i, j] = _covariance_element(obs[i], obs[j])
-1475    cov = cov + cov.T - np.diag(np.diag(cov))
-1476
-1477    corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+            
1439def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
+1440    r'''Calculates the error covariance matrix of a set of observables.
+1441
+1442    WARNING: This function should be used with care, especially for observables with support on multiple
+1443             ensembles with differing autocorrelations. See the notes below for details.
+1444
+1445    The gamma method has to be applied first to all observables.
+1446
+1447    Parameters
+1448    ----------
+1449    obs : list or numpy.ndarray
+1450        List or one dimensional array of Obs
+1451    visualize : bool
+1452        If True plots the corresponding normalized correlation matrix (default False).
+1453    correlation : bool
+1454        If True the correlation matrix instead of the error covariance matrix is returned (default False).
+1455    smooth : None or int
+1456        If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
+1457        smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
+1458        largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
+1459        small ones.
+1460
+1461    Notes
+1462    -----
+1463    The error covariance is defined such that it agrees with the squared standard error for two identical observables
+1464    $$\operatorname{cov}(a,a)=\sum_{s=1}^N\delta_a^s\delta_a^s/N^2=\Gamma_{aa}(0)/N=\operatorname{var}(a)/N=\sigma_a^2$$
+1465    in the absence of autocorrelation.
+1466    The error covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. The covariance at windowsize 0 is guaranteed to be positive semi-definite
+1467    $$\sum_{i,j}v_i\Gamma_{ij}(0)v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i,j}v_i\delta_i^s\delta_j^s v_j=\frac{1}{N}\sum_{s=1}^N\sum_{i}|v_i\delta_i^s|^2\geq 0\,,$$ for every $v\in\mathbb{R}^M$, while such an identity does not hold for larger windows/lags.
+1468    For observables defined on a single ensemble our approximation is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
+1469    $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
+1470    This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
+1471    '''
+1472
+1473    length = len(obs)
+1474
+1475    max_samples = np.max([o.N for o in obs])
+1476    if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
+1477        warnings.warn(f"The dimension of the covariance matrix ({length}) is larger or equal to the number of samples ({max_samples}). This will result in a rank deficient matrix.", RuntimeWarning)
 1478
-1479    if isinstance(smooth, int):
-1480        corr = _smooth_eigenvalues(corr, smooth)
-1481
-1482    if visualize:
-1483        plt.matshow(corr, vmin=-1, vmax=1)
-1484        plt.set_cmap('RdBu')
-1485        plt.colorbar()
-1486        plt.draw()
-1487
-1488    if correlation is True:
-1489        return corr
-1490
-1491    errors = [o.dvalue for o in obs]
-1492    cov = np.diag(errors) @ corr @ np.diag(errors)
-1493
-1494    eigenvalues = np.linalg.eigh(cov)[0]
-1495    if not np.all(eigenvalues >= 0):
-1496        warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
-1497
-1498    return cov
+1479    cov = np.zeros((length, length))
+1480    for i in range(length):
+1481        for j in range(i, length):
+1482            cov[i, j] = _covariance_element(obs[i], obs[j])
+1483    cov = cov + cov.T - np.diag(np.diag(cov))
+1484
+1485    corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+1486
+1487    if isinstance(smooth, int):
+1488        corr = _smooth_eigenvalues(corr, smooth)
+1489
+1490    if visualize:
+1491        plt.matshow(corr, vmin=-1, vmax=1)
+1492        plt.set_cmap('RdBu')
+1493        plt.colorbar()
+1494        plt.draw()
+1495
+1496    if correlation is True:
+1497        return corr
+1498
+1499    errors = [o.dvalue for o in obs]
+1500    cov = np.diag(errors) @ corr @ np.diag(errors)
+1501
+1502    eigenvalues = np.linalg.eigh(cov)[0]
+1503    if not np.all(eigenvalues >= 0):
+1504        warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
+1505
+1506    return cov
 

 
 
@@ -5595,24 +5611,24 @@ This construction ensures that the estimated covariance matrix is positive semi-
- 
1578def import_jackknife(jacks, name, idl=None):
-1579    """Imports jackknife samples and returns an Obs
-1580
-1581    Parameters
-1582    ----------
-1583    jacks : numpy.ndarray
-1584        numpy array containing the mean value as zeroth entry and
-1585        the N jackknife samples as first to Nth entry.
-1586    name : str
-1587        name of the ensemble the samples are defined on.
-1588    """
-1589    length = len(jacks) - 1
-1590    prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
-1591    samples = jacks[1:] @ prj
-1592    mean = np.mean(samples)
-1593    new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
-1594    new_obs._value = jacks[0]
-1595    return new_obs
+            
1586def import_jackknife(jacks, name, idl=None):
+1587    """Imports jackknife samples and returns an Obs
+1588
+1589    Parameters
+1590    ----------
+1591    jacks : numpy.ndarray
+1592        numpy array containing the mean value as zeroth entry and
+1593        the N jackknife samples as first to Nth entry.
+1594    name : str
+1595        name of the ensemble the samples are defined on.
+1596    """
+1597    length = len(jacks) - 1
+1598    prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
+1599    samples = jacks[1:] @ prj
+1600    mean = np.mean(samples)
+1601    new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
+1602    new_obs._value = jacks[0]
+1603    return new_obs
 

 
 
@@ -5642,34 +5658,34 @@ name of the ensemble the samples are defined on.
- 
1598def import_bootstrap(boots, name, random_numbers):
-1599    """Imports bootstrap samples and returns an Obs
-1600
-1601    Parameters
-1602    ----------
-1603    boots : numpy.ndarray
-1604        numpy array containing the mean value as zeroth entry and
-1605        the N bootstrap samples as first to Nth entry.
-1606    name : str
-1607        name of the ensemble the samples are defined on.
-1608    random_numbers : np.ndarray
-1609        Array of shape (samples, length) containing the random numbers to generate the bootstrap samples,
-1610        where samples is the number of bootstrap samples and length is the length of the original Monte Carlo
-1611        chain to be reconstructed.
-1612    """
-1613    samples, length = random_numbers.shape
-1614    if samples != len(boots) - 1:
-1615        raise ValueError("Random numbers do not have the correct shape.")
-1616
-1617    if samples < length:
-1618        raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.")
-1619
-1620    proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
-1621
-1622    samples = scipy.linalg.lstsq(proj, boots[1:])[0]
-1623    ret = Obs([samples], [name])
-1624    ret._value = boots[0]
-1625    return ret
+            
1606def import_bootstrap(boots, name, random_numbers):
+1607    """Imports bootstrap samples and returns an Obs
+1608
+1609    Parameters
+1610    ----------
+1611    boots : numpy.ndarray
+1612        numpy array containing the mean value as zeroth entry and
+1613        the N bootstrap samples as first to Nth entry.
+1614    name : str
+1615        name of the ensemble the samples are defined on.
+1616    random_numbers : np.ndarray
+1617        Array of shape (samples, length) containing the random numbers to generate the bootstrap samples,
+1618        where samples is the number of bootstrap samples and length is the length of the original Monte Carlo
+1619        chain to be reconstructed.
+1620    """
+1621    samples, length = random_numbers.shape
+1622    if samples != len(boots) - 1:
+1623        raise ValueError("Random numbers do not have the correct shape.")
+1624
+1625    if samples < length:
+1626        raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.")
+1627
+1628    proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
+1629
+1630    samples = scipy.linalg.lstsq(proj, boots[1:])[0]
+1631    ret = Obs([samples], [name])
+1632    ret._value = boots[0]
+1633    return ret
 

 
 
@@ -5703,34 +5719,34 @@ chain to be reconstructed.
- 
1628def merge_obs(list_of_obs):
-1629    """Combine all observables in list_of_obs into one new observable
-1630
-1631    Parameters
-1632    ----------
-1633    list_of_obs : list
-1634        list of the Obs object to be combined
-1635
-1636    Notes
-1637    -----
-1638    It is not possible to combine obs which are based on the same replicum
-1639    """
-1640    replist = [item for obs in list_of_obs for item in obs.names]
-1641    if (len(replist) == len(set(replist))) is False:
-1642        raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
-1643    if any([len(o.cov_names) for o in list_of_obs]):
-1644        raise Exception('Not possible to merge data that contains covobs!')
-1645    new_dict = {}
-1646    idl_dict = {}
-1647    for o in list_of_obs:
-1648        new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
-1649                        for key in set(o.deltas) | set(o.r_values)})
-1650        idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
-1651
-1652    names = sorted(new_dict.keys())
-1653    o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
-1654    o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
-1655    return o
+            
1636def merge_obs(list_of_obs):
+1637    """Combine all observables in list_of_obs into one new observable
+1638
+1639    Parameters
+1640    ----------
+1641    list_of_obs : list
+1642        list of the Obs object to be combined
+1643
+1644    Notes
+1645    -----
+1646    It is not possible to combine obs which are based on the same replicum
+1647    """
+1648    replist = [item for obs in list_of_obs for item in obs.names]
+1649    if (len(replist) == len(set(replist))) is False:
+1650        raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
+1651    if any([len(o.cov_names) for o in list_of_obs]):
+1652        raise Exception('Not possible to merge data that contains covobs!')
+1653    new_dict = {}
+1654    idl_dict = {}
+1655    for o in list_of_obs:
+1656        new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
+1657                        for key in set(o.deltas) | set(o.r_values)})
+1658        idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
+1659
+1660    names = sorted(new_dict.keys())
+1661    o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
+1662    o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
+1663    return o
 

 
 
@@ -5761,47 +5777,47 @@ list of the Obs object to be combined
- 
1658def cov_Obs(means, cov, name, grad=None):
-1659    """Create an Obs based on mean(s) and a covariance matrix
-1660
-1661    Parameters
-1662    ----------
-1663    mean : list of floats or float
-1664        N mean value(s) of the new Obs
-1665    cov : list or array
-1666        2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
-1667    name : str
-1668        identifier for the covariance matrix
-1669    grad : list or array
-1670        Gradient of the Covobs wrt. the means belonging to cov.
-1671    """
-1672
-1673    def covobs_to_obs(co):
-1674        """Make an Obs out of a Covobs
-1675
-1676        Parameters
-1677        ----------
-1678        co : Covobs
-1679            Covobs to be embedded into the Obs
-1680        """
-1681        o = Obs([], [], means=[])
-1682        o._value = co.value
-1683        o.names.append(co.name)
-1684        o._covobs[co.name] = co
-1685        o._dvalue = np.sqrt(co.errsq())
-1686        return o
-1687
-1688    ol = []
-1689    if isinstance(means, (float, int)):
-1690        means = [means]
-1691
-1692    for i in range(len(means)):
-1693        ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
-1694    if ol[0].covobs[name].N != len(means):
-1695        raise Exception('You have to provide %d mean values!' % (ol[0].N))
-1696    if len(ol) == 1:
-1697        return ol[0]
-1698    return ol
+            
1666def cov_Obs(means, cov, name, grad=None):
+1667    """Create an Obs based on mean(s) and a covariance matrix
+1668
+1669    Parameters
+1670    ----------
+1671    mean : list of floats or float
+1672        N mean value(s) of the new Obs
+1673    cov : list or array
+1674        2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
+1675    name : str
+1676        identifier for the covariance matrix
+1677    grad : list or array
+1678        Gradient of the Covobs wrt. the means belonging to cov.
+1679    """
+1680
+1681    def covobs_to_obs(co):
+1682        """Make an Obs out of a Covobs
+1683
+1684        Parameters
+1685        ----------
+1686        co : Covobs
+1687            Covobs to be embedded into the Obs
+1688        """
+1689        o = Obs([], [], means=[])
+1690        o._value = co.value
+1691        o.names.append(co.name)
+1692        o._covobs[co.name] = co
+1693        o._dvalue = np.sqrt(co.errsq())
+1694        return o
+1695
+1696    ol = []
+1697    if isinstance(means, (float, int)):
+1698        means = [means]
+1699
+1700    for i in range(len(means)):
+1701        ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
+1702    if ol[0].covobs[name].N != len(means):
+1703        raise Exception('You have to provide %d mean values!' % (ol[0].N))
+1704    if len(ol) == 1:
+1705        return ol[0]
+1706    return ol