From c990fa5e01f052f65f1a82a645fb08da0e2a0486 Mon Sep 17 00:00:00 2001 From: fjosw Date: Tue, 26 Nov 2024 16:53:12 +0000 Subject: [PATCH] Documentation updated --- docs/pyerrors/obs.html | 3316 ++++++++++++++++++++-------------------- 1 file changed, 1655 insertions(+), 1661 deletions(-) diff --git a/docs/pyerrors/obs.html b/docs/pyerrors/obs.html index 57a50f3f..36510f1f 100644 --- a/docs/pyerrors/obs.html +++ b/docs/pyerrors/obs.html @@ -1192,1008 +1192,1005 @@ 856 857 def __pow__(self, y): 858 if isinstance(y, Obs): - 859 return derived_observable(lambda x: x[0] ** x[1], [self, y]) + 859 return derived_observable(lambda x, **kwargs: x[0] ** x[1], [self, y], man_grad=[y.value * self.value ** (y.value - 1), self.value ** y.value * np.log(self.value)]) 860 else: - 861 return derived_observable(lambda x: x[0] ** y, [self]) + 861 return derived_observable(lambda x, **kwargs: x[0] ** y, [self], man_grad=[y * self.value ** (y - 1)]) 862 863 def __rpow__(self, y): - 864 if isinstance(y, Obs): - 865 return derived_observable(lambda x: x[0] ** x[1], [y, self]) - 866 else: - 867 return derived_observable(lambda x: y ** x[0], [self]) + 864 return derived_observable(lambda x, **kwargs: y ** x[0], [self], man_grad=[y ** self.value * np.log(y)]) + 865 + 866 def __abs__(self): + 867 return derived_observable(lambda x: anp.abs(x[0]), [self]) 868 - 869 def __abs__(self): - 870 return derived_observable(lambda x: anp.abs(x[0]), [self]) - 871 - 872 # Overload numpy functions - 873 def sqrt(self): - 874 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)]) + 869 # Overload numpy functions + 870 def sqrt(self): + 871 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)]) + 872 + 873 def log(self): + 874 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value]) 875 - 876 def log(self): - 877 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value]) + 876 def exp(self): + 877 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)]) 878 - 879 def exp(self): - 880 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)]) + 879 def sin(self): + 880 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)]) 881 - 882 def sin(self): - 883 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)]) + 882 def cos(self): + 883 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)]) 884 - 885 def cos(self): - 886 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)]) + 885 def tan(self): + 886 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2]) 887 - 888 def tan(self): - 889 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2]) + 888 def arcsin(self): + 889 return derived_observable(lambda x: anp.arcsin(x[0]), [self]) 890 - 891 def arcsin(self): - 892 return derived_observable(lambda x: anp.arcsin(x[0]), [self]) + 891 def arccos(self): + 892 return derived_observable(lambda x: anp.arccos(x[0]), [self]) 893 - 894 def arccos(self): - 895 return derived_observable(lambda x: anp.arccos(x[0]), [self]) + 894 def arctan(self): + 895 return derived_observable(lambda x: anp.arctan(x[0]), [self]) 896 - 897 def arctan(self): - 898 return derived_observable(lambda x: anp.arctan(x[0]), [self]) + 897 def sinh(self): + 898 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)]) 899 - 900 def sinh(self): - 901 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)]) + 900 def cosh(self): + 901 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)]) 902 - 903 def cosh(self): - 904 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)]) + 903 def tanh(self): + 904 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2]) 905 - 906 def tanh(self): - 907 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2]) + 906 def arcsinh(self): + 907 return derived_observable(lambda x: anp.arcsinh(x[0]), [self]) 908 - 909 def arcsinh(self): - 910 return derived_observable(lambda x: anp.arcsinh(x[0]), [self]) + 909 def arccosh(self): + 910 return derived_observable(lambda x: anp.arccosh(x[0]), [self]) 911 - 912 def arccosh(self): - 913 return derived_observable(lambda x: anp.arccosh(x[0]), [self]) + 912 def arctanh(self): + 913 return derived_observable(lambda x: anp.arctanh(x[0]), [self]) 914 - 915 def arctanh(self): - 916 return derived_observable(lambda x: anp.arctanh(x[0]), [self]) - 917 - 918 - 919class CObs: - 920 """Class for a complex valued observable.""" - 921 __slots__ = ['_real', '_imag', 'tag'] - 922 - 923 def __init__(self, real, imag=0.0): - 924 self._real = real - 925 self._imag = imag - 926 self.tag = None - 927 - 928 @property - 929 def real(self): - 930 return self._real - 931 - 932 @property - 933 def imag(self): - 934 return self._imag - 935 - 936 def gamma_method(self, **kwargs): - 937 """Executes the gamma_method for the real and the imaginary part.""" - 938 if isinstance(self.real, Obs): - 939 self.real.gamma_method(**kwargs) - 940 if isinstance(self.imag, Obs): - 941 self.imag.gamma_method(**kwargs) - 942 - 943 def is_zero(self): - 944 """Checks whether both real and imaginary part are zero within machine precision.""" - 945 return self.real == 0.0 and self.imag == 0.0 + 915 + 916class CObs: + 917 """Class for a complex valued observable.""" + 918 __slots__ = ['_real', '_imag', 'tag'] + 919 + 920 def __init__(self, real, imag=0.0): + 921 self._real = real + 922 self._imag = imag + 923 self.tag = None + 924 + 925 @property + 926 def real(self): + 927 return self._real + 928 + 929 @property + 930 def imag(self): + 931 return self._imag + 932 + 933 def gamma_method(self, **kwargs): + 934 """Executes the gamma_method for the real and the imaginary part.""" + 935 if isinstance(self.real, Obs): + 936 self.real.gamma_method(**kwargs) + 937 if isinstance(self.imag, Obs): + 938 self.imag.gamma_method(**kwargs) + 939 + 940 def is_zero(self): + 941 """Checks whether both real and imaginary part are zero within machine precision.""" + 942 return self.real == 0.0 and self.imag == 0.0 + 943 + 944 def conjugate(self): + 945 return CObs(self.real, -self.imag) 946 - 947 def conjugate(self): - 948 return CObs(self.real, -self.imag) - 949 - 950 def __add__(self, other): - 951 if isinstance(other, np.ndarray): - 952 return other + self - 953 elif hasattr(other, 'real') and hasattr(other, 'imag'): - 954 return CObs(self.real + other.real, - 955 self.imag + other.imag) - 956 else: - 957 return CObs(self.real + other, self.imag) + 947 def __add__(self, other): + 948 if isinstance(other, np.ndarray): + 949 return other + self + 950 elif hasattr(other, 'real') and hasattr(other, 'imag'): + 951 return CObs(self.real + other.real, + 952 self.imag + other.imag) + 953 else: + 954 return CObs(self.real + other, self.imag) + 955 + 956 def __radd__(self, y): + 957 return self + y 958 - 959 def __radd__(self, y): - 960 return self + y - 961 - 962 def __sub__(self, other): - 963 if isinstance(other, np.ndarray): - 964 return -1 * (other - self) - 965 elif hasattr(other, 'real') and hasattr(other, 'imag'): - 966 return CObs(self.real - other.real, self.imag - other.imag) - 967 else: - 968 return CObs(self.real - other, self.imag) + 959 def __sub__(self, other): + 960 if isinstance(other, np.ndarray): + 961 return -1 * (other - self) + 962 elif hasattr(other, 'real') and hasattr(other, 'imag'): + 963 return CObs(self.real - other.real, self.imag - other.imag) + 964 else: + 965 return CObs(self.real - other, self.imag) + 966 + 967 def __rsub__(self, other): + 968 return -1 * (self - other) 969 - 970 def __rsub__(self, other): - 971 return -1 * (self - other) - 972 - 973 def __mul__(self, other): - 974 if isinstance(other, np.ndarray): - 975 return other * self - 976 elif hasattr(other, 'real') and hasattr(other, 'imag'): - 977 if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]): - 978 return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3], + 970 def __mul__(self, other): + 971 if isinstance(other, np.ndarray): + 972 return other * self + 973 elif hasattr(other, 'real') and hasattr(other, 'imag'): + 974 if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]): + 975 return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3], + 976 [self.real, other.real, self.imag, other.imag], + 977 man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]), + 978 derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3], 979 [self.real, other.real, self.imag, other.imag], - 980 man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]), - 981 derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3], - 982 [self.real, other.real, self.imag, other.imag], - 983 man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value])) - 984 elif getattr(other, 'imag', 0) != 0: - 985 return CObs(self.real * other.real - self.imag * other.imag, - 986 self.imag * other.real + self.real * other.imag) - 987 else: - 988 return CObs(self.real * other.real, self.imag * other.real) - 989 else: - 990 return CObs(self.real * other, self.imag * other) + 980 man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value])) + 981 elif getattr(other, 'imag', 0) != 0: + 982 return CObs(self.real * other.real - self.imag * other.imag, + 983 self.imag * other.real + self.real * other.imag) + 984 else: + 985 return CObs(self.real * other.real, self.imag * other.real) + 986 else: + 987 return CObs(self.real * other, self.imag * other) + 988 + 989 def __rmul__(self, other): + 990 return self * other 991 - 992 def __rmul__(self, other): - 993 return self * other - 994 - 995 def __truediv__(self, other): - 996 if isinstance(other, np.ndarray): - 997 return 1 / (other / self) - 998 elif hasattr(other, 'real') and hasattr(other, 'imag'): - 999 r = other.real ** 2 + other.imag ** 2 -1000 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r) -1001 else: -1002 return CObs(self.real / other, self.imag / other) -1003 -1004 def __rtruediv__(self, other): -1005 r = self.real ** 2 + self.imag ** 2 -1006 if hasattr(other, 'real') and hasattr(other, 'imag'): -1007 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r) -1008 else: -1009 return CObs(self.real * other / r, -self.imag * other / r) + 992 def __truediv__(self, other): + 993 if isinstance(other, np.ndarray): + 994 return 1 / (other / self) + 995 elif hasattr(other, 'real') and hasattr(other, 'imag'): + 996 r = other.real ** 2 + other.imag ** 2 + 997 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r) + 998 else: + 999 return CObs(self.real / other, self.imag / other) +1000 +1001 def __rtruediv__(self, other): +1002 r = self.real ** 2 + self.imag ** 2 +1003 if hasattr(other, 'real') and hasattr(other, 'imag'): +1004 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r) +1005 else: +1006 return CObs(self.real * other / r, -self.imag * other / r) +1007 +1008 def __abs__(self): +1009 return np.sqrt(self.real**2 + self.imag**2) 1010 -1011 def __abs__(self): -1012 return np.sqrt(self.real**2 + self.imag**2) +1011 def __pos__(self): +1012 return self 1013 -1014 def __pos__(self): -1015 return self +1014 def __neg__(self): +1015 return -1 * self 1016 -1017 def __neg__(self): -1018 return -1 * self +1017 def __eq__(self, other): +1018 return self.real == other.real and self.imag == other.imag 1019 -1020 def __eq__(self, other): -1021 return self.real == other.real and self.imag == other.imag +1020 def __str__(self): +1021 return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)' 1022 -1023 def __str__(self): -1024 return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)' +1023 def __repr__(self): +1024 return 'CObs[' + str(self) + ']' 1025 -1026 def __repr__(self): -1027 return 'CObs[' + str(self) + ']' -1028 -1029 def __format__(self, format_type): -1030 if format_type == "": -1031 significance = 2 -1032 format_type = "2" -1033 else: -1034 significance = int(float(format_type.replace("+", "").replace("-", ""))) -1035 return f"({self.real:{format_type}}{self.imag:+{significance}}j)" -1036 +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 gamma_method(x, **kwargs): +1036 """Vectorized version of the gamma_method applicable to lists or arrays of Obs. 1037 -1038def gamma_method(x, **kwargs): -1039 """Vectorized version of the gamma_method applicable to lists or arrays of Obs. -1040 -1041 See docstring of pe.Obs.gamma_method for details. -1042 """ -1043 return np.vectorize(lambda o: o.gm(**kwargs))(x) +1038 See docstring of pe.Obs.gamma_method for details. +1039 """ +1040 return np.vectorize(lambda o: o.gm(**kwargs))(x) +1041 +1042 +1043gm = gamma_method 1044 1045 -1046gm = gamma_method -1047 -1048 -1049def _format_uncertainty(value, dvalue, significance=2): -1050 """Creates a string of a value and its error in paranthesis notation, e.g., 13.02(45)""" -1051 if dvalue == 0.0 or (not np.isfinite(dvalue)): -1052 return str(value) -1053 if not isinstance(significance, int): -1054 raise TypeError("significance needs to be an integer.") -1055 if significance < 1: -1056 raise ValueError("significance needs to be larger than zero.") -1057 fexp = np.floor(np.log10(dvalue)) -1058 if fexp < 0.0: -1059 return '{:{form}}({:1.0f})'.format(value, dvalue * 10 ** (-fexp + significance - 1), form='.' + str(-int(fexp) + significance - 1) + 'f') -1060 elif fexp == 0.0: -1061 return f"{value:.{significance - 1}f}({dvalue:1.{significance - 1}f})" -1062 else: -1063 return f"{value:.{max(0, int(significance - fexp - 1))}f}({dvalue:2.{max(0, int(significance - fexp - 1))}f})" -1064 -1065 -1066def _expand_deltas(deltas, idx, shape, gapsize): -1067 """Expand deltas defined on idx to a regular range with spacing gapsize between two -1068 configurations and where holes are filled by 0. -1069 If idx is of type range, the deltas are not changed if the idx.step == gapsize. -1070 -1071 Parameters -1072 ---------- -1073 deltas : list -1074 List of fluctuations -1075 idx : list -1076 List or range of configs on which the deltas are defined, has to be sorted in ascending order. -1077 shape : int -1078 Number of configs in idx. -1079 gapsize : int -1080 The target distance between two configurations. If longer distances -1081 are found in idx, the data is expanded. -1082 """ -1083 if isinstance(idx, range): -1084 if (idx.step == gapsize): -1085 return deltas -1086 ret = np.zeros((idx[-1] - idx[0] + gapsize) // gapsize) -1087 for i in range(shape): -1088 ret[(idx[i] - idx[0]) // gapsize] = deltas[i] -1089 return ret -1090 +1046def _format_uncertainty(value, dvalue, significance=2): +1047 """Creates a string of a value and its error in paranthesis notation, e.g., 13.02(45)""" +1048 if dvalue == 0.0 or (not np.isfinite(dvalue)): +1049 return str(value) +1050 if not isinstance(significance, int): +1051 raise TypeError("significance needs to be an integer.") +1052 if significance < 1: +1053 raise ValueError("significance needs to be larger than zero.") +1054 fexp = np.floor(np.log10(dvalue)) +1055 if fexp < 0.0: +1056 return '{:{form}}({:1.0f})'.format(value, dvalue * 10 ** (-fexp + significance - 1), form='.' + str(-int(fexp) + significance - 1) + 'f') +1057 elif fexp == 0.0: +1058 return f"{value:.{significance - 1}f}({dvalue:1.{significance - 1}f})" +1059 else: +1060 return f"{value:.{max(0, int(significance - fexp - 1))}f}({dvalue:2.{max(0, int(significance - fexp - 1))}f})" +1061 +1062 +1063def _expand_deltas(deltas, idx, shape, gapsize): +1064 """Expand deltas defined on idx to a regular range with spacing gapsize between two +1065 configurations and where holes are filled by 0. +1066 If idx is of type range, the deltas are not changed if the idx.step == gapsize. +1067 +1068 Parameters +1069 ---------- +1070 deltas : list +1071 List of fluctuations +1072 idx : list +1073 List or range of configs on which the deltas are defined, has to be sorted in ascending order. +1074 shape : int +1075 Number of configs in idx. +1076 gapsize : int +1077 The target distance between two configurations. If longer distances +1078 are found in idx, the data is expanded. +1079 """ +1080 if isinstance(idx, range): +1081 if (idx.step == gapsize): +1082 return deltas +1083 ret = np.zeros((idx[-1] - idx[0] + gapsize) // gapsize) +1084 for i in range(shape): +1085 ret[(idx[i] - idx[0]) // gapsize] = deltas[i] +1086 return ret +1087 +1088 +1089def _merge_idx(idl): +1090 """Returns the union of all lists in idl as range or sorted list 1091 -1092def _merge_idx(idl): -1093 """Returns the union of all lists in idl as range or sorted list -1094 -1095 Parameters -1096 ---------- -1097 idl : list -1098 List of lists or ranges. -1099 """ +1092 Parameters +1093 ---------- +1094 idl : list +1095 List of lists or ranges. +1096 """ +1097 +1098 if _check_lists_equal(idl): +1099 return idl[0] 1100 -1101 if _check_lists_equal(idl): -1102 return idl[0] -1103 -1104 idunion = sorted(set().union(*idl)) -1105 -1106 # Check whether idunion can be expressed as range -1107 idrange = range(idunion[0], idunion[-1] + 1, idunion[1] - idunion[0]) -1108 idtest = [list(idrange), idunion] -1109 if _check_lists_equal(idtest): -1110 return idrange +1101 idunion = sorted(set().union(*idl)) +1102 +1103 # Check whether idunion can be expressed as range +1104 idrange = range(idunion[0], idunion[-1] + 1, idunion[1] - idunion[0]) +1105 idtest = [list(idrange), idunion] +1106 if _check_lists_equal(idtest): +1107 return idrange +1108 +1109 return idunion +1110 1111 -1112 return idunion -1113 +1112def _intersection_idx(idl): +1113 """Returns the intersection of all lists in idl as range or sorted list 1114 -1115def _intersection_idx(idl): -1116 """Returns the intersection of all lists in idl as range or sorted list -1117 -1118 Parameters -1119 ---------- -1120 idl : list -1121 List of lists or ranges. -1122 """ +1115 Parameters +1116 ---------- +1117 idl : list +1118 List of lists or ranges. +1119 """ +1120 +1121 if _check_lists_equal(idl): +1122 return idl[0] 1123 -1124 if _check_lists_equal(idl): -1125 return idl[0] -1126 -1127 idinter = sorted(set.intersection(*[set(o) for o in idl])) -1128 -1129 # Check whether idinter can be expressed as range -1130 try: -1131 idrange = range(idinter[0], idinter[-1] + 1, idinter[1] - idinter[0]) -1132 idtest = [list(idrange), idinter] -1133 if _check_lists_equal(idtest): -1134 return idrange -1135 except IndexError: -1136 pass +1124 idinter = sorted(set.intersection(*[set(o) for o in idl])) +1125 +1126 # Check whether idinter can be expressed as range +1127 try: +1128 idrange = range(idinter[0], idinter[-1] + 1, idinter[1] - idinter[0]) +1129 idtest = [list(idrange), idinter] +1130 if _check_lists_equal(idtest): +1131 return idrange +1132 except IndexError: +1133 pass +1134 +1135 return idinter +1136 1137 -1138 return idinter -1139 -1140 -1141def _expand_deltas_for_merge(deltas, idx, shape, new_idx, scalefactor): -1142 """Expand deltas defined on idx to the list of configs that is defined by new_idx. -1143 New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest -1144 common divisor of the step sizes is used as new step size. -1145 -1146 Parameters -1147 ---------- -1148 deltas : list -1149 List of fluctuations -1150 idx : list -1151 List or range of configs on which the deltas are defined. -1152 Has to be a subset of new_idx and has to be sorted in ascending order. -1153 shape : list -1154 Number of configs in idx. -1155 new_idx : list -1156 List of configs that defines the new range, has to be sorted in ascending order. -1157 scalefactor : float -1158 An additional scaling factor that can be applied to scale the fluctuations, -1159 e.g., when Obs with differing numbers of replica are merged. -1160 """ -1161 if type(idx) is range and type(new_idx) is range: -1162 if idx == new_idx: -1163 if scalefactor == 1: -1164 return deltas -1165 else: -1166 return deltas * scalefactor -1167 ret = np.zeros(new_idx[-1] - new_idx[0] + 1) -1168 for i in range(shape): -1169 ret[idx[i] - new_idx[0]] = deltas[i] -1170 return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx) * scalefactor -1171 +1138def _expand_deltas_for_merge(deltas, idx, shape, new_idx, scalefactor): +1139 """Expand deltas defined on idx to the list of configs that is defined by new_idx. +1140 New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest +1141 common divisor of the step sizes is used as new step size. +1142 +1143 Parameters +1144 ---------- +1145 deltas : list +1146 List of fluctuations +1147 idx : list +1148 List or range of configs on which the deltas are defined. +1149 Has to be a subset of new_idx and has to be sorted in ascending order. +1150 shape : list +1151 Number of configs in idx. +1152 new_idx : list +1153 List of configs that defines the new range, has to be sorted in ascending order. +1154 scalefactor : float +1155 An additional scaling factor that can be applied to scale the fluctuations, +1156 e.g., when Obs with differing numbers of replica are merged. +1157 """ +1158 if type(idx) is range and type(new_idx) is range: +1159 if idx == new_idx: +1160 if scalefactor == 1: +1161 return deltas +1162 else: +1163 return deltas * scalefactor +1164 ret = np.zeros(new_idx[-1] - new_idx[0] + 1) +1165 for i in range(shape): +1166 ret[idx[i] - new_idx[0]] = deltas[i] +1167 return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx) * scalefactor +1168 +1169 +1170def derived_observable(func, data, array_mode=False, **kwargs): +1171 """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation. 1172 -1173def derived_observable(func, data, array_mode=False, **kwargs): -1174 """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation. -1175 -1176 Parameters -1177 ---------- -1178 func : object -1179 arbitrary function of the form func(data, **kwargs). For the -1180 automatic differentiation to work, all numpy functions have to have -1181 the autograd wrapper (use 'import autograd.numpy as anp'). -1182 data : list -1183 list of Obs, e.g. [obs1, obs2, obs3]. -1184 num_grad : bool -1185 if True, numerical derivatives are used instead of autograd -1186 (default False). To control the numerical differentiation the -1187 kwargs of numdifftools.step_generators.MaxStepGenerator -1188 can be used. -1189 man_grad : list -1190 manually supply a list or an array which contains the jacobian -1191 of func. Use cautiously, supplying the wrong derivative will -1192 not be intercepted. -1193 -1194 Notes -1195 ----- -1196 For simple mathematical operations it can be practical to use anonymous -1197 functions. For the ratio of two observables one can e.g. use +1173 Parameters +1174 ---------- +1175 func : object +1176 arbitrary function of the form func(data, **kwargs). For the +1177 automatic differentiation to work, all numpy functions have to have +1178 the autograd wrapper (use 'import autograd.numpy as anp'). +1179 data : list +1180 list of Obs, e.g. [obs1, obs2, obs3]. +1181 num_grad : bool +1182 if True, numerical derivatives are used instead of autograd +1183 (default False). To control the numerical differentiation the +1184 kwargs of numdifftools.step_generators.MaxStepGenerator +1185 can be used. +1186 man_grad : list +1187 manually supply a list or an array which contains the jacobian +1188 of func. Use cautiously, supplying the wrong derivative will +1189 not be intercepted. +1190 +1191 Notes +1192 ----- +1193 For simple mathematical operations it can be practical to use anonymous +1194 functions. For the ratio of two observables one can e.g. use +1195 +1196 new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2]) +1197 """ 1198 -1199 new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2]) -1200 """ +1199 data = np.asarray(data) +1200 raveled_data = data.ravel() 1201 -1202 data = np.asarray(data) -1203 raveled_data = data.ravel() -1204 -1205 # Workaround for matrix operations containing non Obs data -1206 if not all(isinstance(x, Obs) for x in raveled_data): -1207 for i in range(len(raveled_data)): -1208 if isinstance(raveled_data[i], (int, float)): -1209 raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###") -1210 -1211 allcov = {} -1212 for o in raveled_data: -1213 for name in o.cov_names: -1214 if name in allcov: -1215 if not np.allclose(allcov[name], o.covobs[name].cov): -1216 raise Exception('Inconsistent covariance matrices for %s!' % (name)) -1217 else: -1218 allcov[name] = o.covobs[name].cov -1219 -1220 n_obs = len(raveled_data) -1221 new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x])) -1222 new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x])) -1223 new_sample_names = sorted(set(new_names) - set(new_cov_names)) -1224 -1225 reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0 -1226 -1227 if data.ndim == 1: -1228 values = np.array([o.value for o in data]) -1229 else: -1230 values = np.vectorize(lambda x: x.value)(data) -1231 -1232 new_values = func(values, **kwargs) -1233 -1234 multi = int(isinstance(new_values, np.ndarray)) -1235 -1236 new_r_values = {} -1237 new_idl_d = {} -1238 for name in new_sample_names: -1239 idl = [] -1240 tmp_values = np.zeros(n_obs) -1241 for i, item in enumerate(raveled_data): -1242 tmp_values[i] = item.r_values.get(name, item.value) -1243 tmp_idl = item.idl.get(name) -1244 if tmp_idl is not None: -1245 idl.append(tmp_idl) -1246 if multi > 0: -1247 tmp_values = np.array(tmp_values).reshape(data.shape) -1248 new_r_values[name] = func(tmp_values, **kwargs) -1249 new_idl_d[name] = _merge_idx(idl) -1250 -1251 def _compute_scalefactor_missing_rep(obs): -1252 """ -1253 Computes the scale factor that is to be multiplied with the deltas -1254 in the case where Obs with different subsets of replica are merged. -1255 Returns a dictionary with the scale factor for each Monte Carlo name. -1256 -1257 Parameters -1258 ---------- -1259 obs : Obs -1260 The observable corresponding to the deltas that are to be scaled -1261 """ -1262 scalef_d = {} -1263 for mc_name in obs.mc_names: -1264 mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')] -1265 new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')] -1266 if len(mc_idl_d) > 0 and len(mc_idl_d) < len(new_mc_idl_d): -1267 scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d]) -1268 return scalef_d -1269 -1270 if 'man_grad' in kwargs: -1271 deriv = np.asarray(kwargs.get('man_grad')) -1272 if new_values.shape + data.shape != deriv.shape: -1273 raise Exception('Manual derivative does not have correct shape.') -1274 elif kwargs.get('num_grad') is True: -1275 if multi > 0: -1276 raise Exception('Multi mode currently not supported for numerical derivative') -1277 options = { -1278 'base_step': 0.1, -1279 'step_ratio': 2.5} -1280 for key in options.keys(): -1281 kwarg = kwargs.get(key) -1282 if kwarg is not None: -1283 options[key] = kwarg -1284 tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs) -1285 if tmp_df.size == 1: -1286 deriv = np.array([tmp_df.real]) -1287 else: -1288 deriv = tmp_df.real -1289 else: -1290 deriv = jacobian(func)(values, **kwargs) -1291 -1292 final_result = np.zeros(new_values.shape, dtype=object) -1293 -1294 if array_mode is True: -1295 -1296 class _Zero_grad(): -1297 def __init__(self, N): -1298 self.grad = np.zeros((N, 1)) -1299 -1300 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])) -1301 d_extracted = {} -1302 g_extracted = {} -1303 for name in new_sample_names: -1304 d_extracted[name] = [] -1305 ens_length = len(new_idl_d[name]) -1306 for i_dat, dat in enumerate(data): -1307 d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name], _compute_scalefactor_missing_rep(o).get(name.split('|')[0], 1)) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, ))) -1308 for name in new_cov_names: -1309 g_extracted[name] = [] -1310 zero_grad = _Zero_grad(new_covobs_lengths[name]) -1311 for i_dat, dat in enumerate(data): -1312 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))) -1313 -1314 for i_val, new_val in np.ndenumerate(new_values): -1315 new_deltas = {} -1316 new_grad = {} -1317 if array_mode is True: -1318 for name in new_sample_names: -1319 ens_length = d_extracted[name][0].shape[-1] -1320 new_deltas[name] = np.zeros(ens_length) -1321 for i_dat, dat in enumerate(d_extracted[name]): -1322 new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat) -1323 for name in new_cov_names: -1324 new_grad[name] = 0 -1325 for i_dat, dat in enumerate(g_extracted[name]): -1326 new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat) -1327 else: -1328 for j_obs, obs in np.ndenumerate(data): -1329 scalef_d = _compute_scalefactor_missing_rep(obs) -1330 for name in obs.names: -1331 if name in obs.cov_names: -1332 new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad -1333 else: -1334 new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name], scalef_d.get(name.split('|')[0], 1)) -1335 -1336 new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad} -1337 -1338 if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()): -1339 raise Exception('The same name has been used for deltas and covobs!') -1340 new_samples = [] -1341 new_means = [] -1342 new_idl = [] -1343 new_names_obs = [] -1344 for name in new_names: -1345 if name not in new_covobs: -1346 new_samples.append(new_deltas[name]) -1347 new_idl.append(new_idl_d[name]) -1348 new_means.append(new_r_values[name][i_val]) -1349 new_names_obs.append(name) -1350 final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl) -1351 for name in new_covobs: -1352 final_result[i_val].names.append(name) -1353 final_result[i_val]._covobs = new_covobs -1354 final_result[i_val]._value = new_val -1355 final_result[i_val].reweighted = reweighted +1202 # Workaround for matrix operations containing non Obs data +1203 if not all(isinstance(x, Obs) for x in raveled_data): +1204 for i in range(len(raveled_data)): +1205 if isinstance(raveled_data[i], (int, float)): +1206 raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###") +1207 +1208 allcov = {} +1209 for o in raveled_data: +1210 for name in o.cov_names: +1211 if name in allcov: +1212 if not np.allclose(allcov[name], o.covobs[name].cov): +1213 raise Exception('Inconsistent covariance matrices for %s!' % (name)) +1214 else: +1215 allcov[name] = o.covobs[name].cov +1216 +1217 n_obs = len(raveled_data) +1218 new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x])) +1219 new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x])) +1220 new_sample_names = sorted(set(new_names) - set(new_cov_names)) +1221 +1222 reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0 +1223 +1224 if data.ndim == 1: +1225 values = np.array([o.value for o in data]) +1226 else: +1227 values = np.vectorize(lambda x: x.value)(data) +1228 +1229 new_values = func(values, **kwargs) +1230 +1231 multi = int(isinstance(new_values, np.ndarray)) +1232 +1233 new_r_values = {} +1234 new_idl_d = {} +1235 for name in new_sample_names: +1236 idl = [] +1237 tmp_values = np.zeros(n_obs) +1238 for i, item in enumerate(raveled_data): +1239 tmp_values[i] = item.r_values.get(name, item.value) +1240 tmp_idl = item.idl.get(name) +1241 if tmp_idl is not None: +1242 idl.append(tmp_idl) +1243 if multi > 0: +1244 tmp_values = np.array(tmp_values).reshape(data.shape) +1245 new_r_values[name] = func(tmp_values, **kwargs) +1246 new_idl_d[name] = _merge_idx(idl) +1247 +1248 def _compute_scalefactor_missing_rep(obs): +1249 """ +1250 Computes the scale factor that is to be multiplied with the deltas +1251 in the case where Obs with different subsets of replica are merged. +1252 Returns a dictionary with the scale factor for each Monte Carlo name. +1253 +1254 Parameters +1255 ---------- +1256 obs : Obs +1257 The observable corresponding to the deltas that are to be scaled +1258 """ +1259 scalef_d = {} +1260 for mc_name in obs.mc_names: +1261 mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')] +1262 new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')] +1263 if len(mc_idl_d) > 0 and len(mc_idl_d) < len(new_mc_idl_d): +1264 scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d]) +1265 return scalef_d +1266 +1267 if 'man_grad' in kwargs: +1268 deriv = np.asarray(kwargs.get('man_grad')) +1269 if new_values.shape + data.shape != deriv.shape: +1270 raise Exception('Manual derivative does not have correct shape.') +1271 elif kwargs.get('num_grad') is True: +1272 if multi > 0: +1273 raise Exception('Multi mode currently not supported for numerical derivative') +1274 options = { +1275 'base_step': 0.1, +1276 'step_ratio': 2.5} +1277 for key in options.keys(): +1278 kwarg = kwargs.get(key) +1279 if kwarg is not None: +1280 options[key] = kwarg +1281 tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs) +1282 if tmp_df.size == 1: +1283 deriv = np.array([tmp_df.real]) +1284 else: +1285 deriv = tmp_df.real +1286 else: +1287 deriv = jacobian(func)(values, **kwargs) +1288 +1289 final_result = np.zeros(new_values.shape, dtype=object) +1290 +1291 if array_mode is True: +1292 +1293 class _Zero_grad(): +1294 def __init__(self, N): +1295 self.grad = np.zeros((N, 1)) +1296 +1297 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])) +1298 d_extracted = {} +1299 g_extracted = {} +1300 for name in new_sample_names: +1301 d_extracted[name] = [] +1302 ens_length = len(new_idl_d[name]) +1303 for i_dat, dat in enumerate(data): +1304 d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name], _compute_scalefactor_missing_rep(o).get(name.split('|')[0], 1)) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, ))) +1305 for name in new_cov_names: +1306 g_extracted[name] = [] +1307 zero_grad = _Zero_grad(new_covobs_lengths[name]) +1308 for i_dat, dat in enumerate(data): +1309 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))) +1310 +1311 for i_val, new_val in np.ndenumerate(new_values): +1312 new_deltas = {} +1313 new_grad = {} +1314 if array_mode is True: +1315 for name in new_sample_names: +1316 ens_length = d_extracted[name][0].shape[-1] +1317 new_deltas[name] = np.zeros(ens_length) +1318 for i_dat, dat in enumerate(d_extracted[name]): +1319 new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat) +1320 for name in new_cov_names: +1321 new_grad[name] = 0 +1322 for i_dat, dat in enumerate(g_extracted[name]): +1323 new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat) +1324 else: +1325 for j_obs, obs in np.ndenumerate(data): +1326 scalef_d = _compute_scalefactor_missing_rep(obs) +1327 for name in obs.names: +1328 if name in obs.cov_names: +1329 new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad +1330 else: +1331 new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name], scalef_d.get(name.split('|')[0], 1)) +1332 +1333 new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad} +1334 +1335 if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()): +1336 raise Exception('The same name has been used for deltas and covobs!') +1337 new_samples = [] +1338 new_means = [] +1339 new_idl = [] +1340 new_names_obs = [] +1341 for name in new_names: +1342 if name not in new_covobs: +1343 new_samples.append(new_deltas[name]) +1344 new_idl.append(new_idl_d[name]) +1345 new_means.append(new_r_values[name][i_val]) +1346 new_names_obs.append(name) +1347 final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl) +1348 for name in new_covobs: +1349 final_result[i_val].names.append(name) +1350 final_result[i_val]._covobs = new_covobs +1351 final_result[i_val]._value = new_val +1352 final_result[i_val].reweighted = reweighted +1353 +1354 if multi == 0: +1355 final_result = final_result.item() 1356 -1357 if multi == 0: -1358 final_result = final_result.item() +1357 return final_result +1358 1359 -1360 return final_result -1361 +1360def _reduce_deltas(deltas, idx_old, idx_new): +1361 """Extract deltas defined on idx_old on all configs of idx_new. 1362 -1363def _reduce_deltas(deltas, idx_old, idx_new): -1364 """Extract deltas defined on idx_old on all configs of idx_new. +1363 Assumes, that idx_old and idx_new are correctly defined idl, i.e., they +1364 are ordered in an ascending order. 1365 -1366 Assumes, that idx_old and idx_new are correctly defined idl, i.e., they -1367 are ordered in an ascending order. -1368 -1369 Parameters -1370 ---------- -1371 deltas : list -1372 List of fluctuations -1373 idx_old : list -1374 List or range of configs on which the deltas are defined -1375 idx_new : list -1376 List of configs for which we want to extract the deltas. -1377 Has to be a subset of idx_old. -1378 """ -1379 if not len(deltas) == len(idx_old): -1380 raise Exception('Length of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old))) -1381 if type(idx_old) is range and type(idx_new) is range: -1382 if idx_old == idx_new: -1383 return deltas -1384 if _check_lists_equal([idx_old, idx_new]): -1385 return deltas -1386 indices = np.intersect1d(idx_old, idx_new, assume_unique=True, return_indices=True)[1] -1387 if len(indices) < len(idx_new): -1388 raise Exception('Error in _reduce_deltas: Config of idx_new not in idx_old') -1389 return np.array(deltas)[indices] -1390 +1366 Parameters +1367 ---------- +1368 deltas : list +1369 List of fluctuations +1370 idx_old : list +1371 List or range of configs on which the deltas are defined +1372 idx_new : list +1373 List of configs for which we want to extract the deltas. +1374 Has to be a subset of idx_old. +1375 """ +1376 if not len(deltas) == len(idx_old): +1377 raise Exception('Length of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old))) +1378 if type(idx_old) is range and type(idx_new) is range: +1379 if idx_old == idx_new: +1380 return deltas +1381 if _check_lists_equal([idx_old, idx_new]): +1382 return deltas +1383 indices = np.intersect1d(idx_old, idx_new, assume_unique=True, return_indices=True)[1] +1384 if len(indices) < len(idx_new): +1385 raise Exception('Error in _reduce_deltas: Config of idx_new not in idx_old') +1386 return np.array(deltas)[indices] +1387 +1388 +1389def reweight(weight, obs, **kwargs): +1390 """Reweight a list of observables. 1391 -1392def reweight(weight, obs, **kwargs): -1393 """Reweight a list of observables. -1394 -1395 Parameters -1396 ---------- -1397 weight : Obs -1398 Reweighting factor. An Observable that has to be defined on a superset of the -1399 configurations in obs[i].idl for all i. -1400 obs : list -1401 list of Obs, e.g. [obs1, obs2, obs3]. -1402 all_configs : bool -1403 if True, the reweighted observables are normalized by the average of -1404 the reweighting factor on all configurations in weight.idl and not -1405 on the configurations in obs[i].idl. Default False. -1406 """ -1407 result = [] -1408 for i in range(len(obs)): -1409 if len(obs[i].cov_names): -1410 raise Exception('Error: Not possible to reweight an Obs that contains covobs!') -1411 if not set(obs[i].names).issubset(weight.names): -1412 raise Exception('Error: Ensembles do not fit') -1413 for name in obs[i].names: -1414 if not set(obs[i].idl[name]).issubset(weight.idl[name]): -1415 raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name)) -1416 new_samples = [] -1417 w_deltas = {} -1418 for name in sorted(obs[i].names): -1419 w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name]) -1420 new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name])) -1421 tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)]) -1422 -1423 if kwargs.get('all_configs'): -1424 new_weight = weight -1425 else: -1426 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)]) +1392 Parameters +1393 ---------- +1394 weight : Obs +1395 Reweighting factor. An Observable that has to be defined on a superset of the +1396 configurations in obs[i].idl for all i. +1397 obs : list +1398 list of Obs, e.g. [obs1, obs2, obs3]. +1399 all_configs : bool +1400 if True, the reweighted observables are normalized by the average of +1401 the reweighting factor on all configurations in weight.idl and not +1402 on the configurations in obs[i].idl. Default False. +1403 """ +1404 result = [] +1405 for i in range(len(obs)): +1406 if len(obs[i].cov_names): +1407 raise Exception('Error: Not possible to reweight an Obs that contains covobs!') +1408 if not set(obs[i].names).issubset(weight.names): +1409 raise Exception('Error: Ensembles do not fit') +1410 for name in obs[i].names: +1411 if not set(obs[i].idl[name]).issubset(weight.idl[name]): +1412 raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name)) +1413 new_samples = [] +1414 w_deltas = {} +1415 for name in sorted(obs[i].names): +1416 w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name]) +1417 new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name])) +1418 tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)]) +1419 +1420 if kwargs.get('all_configs'): +1421 new_weight = weight +1422 else: +1423 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)]) +1424 +1425 result.append(tmp_obs / new_weight) +1426 result[-1].reweighted = True 1427 -1428 result.append(tmp_obs / new_weight) -1429 result[-1].reweighted = True +1428 return result +1429 1430 -1431 return result -1432 +1431def correlate(obs_a, obs_b): +1432 """Correlate two observables. 1433 -1434def correlate(obs_a, obs_b): -1435 """Correlate two observables. -1436 -1437 Parameters -1438 ---------- -1439 obs_a : Obs -1440 First observable -1441 obs_b : Obs -1442 Second observable -1443 -1444 Notes -1445 ----- -1446 Keep in mind to only correlate primary observables which have not been reweighted -1447 yet. The reweighting has to be applied after correlating the observables. -1448 Currently only works if ensembles are identical (this is not strictly necessary). -1449 """ -1450 -1451 if sorted(obs_a.names) != sorted(obs_b.names): -1452 raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}") -1453 if len(obs_a.cov_names) or len(obs_b.cov_names): -1454 raise Exception('Error: Not possible to correlate Obs that contain covobs!') -1455 for name in obs_a.names: -1456 if obs_a.shape[name] != obs_b.shape[name]: -1457 raise Exception('Shapes of ensemble', name, 'do not fit') -1458 if obs_a.idl[name] != obs_b.idl[name]: -1459 raise Exception('idl of ensemble', name, 'do not fit') -1460 -1461 if obs_a.reweighted is True: -1462 warnings.warn("The first observable is already reweighted.", RuntimeWarning) -1463 if obs_b.reweighted is True: -1464 warnings.warn("The second observable is already reweighted.", RuntimeWarning) -1465 -1466 new_samples = [] -1467 new_idl = [] -1468 for name in sorted(obs_a.names): -1469 new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name])) -1470 new_idl.append(obs_a.idl[name]) -1471 -1472 o = Obs(new_samples, sorted(obs_a.names), idl=new_idl) -1473 o.reweighted = obs_a.reweighted or obs_b.reweighted -1474 return o -1475 +1434 Parameters +1435 ---------- +1436 obs_a : Obs +1437 First observable +1438 obs_b : Obs +1439 Second observable +1440 +1441 Notes +1442 ----- +1443 Keep in mind to only correlate primary observables which have not been reweighted +1444 yet. The reweighting has to be applied after correlating the observables. +1445 Currently only works if ensembles are identical (this is not strictly necessary). +1446 """ +1447 +1448 if sorted(obs_a.names) != sorted(obs_b.names): +1449 raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}") +1450 if len(obs_a.cov_names) or len(obs_b.cov_names): +1451 raise Exception('Error: Not possible to correlate Obs that contain covobs!') +1452 for name in obs_a.names: +1453 if obs_a.shape[name] != obs_b.shape[name]: +1454 raise Exception('Shapes of ensemble', name, 'do not fit') +1455 if obs_a.idl[name] != obs_b.idl[name]: +1456 raise Exception('idl of ensemble', name, 'do not fit') +1457 +1458 if obs_a.reweighted is True: +1459 warnings.warn("The first observable is already reweighted.", RuntimeWarning) +1460 if obs_b.reweighted is True: +1461 warnings.warn("The second observable is already reweighted.", RuntimeWarning) +1462 +1463 new_samples = [] +1464 new_idl = [] +1465 for name in sorted(obs_a.names): +1466 new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name])) +1467 new_idl.append(obs_a.idl[name]) +1468 +1469 o = Obs(new_samples, sorted(obs_a.names), idl=new_idl) +1470 o.reweighted = obs_a.reweighted or obs_b.reweighted +1471 return o +1472 +1473 +1474def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs): +1475 r'''Calculates the error covariance matrix of a set of observables. 1476 -1477def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs): -1478 r'''Calculates the error covariance matrix of a set of observables. +1477 WARNING: This function should be used with care, especially for observables with support on multiple +1478 ensembles with differing autocorrelations. See the notes below for details. 1479 -1480 WARNING: This function should be used with care, especially for observables with support on multiple -1481 ensembles with differing autocorrelations. See the notes below for details. -1482 -1483 The gamma method has to be applied first to all observables. -1484 -1485 Parameters -1486 ---------- -1487 obs : list or numpy.ndarray -1488 List or one dimensional array of Obs -1489 visualize : bool -1490 If True plots the corresponding normalized correlation matrix (default False). -1491 correlation : bool -1492 If True the correlation matrix instead of the error covariance matrix is returned (default False). -1493 smooth : None or int -1494 If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue -1495 smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the -1496 largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely -1497 small ones. -1498 -1499 Notes -1500 ----- -1501 The error covariance is defined such that it agrees with the squared standard error for two identical observables -1502 $$\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$$ -1503 in the absence of autocorrelation. -1504 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 -1505 $$\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. -1506 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. -1507 $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$ -1508 This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors). -1509 ''' -1510 -1511 length = len(obs) -1512 -1513 max_samples = np.max([o.N for o in obs]) -1514 if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]: -1515 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) -1516 -1517 cov = np.zeros((length, length)) -1518 for i in range(length): -1519 for j in range(i, length): -1520 cov[i, j] = _covariance_element(obs[i], obs[j]) -1521 cov = cov + cov.T - np.diag(np.diag(cov)) -1522 -1523 corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov))) +1480 The gamma method has to be applied first to all observables. +1481 +1482 Parameters +1483 ---------- +1484 obs : list or numpy.ndarray +1485 List or one dimensional array of Obs +1486 visualize : bool +1487 If True plots the corresponding normalized correlation matrix (default False). +1488 correlation : bool +1489 If True the correlation matrix instead of the error covariance matrix is returned (default False). +1490 smooth : None or int +1491 If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue +1492 smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the +1493 largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely +1494 small ones. +1495 +1496 Notes +1497 ----- +1498 The error covariance is defined such that it agrees with the squared standard error for two identical observables +1499 $$\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$$ +1500 in the absence of autocorrelation. +1501 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 +1502 $$\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. +1503 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. +1504 $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$ +1505 This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors). +1506 ''' +1507 +1508 length = len(obs) +1509 +1510 max_samples = np.max([o.N for o in obs]) +1511 if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]: +1512 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) +1513 +1514 cov = np.zeros((length, length)) +1515 for i in range(length): +1516 for j in range(i, length): +1517 cov[i, j] = _covariance_element(obs[i], obs[j]) +1518 cov = cov + cov.T - np.diag(np.diag(cov)) +1519 +1520 corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov))) +1521 +1522 if isinstance(smooth, int): +1523 corr = _smooth_eigenvalues(corr, smooth) 1524 -1525 if isinstance(smooth, int): -1526 corr = _smooth_eigenvalues(corr, smooth) -1527 -1528 if visualize: -1529 plt.matshow(corr, vmin=-1, vmax=1) -1530 plt.set_cmap('RdBu') -1531 plt.colorbar() -1532 plt.draw() +1525 if visualize: +1526 plt.matshow(corr, vmin=-1, vmax=1) +1527 plt.set_cmap('RdBu') +1528 plt.colorbar() +1529 plt.draw() +1530 +1531 if correlation is True: +1532 return corr 1533 -1534 if correlation is True: -1535 return corr +1534 errors = [o.dvalue for o in obs] +1535 cov = np.diag(errors) @ corr @ np.diag(errors) 1536 -1537 errors = [o.dvalue for o in obs] -1538 cov = np.diag(errors) @ corr @ np.diag(errors) -1539 -1540 eigenvalues = np.linalg.eigh(cov)[0] -1541 if not np.all(eigenvalues >= 0): -1542 warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning) +1537 eigenvalues = np.linalg.eigh(cov)[0] +1538 if not np.all(eigenvalues >= 0): +1539 warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning) +1540 +1541 return cov +1542 1543 -1544 return cov -1545 -1546 -1547def invert_corr_cov_cholesky(corr, inverrdiag): -1548 """Constructs a lower triangular matrix `chol` via the Cholesky decomposition of the correlation matrix `corr` -1549 and then returns the inverse covariance matrix `chol_inv` as a lower triangular matrix by solving `chol * x = inverrdiag`. -1550 -1551 Parameters -1552 ---------- -1553 corr : np.ndarray -1554 correlation matrix -1555 inverrdiag : np.ndarray -1556 diagonal matrix, the entries are the inverse errors of the data points considered -1557 """ -1558 -1559 condn = np.linalg.cond(corr) -1560 if condn > 0.1 / np.finfo(float).eps: -1561 raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})") -1562 if condn > 1e13: -1563 warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning) -1564 chol = np.linalg.cholesky(corr) -1565 chol_inv = scipy.linalg.solve_triangular(chol, inverrdiag, lower=True) +1544def invert_corr_cov_cholesky(corr, inverrdiag): +1545 """Constructs a lower triangular matrix `chol` via the Cholesky decomposition of the correlation matrix `corr` +1546 and then returns the inverse covariance matrix `chol_inv` as a lower triangular matrix by solving `chol * x = inverrdiag`. +1547 +1548 Parameters +1549 ---------- +1550 corr : np.ndarray +1551 correlation matrix +1552 inverrdiag : np.ndarray +1553 diagonal matrix, the entries are the inverse errors of the data points considered +1554 """ +1555 +1556 condn = np.linalg.cond(corr) +1557 if condn > 0.1 / np.finfo(float).eps: +1558 raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})") +1559 if condn > 1e13: +1560 warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning) +1561 chol = np.linalg.cholesky(corr) +1562 chol_inv = scipy.linalg.solve_triangular(chol, inverrdiag, lower=True) +1563 +1564 return chol_inv +1565 1566 -1567 return chol_inv -1568 +1567def sort_corr(corr, kl, yd): +1568 """ Reorders a correlation matrix to match the alphabetical order of its underlying y data. 1569 -1570def sort_corr(corr, kl, yd): -1571 """ Reorders a correlation matrix to match the alphabetical order of its underlying y data. -1572 -1573 The ordering of the input correlation matrix `corr` is given by the list of keys `kl`. -1574 The input dictionary `yd` (with the same keys `kl`) must contain the corresponding y data -1575 that the correlation matrix is based on. -1576 This function sorts the list of keys `kl` alphabetically and sorts the matrix `corr` -1577 according to this alphabetical order such that the sorted matrix `corr_sorted` corresponds -1578 to the y data `yd` when arranged in an alphabetical order by its keys. -1579 -1580 Parameters -1581 ---------- -1582 corr : np.ndarray -1583 A square correlation matrix constructed using the order of the y data specified by `kl`. -1584 The dimensions of `corr` should match the total number of y data points in `yd` combined. -1585 kl : list of str -1586 A list of keys that denotes the order in which the y data from `yd` was used to build the -1587 input correlation matrix `corr`. -1588 yd : dict of list -1589 A dictionary where each key corresponds to a unique identifier, and its value is a list of -1590 y data points. The total number of y data points across all keys must match the dimensions -1591 of `corr`. The lists in the dictionary can be lists of Obs. -1592 -1593 Returns -1594 ------- -1595 np.ndarray -1596 A new, sorted correlation matrix that corresponds to the y data from `yd` when arranged alphabetically by its keys. -1597 -1598 Example -1599 ------- -1600 >>> import numpy as np -1601 >>> import pyerrors as pe -1602 >>> corr = np.array([[1, 0.2, 0.3], [0.2, 1, 0.4], [0.3, 0.4, 1]]) -1603 >>> kl = ['b', 'a'] -1604 >>> yd = {'a': [1, 2], 'b': [3]} -1605 >>> sorted_corr = pe.obs.sort_corr(corr, kl, yd) -1606 >>> print(sorted_corr) -1607 array([[1. , 0.3, 0.4], -1608 [0.3, 1. , 0.2], -1609 [0.4, 0.2, 1. ]]) +1570 The ordering of the input correlation matrix `corr` is given by the list of keys `kl`. +1571 The input dictionary `yd` (with the same keys `kl`) must contain the corresponding y data +1572 that the correlation matrix is based on. +1573 This function sorts the list of keys `kl` alphabetically and sorts the matrix `corr` +1574 according to this alphabetical order such that the sorted matrix `corr_sorted` corresponds +1575 to the y data `yd` when arranged in an alphabetical order by its keys. +1576 +1577 Parameters +1578 ---------- +1579 corr : np.ndarray +1580 A square correlation matrix constructed using the order of the y data specified by `kl`. +1581 The dimensions of `corr` should match the total number of y data points in `yd` combined. +1582 kl : list of str +1583 A list of keys that denotes the order in which the y data from `yd` was used to build the +1584 input correlation matrix `corr`. +1585 yd : dict of list +1586 A dictionary where each key corresponds to a unique identifier, and its value is a list of +1587 y data points. The total number of y data points across all keys must match the dimensions +1588 of `corr`. The lists in the dictionary can be lists of Obs. +1589 +1590 Returns +1591 ------- +1592 np.ndarray +1593 A new, sorted correlation matrix that corresponds to the y data from `yd` when arranged alphabetically by its keys. +1594 +1595 Example +1596 ------- +1597 >>> import numpy as np +1598 >>> import pyerrors as pe +1599 >>> corr = np.array([[1, 0.2, 0.3], [0.2, 1, 0.4], [0.3, 0.4, 1]]) +1600 >>> kl = ['b', 'a'] +1601 >>> yd = {'a': [1, 2], 'b': [3]} +1602 >>> sorted_corr = pe.obs.sort_corr(corr, kl, yd) +1603 >>> print(sorted_corr) +1604 array([[1. , 0.3, 0.4], +1605 [0.3, 1. , 0.2], +1606 [0.4, 0.2, 1. ]]) +1607 +1608 """ +1609 kl_sorted = sorted(kl) 1610 -1611 """ -1612 kl_sorted = sorted(kl) -1613 -1614 posd = {} -1615 ofs = 0 -1616 for ki, k in enumerate(kl): -1617 posd[k] = [i + ofs for i in range(len(yd[k]))] -1618 ofs += len(posd[k]) -1619 -1620 mapping = [] -1621 for k in kl_sorted: -1622 for i in range(len(yd[k])): -1623 mapping.append(posd[k][i]) -1624 -1625 corr_sorted = np.zeros_like(corr) -1626 for i in range(corr.shape[0]): -1627 for j in range(corr.shape[0]): -1628 corr_sorted[i][j] = corr[mapping[i]][mapping[j]] +1611 posd = {} +1612 ofs = 0 +1613 for ki, k in enumerate(kl): +1614 posd[k] = [i + ofs for i in range(len(yd[k]))] +1615 ofs += len(posd[k]) +1616 +1617 mapping = [] +1618 for k in kl_sorted: +1619 for i in range(len(yd[k])): +1620 mapping.append(posd[k][i]) +1621 +1622 corr_sorted = np.zeros_like(corr) +1623 for i in range(corr.shape[0]): +1624 for j in range(corr.shape[0]): +1625 corr_sorted[i][j] = corr[mapping[i]][mapping[j]] +1626 +1627 return corr_sorted +1628 1629 -1630 return corr_sorted -1631 +1630def _smooth_eigenvalues(corr, E): +1631 """Eigenvalue smoothing as described in hep-lat/9412087 1632 -1633def _smooth_eigenvalues(corr, E): -1634 """Eigenvalue smoothing as described in hep-lat/9412087 -1635 -1636 corr : np.ndarray -1637 correlation matrix -1638 E : integer -1639 Number of eigenvalues to be left substantially unchanged -1640 """ -1641 if not (2 < E < corr.shape[0] - 1): -1642 raise Exception(f"'E' has to be between 2 and the dimension of the correlation matrix minus 1 ({corr.shape[0] - 1}).") -1643 vals, vec = np.linalg.eigh(corr) -1644 lambda_min = np.mean(vals[:-E]) -1645 vals[vals < lambda_min] = lambda_min -1646 vals /= np.mean(vals) -1647 return vec @ np.diag(vals) @ vec.T -1648 +1633 corr : np.ndarray +1634 correlation matrix +1635 E : integer +1636 Number of eigenvalues to be left substantially unchanged +1637 """ +1638 if not (2 < E < corr.shape[0] - 1): +1639 raise Exception(f"'E' has to be between 2 and the dimension of the correlation matrix minus 1 ({corr.shape[0] - 1}).") +1640 vals, vec = np.linalg.eigh(corr) +1641 lambda_min = np.mean(vals[:-E]) +1642 vals[vals < lambda_min] = lambda_min +1643 vals /= np.mean(vals) +1644 return vec @ np.diag(vals) @ vec.T +1645 +1646 +1647def _covariance_element(obs1, obs2): +1648 """Estimates the covariance of two Obs objects, neglecting autocorrelations.""" 1649 -1650def _covariance_element(obs1, obs2): -1651 """Estimates the covariance of two Obs objects, neglecting autocorrelations.""" -1652 -1653 def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx): -1654 deltas1 = _reduce_deltas(deltas1, idx1, new_idx) -1655 deltas2 = _reduce_deltas(deltas2, idx2, new_idx) -1656 return np.sum(deltas1 * deltas2) +1650 def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx): +1651 deltas1 = _reduce_deltas(deltas1, idx1, new_idx) +1652 deltas2 = _reduce_deltas(deltas2, idx2, new_idx) +1653 return np.sum(deltas1 * deltas2) +1654 +1655 if set(obs1.names).isdisjoint(set(obs2.names)): +1656 return 0.0 1657 -1658 if set(obs1.names).isdisjoint(set(obs2.names)): -1659 return 0.0 +1658 if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'): +1659 raise Exception('The gamma method has to be applied to both Obs first.') 1660 -1661 if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'): -1662 raise Exception('The gamma method has to be applied to both Obs first.') -1663 -1664 dvalue = 0.0 -1665 -1666 for e_name in obs1.mc_names: +1661 dvalue = 0.0 +1662 +1663 for e_name in obs1.mc_names: +1664 +1665 if e_name not in obs2.mc_names: +1666 continue 1667 -1668 if e_name not in obs2.mc_names: -1669 continue -1670 -1671 idl_d = {} -1672 for r_name in obs1.e_content[e_name]: -1673 if r_name not in obs2.e_content[e_name]: -1674 continue -1675 idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]]) -1676 -1677 gamma = 0.0 -1678 -1679 for r_name in obs1.e_content[e_name]: -1680 if r_name not in obs2.e_content[e_name]: -1681 continue -1682 if len(idl_d[r_name]) == 0: -1683 continue -1684 gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name]) +1668 idl_d = {} +1669 for r_name in obs1.e_content[e_name]: +1670 if r_name not in obs2.e_content[e_name]: +1671 continue +1672 idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]]) +1673 +1674 gamma = 0.0 +1675 +1676 for r_name in obs1.e_content[e_name]: +1677 if r_name not in obs2.e_content[e_name]: +1678 continue +1679 if len(idl_d[r_name]) == 0: +1680 continue +1681 gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name]) +1682 +1683 if gamma == 0.0: +1684 continue 1685 -1686 if gamma == 0.0: -1687 continue -1688 -1689 gamma_div = 0.0 -1690 for r_name in obs1.e_content[e_name]: -1691 if r_name not in obs2.e_content[e_name]: -1692 continue -1693 if len(idl_d[r_name]) == 0: -1694 continue -1695 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])) -1696 gamma /= gamma_div -1697 -1698 dvalue += gamma -1699 -1700 for e_name in obs1.cov_names: +1686 gamma_div = 0.0 +1687 for r_name in obs1.e_content[e_name]: +1688 if r_name not in obs2.e_content[e_name]: +1689 continue +1690 if len(idl_d[r_name]) == 0: +1691 continue +1692 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])) +1693 gamma /= gamma_div +1694 +1695 dvalue += gamma +1696 +1697 for e_name in obs1.cov_names: +1698 +1699 if e_name not in obs2.cov_names: +1700 continue 1701 -1702 if e_name not in obs2.cov_names: -1703 continue -1704 -1705 dvalue += np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad)).item() +1702 dvalue += np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad)).item() +1703 +1704 return dvalue +1705 1706 -1707 return dvalue -1708 +1707def import_jackknife(jacks, name, idl=None): +1708 """Imports jackknife samples and returns an Obs 1709 -1710def import_jackknife(jacks, name, idl=None): -1711 """Imports jackknife samples and returns an Obs -1712 -1713 Parameters -1714 ---------- -1715 jacks : numpy.ndarray -1716 numpy array containing the mean value as zeroth entry and -1717 the N jackknife samples as first to Nth entry. -1718 name : str -1719 name of the ensemble the samples are defined on. -1720 """ -1721 length = len(jacks) - 1 -1722 prj = (np.ones((length, length)) - (length - 1) * np.identity(length)) -1723 samples = jacks[1:] @ prj -1724 mean = np.mean(samples) -1725 new_obs = Obs([samples - mean], [name], idl=idl, means=[mean]) -1726 new_obs._value = jacks[0] -1727 return new_obs -1728 +1710 Parameters +1711 ---------- +1712 jacks : numpy.ndarray +1713 numpy array containing the mean value as zeroth entry and +1714 the N jackknife samples as first to Nth entry. +1715 name : str +1716 name of the ensemble the samples are defined on. +1717 """ +1718 length = len(jacks) - 1 +1719 prj = (np.ones((length, length)) - (length - 1) * np.identity(length)) +1720 samples = jacks[1:] @ prj +1721 mean = np.mean(samples) +1722 new_obs = Obs([samples - mean], [name], idl=idl, means=[mean]) +1723 new_obs._value = jacks[0] +1724 return new_obs +1725 +1726 +1727def import_bootstrap(boots, name, random_numbers): +1728 """Imports bootstrap samples and returns an Obs 1729 -1730def import_bootstrap(boots, name, random_numbers): -1731 """Imports bootstrap samples and returns an Obs -1732 -1733 Parameters -1734 ---------- -1735 boots : numpy.ndarray -1736 numpy array containing the mean value as zeroth entry and -1737 the N bootstrap samples as first to Nth entry. -1738 name : str -1739 name of the ensemble the samples are defined on. -1740 random_numbers : np.ndarray -1741 Array of shape (samples, length) containing the random numbers to generate the bootstrap samples, -1742 where samples is the number of bootstrap samples and length is the length of the original Monte Carlo -1743 chain to be reconstructed. -1744 """ -1745 samples, length = random_numbers.shape -1746 if samples != len(boots) - 1: -1747 raise ValueError("Random numbers do not have the correct shape.") +1730 Parameters +1731 ---------- +1732 boots : numpy.ndarray +1733 numpy array containing the mean value as zeroth entry and +1734 the N bootstrap samples as first to Nth entry. +1735 name : str +1736 name of the ensemble the samples are defined on. +1737 random_numbers : np.ndarray +1738 Array of shape (samples, length) containing the random numbers to generate the bootstrap samples, +1739 where samples is the number of bootstrap samples and length is the length of the original Monte Carlo +1740 chain to be reconstructed. +1741 """ +1742 samples, length = random_numbers.shape +1743 if samples != len(boots) - 1: +1744 raise ValueError("Random numbers do not have the correct shape.") +1745 +1746 if samples < length: +1747 raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.") 1748 -1749 if samples < length: -1750 raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.") -1751 -1752 proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length -1753 -1754 samples = scipy.linalg.lstsq(proj, boots[1:])[0] -1755 ret = Obs([samples], [name]) -1756 ret._value = boots[0] -1757 return ret -1758 +1749 proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length +1750 +1751 samples = scipy.linalg.lstsq(proj, boots[1:])[0] +1752 ret = Obs([samples], [name]) +1753 ret._value = boots[0] +1754 return ret +1755 +1756 +1757def merge_obs(list_of_obs): +1758 """Combine all observables in list_of_obs into one new observable 1759 -1760def merge_obs(list_of_obs): -1761 """Combine all observables in list_of_obs into one new observable -1762 -1763 Parameters -1764 ---------- -1765 list_of_obs : list -1766 list of the Obs object to be combined -1767 -1768 Notes -1769 ----- -1770 It is not possible to combine obs which are based on the same replicum -1771 """ -1772 replist = [item for obs in list_of_obs for item in obs.names] -1773 if (len(replist) == len(set(replist))) is False: -1774 raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist))) -1775 if any([len(o.cov_names) for o in list_of_obs]): -1776 raise Exception('Not possible to merge data that contains covobs!') -1777 new_dict = {} -1778 idl_dict = {} -1779 for o in list_of_obs: -1780 new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0) -1781 for key in set(o.deltas) | set(o.r_values)}) -1782 idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)}) -1783 -1784 names = sorted(new_dict.keys()) -1785 o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names]) -1786 o.reweighted = np.max([oi.reweighted for oi in list_of_obs]) -1787 return o -1788 +1760 Parameters +1761 ---------- +1762 list_of_obs : list +1763 list of the Obs object to be combined +1764 +1765 Notes +1766 ----- +1767 It is not possible to combine obs which are based on the same replicum +1768 """ +1769 replist = [item for obs in list_of_obs for item in obs.names] +1770 if (len(replist) == len(set(replist))) is False: +1771 raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist))) +1772 if any([len(o.cov_names) for o in list_of_obs]): +1773 raise Exception('Not possible to merge data that contains covobs!') +1774 new_dict = {} +1775 idl_dict = {} +1776 for o in list_of_obs: +1777 new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0) +1778 for key in set(o.deltas) | set(o.r_values)}) +1779 idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)}) +1780 +1781 names = sorted(new_dict.keys()) +1782 o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names]) +1783 o.reweighted = np.max([oi.reweighted for oi in list_of_obs]) +1784 return o +1785 +1786 +1787def cov_Obs(means, cov, name, grad=None): +1788 """Create an Obs based on mean(s) and a covariance matrix 1789 -1790def cov_Obs(means, cov, name, grad=None): -1791 """Create an Obs based on mean(s) and a covariance matrix -1792 -1793 Parameters -1794 ---------- -1795 mean : list of floats or float -1796 N mean value(s) of the new Obs -1797 cov : list or array -1798 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance -1799 name : str -1800 identifier for the covariance matrix -1801 grad : list or array -1802 Gradient of the Covobs wrt. the means belonging to cov. -1803 """ +1790 Parameters +1791 ---------- +1792 mean : list of floats or float +1793 N mean value(s) of the new Obs +1794 cov : list or array +1795 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance +1796 name : str +1797 identifier for the covariance matrix +1798 grad : list or array +1799 Gradient of the Covobs wrt. the means belonging to cov. +1800 """ +1801 +1802 def covobs_to_obs(co): +1803 """Make an Obs out of a Covobs 1804 -1805 def covobs_to_obs(co): -1806 """Make an Obs out of a Covobs -1807 -1808 Parameters -1809 ---------- -1810 co : Covobs -1811 Covobs to be embedded into the Obs -1812 """ -1813 o = Obs([], [], means=[]) -1814 o._value = co.value -1815 o.names.append(co.name) -1816 o._covobs[co.name] = co -1817 o._dvalue = np.sqrt(co.errsq()) -1818 return o -1819 -1820 ol = [] -1821 if isinstance(means, (float, int)): -1822 means = [means] -1823 -1824 for i in range(len(means)): -1825 ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad))) -1826 if ol[0].covobs[name].N != len(means): -1827 raise Exception('You have to provide %d mean values!' % (ol[0].N)) -1828 if len(ol) == 1: -1829 return ol[0] -1830 return ol -1831 -1832 -1833def _determine_gap(o, e_content, e_name): -1834 gaps = [] -1835 for r_name in e_content[e_name]: -1836 if isinstance(o.idl[r_name], range): -1837 gaps.append(o.idl[r_name].step) -1838 else: -1839 gaps.append(np.min(np.diff(o.idl[r_name]))) -1840 -1841 gap = min(gaps) -1842 if not np.all([gi % gap == 0 for gi in gaps]): -1843 raise Exception(f"Replica for ensemble {e_name} do not have a common spacing.", gaps) +1805 Parameters +1806 ---------- +1807 co : Covobs +1808 Covobs to be embedded into the Obs +1809 """ +1810 o = Obs([], [], means=[]) +1811 o._value = co.value +1812 o.names.append(co.name) +1813 o._covobs[co.name] = co +1814 o._dvalue = np.sqrt(co.errsq()) +1815 return o +1816 +1817 ol = [] +1818 if isinstance(means, (float, int)): +1819 means = [means] +1820 +1821 for i in range(len(means)): +1822 ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad))) +1823 if ol[0].covobs[name].N != len(means): +1824 raise Exception('You have to provide %d mean values!' % (ol[0].N)) +1825 if len(ol) == 1: +1826 return ol[0] +1827 return ol +1828 +1829 +1830def _determine_gap(o, e_content, e_name): +1831 gaps = [] +1832 for r_name in e_content[e_name]: +1833 if isinstance(o.idl[r_name], range): +1834 gaps.append(o.idl[r_name].step) +1835 else: +1836 gaps.append(np.min(np.diff(o.idl[r_name]))) +1837 +1838 gap = min(gaps) +1839 if not np.all([gi % gap == 0 for gi in gaps]): +1840 raise Exception(f"Replica for ensemble {e_name} do not have a common spacing.", gaps) +1841 +1842 return gap +1843 1844 -1845 return gap -1846 -1847 -1848def _check_lists_equal(idl): -1849 ''' -1850 Use groupby to efficiently check whether all elements of idl are identical. -1851 Returns True if all elements are equal, otherwise False. -1852 -1853 Parameters -1854 ---------- -1855 idl : list of lists, ranges or np.ndarrays -1856 ''' -1857 g = groupby([np.nditer(el) if isinstance(el, np.ndarray) else el for el in idl]) -1858 if next(g, True) and not next(g, False): -1859 return True -1860 return False +1845def _check_lists_equal(idl): +1846 ''' +1847 Use groupby to efficiently check whether all elements of idl are identical. +1848 Returns True if all elements are equal, otherwise False. +1849 +1850 Parameters +1851 ---------- +1852 idl : list of lists, ranges or np.ndarrays +1853 ''' +1854 g = groupby([np.nditer(el) if isinstance(el, np.ndarray) else el for el in idl]) +1855 if next(g, True) and not next(g, False): +1856 return True +1857 return False @@ -3050,64 +3047,61 @@ 857 858 def __pow__(self, y): 859 if isinstance(y, Obs): -860 return derived_observable(lambda x: x[0] ** x[1], [self, y]) +860 return derived_observable(lambda x, **kwargs: x[0] ** x[1], [self, y], man_grad=[y.value * self.value ** (y.value - 1), self.value ** y.value * np.log(self.value)]) 861 else: -862 return derived_observable(lambda x: x[0] ** y, [self]) +862 return derived_observable(lambda x, **kwargs: x[0] ** y, [self], man_grad=[y * self.value ** (y - 1)]) 863 864 def __rpow__(self, y): -865 if isinstance(y, Obs): -866 return derived_observable(lambda x: x[0] ** x[1], [y, self]) -867 else: -868 return derived_observable(lambda x: y ** x[0], [self]) +865 return derived_observable(lambda x, **kwargs: y ** x[0], [self], man_grad=[y ** self.value * np.log(y)]) +866 +867 def __abs__(self): +868 return derived_observable(lambda x: anp.abs(x[0]), [self]) 869 -870 def __abs__(self): -871 return derived_observable(lambda x: anp.abs(x[0]), [self]) -872 -873 # Overload numpy functions -874 def sqrt(self): -875 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)]) +870 # Overload numpy functions +871 def sqrt(self): +872 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)]) +873 +874 def log(self): +875 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value]) 876 -877 def log(self): -878 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value]) +877 def exp(self): +878 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)]) 879 -880 def exp(self): -881 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)]) +880 def sin(self): +881 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)]) 882 -883 def sin(self): -884 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)]) +883 def cos(self): +884 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)]) 885 -886 def cos(self): -887 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)]) +886 def tan(self): +887 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2]) 888 -889 def tan(self): -890 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2]) +889 def arcsin(self): +890 return derived_observable(lambda x: anp.arcsin(x[0]), [self]) 891 -892 def arcsin(self): -893 return derived_observable(lambda x: anp.arcsin(x[0]), [self]) +892 def arccos(self): +893 return derived_observable(lambda x: anp.arccos(x[0]), [self]) 894 -895 def arccos(self): -896 return derived_observable(lambda x: anp.arccos(x[0]), [self]) +895 def arctan(self): +896 return derived_observable(lambda x: anp.arctan(x[0]), [self]) 897 -898 def arctan(self): -899 return derived_observable(lambda x: anp.arctan(x[0]), [self]) +898 def sinh(self): +899 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)]) 900 -901 def sinh(self): -902 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)]) +901 def cosh(self): +902 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)]) 903 -904 def cosh(self): -905 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)]) +904 def tanh(self): +905 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2]) 906 -907 def tanh(self): -908 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2]) +907 def arcsinh(self): +908 return derived_observable(lambda x: anp.arcsinh(x[0]), [self]) 909 -910 def arcsinh(self): -911 return derived_observable(lambda x: anp.arcsinh(x[0]), [self]) +910 def arccosh(self): +911 return derived_observable(lambda x: anp.arccosh(x[0]), [self]) 912 -913 def arccosh(self): -914 return derived_observable(lambda x: anp.arccosh(x[0]), [self]) -915 -916 def arctanh(self): -917 return derived_observable(lambda x: anp.arctanh(x[0]), [self]) +913 def arctanh(self): +914 return derived_observable(lambda x: anp.arctanh(x[0]), [self]) @@ -4653,8 +4647,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N). - 
874    def sqrt(self):
-875        return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
+            
871    def sqrt(self):
+872        return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
 

 
 
@@ -4672,8 +4666,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
877    def log(self):
-878        return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
+            
874    def log(self):
+875        return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
 

 
 
@@ -4691,8 +4685,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
880    def exp(self):
-881        return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
+            
877    def exp(self):
+878        return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
 

 
 
@@ -4710,8 +4704,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
883    def sin(self):
-884        return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
+            
880    def sin(self):
+881        return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
 

 
 
@@ -4729,8 +4723,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
886    def cos(self):
-887        return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
+            
883    def cos(self):
+884        return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
 

 
 
@@ -4748,8 +4742,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
889    def tan(self):
-890        return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
+            
886    def tan(self):
+887        return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
 

 
 
@@ -4767,8 +4761,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
892    def arcsin(self):
-893        return derived_observable(lambda x: anp.arcsin(x[0]), [self])
+            
889    def arcsin(self):
+890        return derived_observable(lambda x: anp.arcsin(x[0]), [self])
 

 
 
@@ -4786,8 +4780,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
895    def arccos(self):
-896        return derived_observable(lambda x: anp.arccos(x[0]), [self])
+            
892    def arccos(self):
+893        return derived_observable(lambda x: anp.arccos(x[0]), [self])
 

 
 
@@ -4805,8 +4799,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
898    def arctan(self):
-899        return derived_observable(lambda x: anp.arctan(x[0]), [self])
+            
895    def arctan(self):
+896        return derived_observable(lambda x: anp.arctan(x[0]), [self])
 

 
 
@@ -4824,8 +4818,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
901    def sinh(self):
-902        return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
+            
898    def sinh(self):
+899        return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
 

 
 
@@ -4843,8 +4837,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
904    def cosh(self):
-905        return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
+            
901    def cosh(self):
+902        return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
 

 
 
@@ -4862,8 +4856,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
907    def tanh(self):
-908        return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
+            
904    def tanh(self):
+905        return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
 

 
 
@@ -4881,8 +4875,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
910    def arcsinh(self):
-911        return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
+            
907    def arcsinh(self):
+908        return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
 

 
 
@@ -4900,8 +4894,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
913    def arccosh(self):
-914        return derived_observable(lambda x: anp.arccosh(x[0]), [self])
+            
910    def arccosh(self):
+911        return derived_observable(lambda x: anp.arccosh(x[0]), [self])
 

 
 
@@ -4919,8 +4913,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
916    def arctanh(self):
-917        return derived_observable(lambda x: anp.arctanh(x[0]), [self])
+            
913    def arctanh(self):
+914        return derived_observable(lambda x: anp.arctanh(x[0]), [self])
 

 
 
@@ -5071,123 +5065,123 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
 920class CObs:
- 921    """Class for a complex valued observable."""
- 922    __slots__ = ['_real', '_imag', 'tag']
- 923
- 924    def __init__(self, real, imag=0.0):
- 925        self._real = real
- 926        self._imag = imag
- 927        self.tag = None
- 928
- 929    @property
- 930    def real(self):
- 931        return self._real
- 932
- 933    @property
- 934    def imag(self):
- 935        return self._imag
- 936
- 937    def gamma_method(self, **kwargs):
- 938        """Executes the gamma_method for the real and the imaginary part."""
- 939        if isinstance(self.real, Obs):
- 940            self.real.gamma_method(**kwargs)
- 941        if isinstance(self.imag, Obs):
- 942            self.imag.gamma_method(**kwargs)
- 943
- 944    def is_zero(self):
- 945        """Checks whether both real and imaginary part are zero within machine precision."""
- 946        return self.real == 0.0 and self.imag == 0.0
+            
 917class CObs:
+ 918    """Class for a complex valued observable."""
+ 919    __slots__ = ['_real', '_imag', 'tag']
+ 920
+ 921    def __init__(self, real, imag=0.0):
+ 922        self._real = real
+ 923        self._imag = imag
+ 924        self.tag = None
+ 925
+ 926    @property
+ 927    def real(self):
+ 928        return self._real
+ 929
+ 930    @property
+ 931    def imag(self):
+ 932        return self._imag
+ 933
+ 934    def gamma_method(self, **kwargs):
+ 935        """Executes the gamma_method for the real and the imaginary part."""
+ 936        if isinstance(self.real, Obs):
+ 937            self.real.gamma_method(**kwargs)
+ 938        if isinstance(self.imag, Obs):
+ 939            self.imag.gamma_method(**kwargs)
+ 940
+ 941    def is_zero(self):
+ 942        """Checks whether both real and imaginary part are zero within machine precision."""
+ 943        return self.real == 0.0 and self.imag == 0.0
+ 944
+ 945    def conjugate(self):
+ 946        return CObs(self.real, -self.imag)
  947
- 948    def conjugate(self):
- 949        return CObs(self.real, -self.imag)
- 950
- 951    def __add__(self, other):
- 952        if isinstance(other, np.ndarray):
- 953            return other + self
- 954        elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 955            return CObs(self.real + other.real,
- 956                        self.imag + other.imag)
- 957        else:
- 958            return CObs(self.real + other, self.imag)
+ 948    def __add__(self, other):
+ 949        if isinstance(other, np.ndarray):
+ 950            return other + self
+ 951        elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 952            return CObs(self.real + other.real,
+ 953                        self.imag + other.imag)
+ 954        else:
+ 955            return CObs(self.real + other, self.imag)
+ 956
+ 957    def __radd__(self, y):
+ 958        return self + y
  959
- 960    def __radd__(self, y):
- 961        return self + y
- 962
- 963    def __sub__(self, other):
- 964        if isinstance(other, np.ndarray):
- 965            return -1 * (other - self)
- 966        elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 967            return CObs(self.real - other.real, self.imag - other.imag)
- 968        else:
- 969            return CObs(self.real - other, self.imag)
+ 960    def __sub__(self, other):
+ 961        if isinstance(other, np.ndarray):
+ 962            return -1 * (other - self)
+ 963        elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 964            return CObs(self.real - other.real, self.imag - other.imag)
+ 965        else:
+ 966            return CObs(self.real - other, self.imag)
+ 967
+ 968    def __rsub__(self, other):
+ 969        return -1 * (self - other)
  970
- 971    def __rsub__(self, other):
- 972        return -1 * (self - other)
- 973
- 974    def __mul__(self, other):
- 975        if isinstance(other, np.ndarray):
- 976            return other * self
- 977        elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 978            if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
- 979                return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
+ 971    def __mul__(self, other):
+ 972        if isinstance(other, np.ndarray):
+ 973            return other * self
+ 974        elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 975            if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
+ 976                return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
+ 977                                               [self.real, other.real, self.imag, other.imag],
+ 978                                               man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
+ 979                            derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
  980                                               [self.real, other.real, self.imag, other.imag],
- 981                                               man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
- 982                            derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
- 983                                               [self.real, other.real, self.imag, other.imag],
- 984                                               man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
- 985            elif getattr(other, 'imag', 0) != 0:
- 986                return CObs(self.real * other.real - self.imag * other.imag,
- 987                            self.imag * other.real + self.real * other.imag)
- 988            else:
- 989                return CObs(self.real * other.real, self.imag * other.real)
- 990        else:
- 991            return CObs(self.real * other, self.imag * other)
+ 981                                               man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
+ 982            elif getattr(other, 'imag', 0) != 0:
+ 983                return CObs(self.real * other.real - self.imag * other.imag,
+ 984                            self.imag * other.real + self.real * other.imag)
+ 985            else:
+ 986                return CObs(self.real * other.real, self.imag * other.real)
+ 987        else:
+ 988            return CObs(self.real * other, self.imag * other)
+ 989
+ 990    def __rmul__(self, other):
+ 991        return self * other
  992
- 993    def __rmul__(self, other):
- 994        return self * other
- 995
- 996    def __truediv__(self, other):
- 997        if isinstance(other, np.ndarray):
- 998            return 1 / (other / self)
- 999        elif hasattr(other, 'real') and hasattr(other, 'imag'):
-1000            r = other.real ** 2 + other.imag ** 2
-1001            return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
-1002        else:
-1003            return CObs(self.real / other, self.imag / other)
-1004
-1005    def __rtruediv__(self, other):
-1006        r = self.real ** 2 + self.imag ** 2
-1007        if hasattr(other, 'real') and hasattr(other, 'imag'):
-1008            return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
-1009        else:
-1010            return CObs(self.real * other / r, -self.imag * other / r)
+ 993    def __truediv__(self, other):
+ 994        if isinstance(other, np.ndarray):
+ 995            return 1 / (other / self)
+ 996        elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 997            r = other.real ** 2 + other.imag ** 2
+ 998            return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
+ 999        else:
+1000            return CObs(self.real / other, self.imag / other)
+1001
+1002    def __rtruediv__(self, other):
+1003        r = self.real ** 2 + self.imag ** 2
+1004        if hasattr(other, 'real') and hasattr(other, 'imag'):
+1005            return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
+1006        else:
+1007            return CObs(self.real * other / r, -self.imag * other / r)
+1008
+1009    def __abs__(self):
+1010        return np.sqrt(self.real**2 + self.imag**2)
 1011
-1012    def __abs__(self):
-1013        return np.sqrt(self.real**2 + self.imag**2)
+1012    def __pos__(self):
+1013        return self
 1014
-1015    def __pos__(self):
-1016        return self
+1015    def __neg__(self):
+1016        return -1 * self
 1017
-1018    def __neg__(self):
-1019        return -1 * self
+1018    def __eq__(self, other):
+1019        return self.real == other.real and self.imag == other.imag
 1020
-1021    def __eq__(self, other):
-1022        return self.real == other.real and self.imag == other.imag
+1021    def __str__(self):
+1022        return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
 1023
-1024    def __str__(self):
-1025        return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
+1024    def __repr__(self):
+1025        return 'CObs[' + str(self) + ']'
 1026
-1027    def __repr__(self):
-1028        return 'CObs[' + str(self) + ']'
-1029
-1030    def __format__(self, format_type):
-1031        if format_type == "":
-1032            significance = 2
-1033            format_type = "2"
-1034        else:
-1035            significance = int(float(format_type.replace("+", "").replace("-", "")))
-1036        return f"({self.real:{format_type}}{self.imag:+{significance}}j)"
+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)"
 

 
 
@@ -5205,10 +5199,10 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
924    def __init__(self, real, imag=0.0):
-925        self._real = real
-926        self._imag = imag
-927        self.tag = None
+            
921    def __init__(self, real, imag=0.0):
+922        self._real = real
+923        self._imag = imag
+924        self.tag = None
 

 
 
@@ -5235,9 +5229,9 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
929    @property
-930    def real(self):
-931        return self._real
+            
926    @property
+927    def real(self):
+928        return self._real
 

 
 
@@ -5253,9 +5247,9 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
933    @property
-934    def imag(self):
-935        return self._imag
+            
930    @property
+931    def imag(self):
+932        return self._imag
 

 
 
@@ -5273,12 +5267,12 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
937    def gamma_method(self, **kwargs):
-938        """Executes the gamma_method for the real and the imaginary part."""
-939        if isinstance(self.real, Obs):
-940            self.real.gamma_method(**kwargs)
-941        if isinstance(self.imag, Obs):
-942            self.imag.gamma_method(**kwargs)
+            
934    def gamma_method(self, **kwargs):
+935        """Executes the gamma_method for the real and the imaginary part."""
+936        if isinstance(self.real, Obs):
+937            self.real.gamma_method(**kwargs)
+938        if isinstance(self.imag, Obs):
+939            self.imag.gamma_method(**kwargs)
 

 
 
@@ -5298,9 +5292,9 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
944    def is_zero(self):
-945        """Checks whether both real and imaginary part are zero within machine precision."""
-946        return self.real == 0.0 and self.imag == 0.0
+            
941    def is_zero(self):
+942        """Checks whether both real and imaginary part are zero within machine precision."""
+943        return self.real == 0.0 and self.imag == 0.0
 

 
 
@@ -5320,8 +5314,8 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
948    def conjugate(self):
-949        return CObs(self.real, -self.imag)
+            
945    def conjugate(self):
+946        return CObs(self.real, -self.imag)
 

 
 
@@ -5340,12 +5334,12 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
1039def gamma_method(x, **kwargs):
-1040    """Vectorized version of the gamma_method applicable to lists or arrays of Obs.
-1041
-1042    See docstring of pe.Obs.gamma_method for details.
-1043    """
-1044    return np.vectorize(lambda o: o.gm(**kwargs))(x)
+            
1036def gamma_method(x, **kwargs):
+1037    """Vectorized version of the gamma_method applicable to lists or arrays of Obs.
+1038
+1039    See docstring of pe.Obs.gamma_method for details.
+1040    """
+1041    return np.vectorize(lambda o: o.gm(**kwargs))(x)
 

 
 
@@ -5367,12 +5361,12 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
1039def gamma_method(x, **kwargs):
-1040    """Vectorized version of the gamma_method applicable to lists or arrays of Obs.
-1041
-1042    See docstring of pe.Obs.gamma_method for details.
-1043    """
-1044    return np.vectorize(lambda o: o.gm(**kwargs))(x)
+            
1036def gamma_method(x, **kwargs):
+1037    """Vectorized version of the gamma_method applicable to lists or arrays of Obs.
+1038
+1039    See docstring of pe.Obs.gamma_method for details.
+1040    """
+1041    return np.vectorize(lambda o: o.gm(**kwargs))(x)
 

 
 
@@ -5394,194 +5388,194 @@ should agree with samples from a full bootstrap analysis up to O(1/N).
- 
1174def derived_observable(func, data, array_mode=False, **kwargs):
-1175    """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
-1176
-1177    Parameters
-1178    ----------
-1179    func : object
-1180        arbitrary function of the form func(data, **kwargs). For the
-1181        automatic differentiation to work, all numpy functions have to have
-1182        the autograd wrapper (use 'import autograd.numpy as anp').
-1183    data : list
-1184        list of Obs, e.g. [obs1, obs2, obs3].
-1185    num_grad : bool
-1186        if True, numerical derivatives are used instead of autograd
-1187        (default False). To control the numerical differentiation the
-1188        kwargs of numdifftools.step_generators.MaxStepGenerator
-1189        can be used.
-1190    man_grad : list
-1191        manually supply a list or an array which contains the jacobian
-1192        of func. Use cautiously, supplying the wrong derivative will
-1193        not be intercepted.
-1194
-1195    Notes
-1196    -----
-1197    For simple mathematical operations it can be practical to use anonymous
-1198    functions. For the ratio of two observables one can e.g. use
+            
1171def derived_observable(func, data, array_mode=False, **kwargs):
+1172    """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
+1173
+1174    Parameters
+1175    ----------
+1176    func : object
+1177        arbitrary function of the form func(data, **kwargs). For the
+1178        automatic differentiation to work, all numpy functions have to have
+1179        the autograd wrapper (use 'import autograd.numpy as anp').
+1180    data : list
+1181        list of Obs, e.g. [obs1, obs2, obs3].
+1182    num_grad : bool
+1183        if True, numerical derivatives are used instead of autograd
+1184        (default False). To control the numerical differentiation the
+1185        kwargs of numdifftools.step_generators.MaxStepGenerator
+1186        can be used.
+1187    man_grad : list
+1188        manually supply a list or an array which contains the jacobian
+1189        of func. Use cautiously, supplying the wrong derivative will
+1190        not be intercepted.
+1191
+1192    Notes
+1193    -----
+1194    For simple mathematical operations it can be practical to use anonymous
+1195    functions. For the ratio of two observables one can e.g. use
+1196
+1197    new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
+1198    """
 1199
-1200    new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
-1201    """
+1200    data = np.asarray(data)
+1201    raveled_data = data.ravel()
 1202
-1203    data = np.asarray(data)
-1204    raveled_data = data.ravel()
-1205
-1206    # Workaround for matrix operations containing non Obs data
-1207    if not all(isinstance(x, Obs) for x in raveled_data):
-1208        for i in range(len(raveled_data)):
-1209            if isinstance(raveled_data[i], (int, float)):
-1210                raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
-1211
-1212    allcov = {}
-1213    for o in raveled_data:
-1214        for name in o.cov_names:
-1215            if name in allcov:
-1216                if not np.allclose(allcov[name], o.covobs[name].cov):
-1217                    raise Exception('Inconsistent covariance matrices for %s!' % (name))
-1218            else:
-1219                allcov[name] = o.covobs[name].cov
-1220
-1221    n_obs = len(raveled_data)
-1222    new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
-1223    new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
-1224    new_sample_names = sorted(set(new_names) - set(new_cov_names))
-1225
-1226    reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
-1227
-1228    if data.ndim == 1:
-1229        values = np.array([o.value for o in data])
-1230    else:
-1231        values = np.vectorize(lambda x: x.value)(data)
-1232
-1233    new_values = func(values, **kwargs)
-1234
-1235    multi = int(isinstance(new_values, np.ndarray))
-1236
-1237    new_r_values = {}
-1238    new_idl_d = {}
-1239    for name in new_sample_names:
-1240        idl = []
-1241        tmp_values = np.zeros(n_obs)
-1242        for i, item in enumerate(raveled_data):
-1243            tmp_values[i] = item.r_values.get(name, item.value)
-1244            tmp_idl = item.idl.get(name)
-1245            if tmp_idl is not None:
-1246                idl.append(tmp_idl)
-1247        if multi > 0:
-1248            tmp_values = np.array(tmp_values).reshape(data.shape)
-1249        new_r_values[name] = func(tmp_values, **kwargs)
-1250        new_idl_d[name] = _merge_idx(idl)
-1251
-1252    def _compute_scalefactor_missing_rep(obs):
-1253        """
-1254        Computes the scale factor that is to be multiplied with the deltas
-1255        in the case where Obs with different subsets of replica are merged.
-1256        Returns a dictionary with the scale factor for each Monte Carlo name.
-1257
-1258        Parameters
-1259        ----------
-1260        obs : Obs
-1261            The observable corresponding to the deltas that are to be scaled
-1262        """
-1263        scalef_d = {}
-1264        for mc_name in obs.mc_names:
-1265            mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')]
-1266            new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')]
-1267            if len(mc_idl_d) > 0 and len(mc_idl_d) < len(new_mc_idl_d):
-1268                scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d])
-1269        return scalef_d
-1270
-1271    if 'man_grad' in kwargs:
-1272        deriv = np.asarray(kwargs.get('man_grad'))
-1273        if new_values.shape + data.shape != deriv.shape:
-1274            raise Exception('Manual derivative does not have correct shape.')
-1275    elif kwargs.get('num_grad') is True:
-1276        if multi > 0:
-1277            raise Exception('Multi mode currently not supported for numerical derivative')
-1278        options = {
-1279            'base_step': 0.1,
-1280            'step_ratio': 2.5}
-1281        for key in options.keys():
-1282            kwarg = kwargs.get(key)
-1283            if kwarg is not None:
-1284                options[key] = kwarg
-1285        tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
-1286        if tmp_df.size == 1:
-1287            deriv = np.array([tmp_df.real])
-1288        else:
-1289            deriv = tmp_df.real
-1290    else:
-1291        deriv = jacobian(func)(values, **kwargs)
-1292
-1293    final_result = np.zeros(new_values.shape, dtype=object)
-1294
-1295    if array_mode is True:
-1296
-1297        class _Zero_grad():
-1298            def __init__(self, N):
-1299                self.grad = np.zeros((N, 1))
-1300
-1301        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]))
-1302        d_extracted = {}
-1303        g_extracted = {}
-1304        for name in new_sample_names:
-1305            d_extracted[name] = []
-1306            ens_length = len(new_idl_d[name])
-1307            for i_dat, dat in enumerate(data):
-1308                d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name], _compute_scalefactor_missing_rep(o).get(name.split('|')[0], 1)) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
-1309        for name in new_cov_names:
-1310            g_extracted[name] = []
-1311            zero_grad = _Zero_grad(new_covobs_lengths[name])
-1312            for i_dat, dat in enumerate(data):
-1313                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)))
-1314
-1315    for i_val, new_val in np.ndenumerate(new_values):
-1316        new_deltas = {}
-1317        new_grad = {}
-1318        if array_mode is True:
-1319            for name in new_sample_names:
-1320                ens_length = d_extracted[name][0].shape[-1]
-1321                new_deltas[name] = np.zeros(ens_length)
-1322                for i_dat, dat in enumerate(d_extracted[name]):
-1323                    new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1324            for name in new_cov_names:
-1325                new_grad[name] = 0
-1326                for i_dat, dat in enumerate(g_extracted[name]):
-1327                    new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1328        else:
-1329            for j_obs, obs in np.ndenumerate(data):
-1330                scalef_d = _compute_scalefactor_missing_rep(obs)
-1331                for name in obs.names:
-1332                    if name in obs.cov_names:
-1333                        new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
-1334                    else:
-1335                        new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name], scalef_d.get(name.split('|')[0], 1))
-1336
-1337        new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
-1338
-1339        if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
-1340            raise Exception('The same name has been used for deltas and covobs!')
-1341        new_samples = []
-1342        new_means = []
-1343        new_idl = []
-1344        new_names_obs = []
-1345        for name in new_names:
-1346            if name not in new_covobs:
-1347                new_samples.append(new_deltas[name])
-1348                new_idl.append(new_idl_d[name])
-1349                new_means.append(new_r_values[name][i_val])
-1350                new_names_obs.append(name)
-1351        final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
-1352        for name in new_covobs:
-1353            final_result[i_val].names.append(name)
-1354        final_result[i_val]._covobs = new_covobs
-1355        final_result[i_val]._value = new_val
-1356        final_result[i_val].reweighted = reweighted
+1203    # Workaround for matrix operations containing non Obs data
+1204    if not all(isinstance(x, Obs) for x in raveled_data):
+1205        for i in range(len(raveled_data)):
+1206            if isinstance(raveled_data[i], (int, float)):
+1207                raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
+1208
+1209    allcov = {}
+1210    for o in raveled_data:
+1211        for name in o.cov_names:
+1212            if name in allcov:
+1213                if not np.allclose(allcov[name], o.covobs[name].cov):
+1214                    raise Exception('Inconsistent covariance matrices for %s!' % (name))
+1215            else:
+1216                allcov[name] = o.covobs[name].cov
+1217
+1218    n_obs = len(raveled_data)
+1219    new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
+1220    new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
+1221    new_sample_names = sorted(set(new_names) - set(new_cov_names))
+1222
+1223    reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
+1224
+1225    if data.ndim == 1:
+1226        values = np.array([o.value for o in data])
+1227    else:
+1228        values = np.vectorize(lambda x: x.value)(data)
+1229
+1230    new_values = func(values, **kwargs)
+1231
+1232    multi = int(isinstance(new_values, np.ndarray))
+1233
+1234    new_r_values = {}
+1235    new_idl_d = {}
+1236    for name in new_sample_names:
+1237        idl = []
+1238        tmp_values = np.zeros(n_obs)
+1239        for i, item in enumerate(raveled_data):
+1240            tmp_values[i] = item.r_values.get(name, item.value)
+1241            tmp_idl = item.idl.get(name)
+1242            if tmp_idl is not None:
+1243                idl.append(tmp_idl)
+1244        if multi > 0:
+1245            tmp_values = np.array(tmp_values).reshape(data.shape)
+1246        new_r_values[name] = func(tmp_values, **kwargs)
+1247        new_idl_d[name] = _merge_idx(idl)
+1248
+1249    def _compute_scalefactor_missing_rep(obs):
+1250        """
+1251        Computes the scale factor that is to be multiplied with the deltas
+1252        in the case where Obs with different subsets of replica are merged.
+1253        Returns a dictionary with the scale factor for each Monte Carlo name.
+1254
+1255        Parameters
+1256        ----------
+1257        obs : Obs
+1258            The observable corresponding to the deltas that are to be scaled
+1259        """
+1260        scalef_d = {}
+1261        for mc_name in obs.mc_names:
+1262            mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')]
+1263            new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')]
+1264            if len(mc_idl_d) > 0 and len(mc_idl_d) < len(new_mc_idl_d):
+1265                scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d])
+1266        return scalef_d
+1267
+1268    if 'man_grad' in kwargs:
+1269        deriv = np.asarray(kwargs.get('man_grad'))
+1270        if new_values.shape + data.shape != deriv.shape:
+1271            raise Exception('Manual derivative does not have correct shape.')
+1272    elif kwargs.get('num_grad') is True:
+1273        if multi > 0:
+1274            raise Exception('Multi mode currently not supported for numerical derivative')
+1275        options = {
+1276            'base_step': 0.1,
+1277            'step_ratio': 2.5}
+1278        for key in options.keys():
+1279            kwarg = kwargs.get(key)
+1280            if kwarg is not None:
+1281                options[key] = kwarg
+1282        tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
+1283        if tmp_df.size == 1:
+1284            deriv = np.array([tmp_df.real])
+1285        else:
+1286            deriv = tmp_df.real
+1287    else:
+1288        deriv = jacobian(func)(values, **kwargs)
+1289
+1290    final_result = np.zeros(new_values.shape, dtype=object)
+1291
+1292    if array_mode is True:
+1293
+1294        class _Zero_grad():
+1295            def __init__(self, N):
+1296                self.grad = np.zeros((N, 1))
+1297
+1298        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]))
+1299        d_extracted = {}
+1300        g_extracted = {}
+1301        for name in new_sample_names:
+1302            d_extracted[name] = []
+1303            ens_length = len(new_idl_d[name])
+1304            for i_dat, dat in enumerate(data):
+1305                d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name], _compute_scalefactor_missing_rep(o).get(name.split('|')[0], 1)) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
+1306        for name in new_cov_names:
+1307            g_extracted[name] = []
+1308            zero_grad = _Zero_grad(new_covobs_lengths[name])
+1309            for i_dat, dat in enumerate(data):
+1310                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)))
+1311
+1312    for i_val, new_val in np.ndenumerate(new_values):
+1313        new_deltas = {}
+1314        new_grad = {}
+1315        if array_mode is True:
+1316            for name in new_sample_names:
+1317                ens_length = d_extracted[name][0].shape[-1]
+1318                new_deltas[name] = np.zeros(ens_length)
+1319                for i_dat, dat in enumerate(d_extracted[name]):
+1320                    new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1321            for name in new_cov_names:
+1322                new_grad[name] = 0
+1323                for i_dat, dat in enumerate(g_extracted[name]):
+1324                    new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1325        else:
+1326            for j_obs, obs in np.ndenumerate(data):
+1327                scalef_d = _compute_scalefactor_missing_rep(obs)
+1328                for name in obs.names:
+1329                    if name in obs.cov_names:
+1330                        new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
+1331                    else:
+1332                        new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name], scalef_d.get(name.split('|')[0], 1))
+1333
+1334        new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
+1335
+1336        if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
+1337            raise Exception('The same name has been used for deltas and covobs!')
+1338        new_samples = []
+1339        new_means = []
+1340        new_idl = []
+1341        new_names_obs = []
+1342        for name in new_names:
+1343            if name not in new_covobs:
+1344                new_samples.append(new_deltas[name])
+1345                new_idl.append(new_idl_d[name])
+1346                new_means.append(new_r_values[name][i_val])
+1347                new_names_obs.append(name)
+1348        final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
+1349        for name in new_covobs:
+1350            final_result[i_val].names.append(name)
+1351        final_result[i_val]._covobs = new_covobs
+1352        final_result[i_val]._value = new_val
+1353        final_result[i_val].reweighted = reweighted
+1354
+1355    if multi == 0:
+1356        final_result = final_result.item()
 1357
-1358    if multi == 0:
-1359        final_result = final_result.item()
-1360
-1361    return final_result
+1358    return final_result
 

 
 
@@ -5628,46 +5622,46 @@ functions. For the ratio of two observables one can e.g. use
- 
1393def reweight(weight, obs, **kwargs):
-1394    """Reweight a list of observables.
-1395
-1396    Parameters
-1397    ----------
-1398    weight : Obs
-1399        Reweighting factor. An Observable that has to be defined on a superset of the
-1400        configurations in obs[i].idl for all i.
-1401    obs : list
-1402        list of Obs, e.g. [obs1, obs2, obs3].
-1403    all_configs : bool
-1404        if True, the reweighted observables are normalized by the average of
-1405        the reweighting factor on all configurations in weight.idl and not
-1406        on the configurations in obs[i].idl. Default False.
-1407    """
-1408    result = []
-1409    for i in range(len(obs)):
-1410        if len(obs[i].cov_names):
-1411            raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
-1412        if not set(obs[i].names).issubset(weight.names):
-1413            raise Exception('Error: Ensembles do not fit')
-1414        for name in obs[i].names:
-1415            if not set(obs[i].idl[name]).issubset(weight.idl[name]):
-1416                raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
-1417        new_samples = []
-1418        w_deltas = {}
-1419        for name in sorted(obs[i].names):
-1420            w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
-1421            new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
-1422        tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
-1423
-1424        if kwargs.get('all_configs'):
-1425            new_weight = weight
-1426        else:
-1427            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)])
+            
1390def reweight(weight, obs, **kwargs):
+1391    """Reweight a list of observables.
+1392
+1393    Parameters
+1394    ----------
+1395    weight : Obs
+1396        Reweighting factor. An Observable that has to be defined on a superset of the
+1397        configurations in obs[i].idl for all i.
+1398    obs : list
+1399        list of Obs, e.g. [obs1, obs2, obs3].
+1400    all_configs : bool
+1401        if True, the reweighted observables are normalized by the average of
+1402        the reweighting factor on all configurations in weight.idl and not
+1403        on the configurations in obs[i].idl. Default False.
+1404    """
+1405    result = []
+1406    for i in range(len(obs)):
+1407        if len(obs[i].cov_names):
+1408            raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
+1409        if not set(obs[i].names).issubset(weight.names):
+1410            raise Exception('Error: Ensembles do not fit')
+1411        for name in obs[i].names:
+1412            if not set(obs[i].idl[name]).issubset(weight.idl[name]):
+1413                raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
+1414        new_samples = []
+1415        w_deltas = {}
+1416        for name in sorted(obs[i].names):
+1417            w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
+1418            new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
+1419        tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
+1420
+1421        if kwargs.get('all_configs'):
+1422            new_weight = weight
+1423        else:
+1424            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)])
+1425
+1426        result.append(tmp_obs / new_weight)
+1427        result[-1].reweighted = True
 1428
-1429        result.append(tmp_obs / new_weight)
-1430        result[-1].reweighted = True
-1431
-1432    return result
+1429    return result
 

 
 
@@ -5701,47 +5695,47 @@ on the configurations in obs[i].idl. Default False.
- 
1435def correlate(obs_a, obs_b):
-1436    """Correlate two observables.
-1437
-1438    Parameters
-1439    ----------
-1440    obs_a : Obs
-1441        First observable
-1442    obs_b : Obs
-1443        Second observable
-1444
-1445    Notes
-1446    -----
-1447    Keep in mind to only correlate primary observables which have not been reweighted
-1448    yet. The reweighting has to be applied after correlating the observables.
-1449    Currently only works if ensembles are identical (this is not strictly necessary).
-1450    """
-1451
-1452    if sorted(obs_a.names) != sorted(obs_b.names):
-1453        raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
-1454    if len(obs_a.cov_names) or len(obs_b.cov_names):
-1455        raise Exception('Error: Not possible to correlate Obs that contain covobs!')
-1456    for name in obs_a.names:
-1457        if obs_a.shape[name] != obs_b.shape[name]:
-1458            raise Exception('Shapes of ensemble', name, 'do not fit')
-1459        if obs_a.idl[name] != obs_b.idl[name]:
-1460            raise Exception('idl of ensemble', name, 'do not fit')
-1461
-1462    if obs_a.reweighted is True:
-1463        warnings.warn("The first observable is already reweighted.", RuntimeWarning)
-1464    if obs_b.reweighted is True:
-1465        warnings.warn("The second observable is already reweighted.", RuntimeWarning)
-1466
-1467    new_samples = []
-1468    new_idl = []
-1469    for name in sorted(obs_a.names):
-1470        new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
-1471        new_idl.append(obs_a.idl[name])
-1472
-1473    o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
-1474    o.reweighted = obs_a.reweighted or obs_b.reweighted
-1475    return o
+            
1432def correlate(obs_a, obs_b):
+1433    """Correlate two observables.
+1434
+1435    Parameters
+1436    ----------
+1437    obs_a : Obs
+1438        First observable
+1439    obs_b : Obs
+1440        Second observable
+1441
+1442    Notes
+1443    -----
+1444    Keep in mind to only correlate primary observables which have not been reweighted
+1445    yet. The reweighting has to be applied after correlating the observables.
+1446    Currently only works if ensembles are identical (this is not strictly necessary).
+1447    """
+1448
+1449    if sorted(obs_a.names) != sorted(obs_b.names):
+1450        raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
+1451    if len(obs_a.cov_names) or len(obs_b.cov_names):
+1452        raise Exception('Error: Not possible to correlate Obs that contain covobs!')
+1453    for name in obs_a.names:
+1454        if obs_a.shape[name] != obs_b.shape[name]:
+1455            raise Exception('Shapes of ensemble', name, 'do not fit')
+1456        if obs_a.idl[name] != obs_b.idl[name]:
+1457            raise Exception('idl of ensemble', name, 'do not fit')
+1458
+1459    if obs_a.reweighted is True:
+1460        warnings.warn("The first observable is already reweighted.", RuntimeWarning)
+1461    if obs_b.reweighted is True:
+1462        warnings.warn("The second observable is already reweighted.", RuntimeWarning)
+1463
+1464    new_samples = []
+1465    new_idl = []
+1466    for name in sorted(obs_a.names):
+1467        new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
+1468        new_idl.append(obs_a.idl[name])
+1469
+1470    o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
+1471    o.reweighted = obs_a.reweighted or obs_b.reweighted
+1472    return o
 

 
 
@@ -5776,74 +5770,74 @@ Currently only works if ensembles are identical (this is not strictly necessary)
- 
1478def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
-1479    r'''Calculates the error covariance matrix of a set of observables.
+            
1475def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
+1476    r'''Calculates the error covariance matrix of a set of observables.
+1477
+1478    WARNING: This function should be used with care, especially for observables with support on multiple
+1479             ensembles with differing autocorrelations. See the notes below for details.
 1480
-1481    WARNING: This function should be used with care, especially for observables with support on multiple
-1482             ensembles with differing autocorrelations. See the notes below for details.
-1483
-1484    The gamma method has to be applied first to all observables.
-1485
-1486    Parameters
-1487    ----------
-1488    obs : list or numpy.ndarray
-1489        List or one dimensional array of Obs
-1490    visualize : bool
-1491        If True plots the corresponding normalized correlation matrix (default False).
-1492    correlation : bool
-1493        If True the correlation matrix instead of the error covariance matrix is returned (default False).
-1494    smooth : None or int
-1495        If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
-1496        smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
-1497        largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
-1498        small ones.
-1499
-1500    Notes
-1501    -----
-1502    The error covariance is defined such that it agrees with the squared standard error for two identical observables
-1503    $$\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$$
-1504    in the absence of autocorrelation.
-1505    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
-1506    $$\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.
-1507    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.
-1508    $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
-1509    This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
-1510    '''
-1511
-1512    length = len(obs)
-1513
-1514    max_samples = np.max([o.N for o in obs])
-1515    if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
-1516        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)
-1517
-1518    cov = np.zeros((length, length))
-1519    for i in range(length):
-1520        for j in range(i, length):
-1521            cov[i, j] = _covariance_element(obs[i], obs[j])
-1522    cov = cov + cov.T - np.diag(np.diag(cov))
-1523
-1524    corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+1481    The gamma method has to be applied first to all observables.
+1482
+1483    Parameters
+1484    ----------
+1485    obs : list or numpy.ndarray
+1486        List or one dimensional array of Obs
+1487    visualize : bool
+1488        If True plots the corresponding normalized correlation matrix (default False).
+1489    correlation : bool
+1490        If True the correlation matrix instead of the error covariance matrix is returned (default False).
+1491    smooth : None or int
+1492        If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
+1493        smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
+1494        largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
+1495        small ones.
+1496
+1497    Notes
+1498    -----
+1499    The error covariance is defined such that it agrees with the squared standard error for two identical observables
+1500    $$\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$$
+1501    in the absence of autocorrelation.
+1502    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
+1503    $$\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.
+1504    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.
+1505    $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
+1506    This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
+1507    '''
+1508
+1509    length = len(obs)
+1510
+1511    max_samples = np.max([o.N for o in obs])
+1512    if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
+1513        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)
+1514
+1515    cov = np.zeros((length, length))
+1516    for i in range(length):
+1517        for j in range(i, length):
+1518            cov[i, j] = _covariance_element(obs[i], obs[j])
+1519    cov = cov + cov.T - np.diag(np.diag(cov))
+1520
+1521    corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+1522
+1523    if isinstance(smooth, int):
+1524        corr = _smooth_eigenvalues(corr, smooth)
 1525
-1526    if isinstance(smooth, int):
-1527        corr = _smooth_eigenvalues(corr, smooth)
-1528
-1529    if visualize:
-1530        plt.matshow(corr, vmin=-1, vmax=1)
-1531        plt.set_cmap('RdBu')
-1532        plt.colorbar()
-1533        plt.draw()
+1526    if visualize:
+1527        plt.matshow(corr, vmin=-1, vmax=1)
+1528        plt.set_cmap('RdBu')
+1529        plt.colorbar()
+1530        plt.draw()
+1531
+1532    if correlation is True:
+1533        return corr
 1534
-1535    if correlation is True:
-1536        return corr
+1535    errors = [o.dvalue for o in obs]
+1536    cov = np.diag(errors) @ corr @ np.diag(errors)
 1537
-1538    errors = [o.dvalue for o in obs]
-1539    cov = np.diag(errors) @ corr @ np.diag(errors)
-1540
-1541    eigenvalues = np.linalg.eigh(cov)[0]
-1542    if not np.all(eigenvalues >= 0):
-1543        warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
-1544
-1545    return cov
+1538    eigenvalues = np.linalg.eigh(cov)[0]
+1539    if not np.all(eigenvalues >= 0):
+1540        warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
+1541
+1542    return cov
 

 
 
@@ -5895,27 +5889,27 @@ This construction ensures that the estimated covariance matrix is positive semi-
- 
1548def invert_corr_cov_cholesky(corr, inverrdiag):
-1549    """Constructs a lower triangular matrix `chol` via the Cholesky decomposition of the correlation matrix `corr`
-1550       and then returns the inverse covariance matrix `chol_inv` as a lower triangular matrix by solving `chol * x = inverrdiag`.
-1551
-1552    Parameters
-1553    ----------
-1554    corr : np.ndarray
-1555           correlation matrix
-1556    inverrdiag : np.ndarray
-1557              diagonal matrix, the entries are the inverse errors of the data points considered
-1558    """
-1559
-1560    condn = np.linalg.cond(corr)
-1561    if condn > 0.1 / np.finfo(float).eps:
-1562        raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
-1563    if condn > 1e13:
-1564        warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
-1565    chol = np.linalg.cholesky(corr)
-1566    chol_inv = scipy.linalg.solve_triangular(chol, inverrdiag, lower=True)
-1567
-1568    return chol_inv
+            
1545def invert_corr_cov_cholesky(corr, inverrdiag):
+1546    """Constructs a lower triangular matrix `chol` via the Cholesky decomposition of the correlation matrix `corr`
+1547       and then returns the inverse covariance matrix `chol_inv` as a lower triangular matrix by solving `chol * x = inverrdiag`.
+1548
+1549    Parameters
+1550    ----------
+1551    corr : np.ndarray
+1552           correlation matrix
+1553    inverrdiag : np.ndarray
+1554              diagonal matrix, the entries are the inverse errors of the data points considered
+1555    """
+1556
+1557    condn = np.linalg.cond(corr)
+1558    if condn > 0.1 / np.finfo(float).eps:
+1559        raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
+1560    if condn > 1e13:
+1561        warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
+1562    chol = np.linalg.cholesky(corr)
+1563    chol_inv = scipy.linalg.solve_triangular(chol, inverrdiag, lower=True)
+1564
+1565    return chol_inv
 

 
 
@@ -5945,67 +5939,67 @@ diagonal matrix, the entries are the inverse errors of the data points considere
- 
1571def sort_corr(corr, kl, yd):
-1572    """ Reorders a correlation matrix to match the alphabetical order of its underlying y data.
-1573
-1574    The ordering of the input correlation matrix `corr` is given by the list of keys `kl`.
-1575    The input dictionary `yd` (with the same keys `kl`) must contain the corresponding y data
-1576    that the correlation matrix is based on.
-1577    This function sorts the list of keys `kl` alphabetically and sorts the matrix `corr`
-1578    according to this alphabetical order such that the sorted matrix `corr_sorted` corresponds
-1579    to the y data `yd` when arranged in an alphabetical order by its keys.
-1580
-1581    Parameters
-1582    ----------
-1583    corr : np.ndarray
-1584        A square correlation matrix constructed using the order of the y data specified by `kl`.
-1585        The dimensions of `corr` should match the total number of y data points in `yd` combined.
-1586    kl : list of str
-1587        A list of keys that denotes the order in which the y data from `yd` was used to build the
-1588        input correlation matrix `corr`.
-1589    yd : dict of list
-1590        A dictionary where each key corresponds to a unique identifier, and its value is a list of
-1591        y data points. The total number of y data points across all keys must match the dimensions
-1592        of `corr`. The lists in the dictionary can be lists of Obs.
-1593
-1594    Returns
-1595    -------
-1596    np.ndarray
-1597        A new, sorted correlation matrix that corresponds to the y data from `yd` when arranged alphabetically by its keys.
-1598
-1599    Example
-1600    -------
-1601    >>> import numpy as np
-1602    >>> import pyerrors as pe
-1603    >>> corr = np.array([[1, 0.2, 0.3], [0.2, 1, 0.4], [0.3, 0.4, 1]])
-1604    >>> kl = ['b', 'a']
-1605    >>> yd = {'a': [1, 2], 'b': [3]}
-1606    >>> sorted_corr = pe.obs.sort_corr(corr, kl, yd)
-1607    >>> print(sorted_corr)
-1608    array([[1. , 0.3, 0.4],
-1609           [0.3, 1. , 0.2],
-1610           [0.4, 0.2, 1. ]])
+            
1568def sort_corr(corr, kl, yd):
+1569    """ Reorders a correlation matrix to match the alphabetical order of its underlying y data.
+1570
+1571    The ordering of the input correlation matrix `corr` is given by the list of keys `kl`.
+1572    The input dictionary `yd` (with the same keys `kl`) must contain the corresponding y data
+1573    that the correlation matrix is based on.
+1574    This function sorts the list of keys `kl` alphabetically and sorts the matrix `corr`
+1575    according to this alphabetical order such that the sorted matrix `corr_sorted` corresponds
+1576    to the y data `yd` when arranged in an alphabetical order by its keys.
+1577
+1578    Parameters
+1579    ----------
+1580    corr : np.ndarray
+1581        A square correlation matrix constructed using the order of the y data specified by `kl`.
+1582        The dimensions of `corr` should match the total number of y data points in `yd` combined.
+1583    kl : list of str
+1584        A list of keys that denotes the order in which the y data from `yd` was used to build the
+1585        input correlation matrix `corr`.
+1586    yd : dict of list
+1587        A dictionary where each key corresponds to a unique identifier, and its value is a list of
+1588        y data points. The total number of y data points across all keys must match the dimensions
+1589        of `corr`. The lists in the dictionary can be lists of Obs.
+1590
+1591    Returns
+1592    -------
+1593    np.ndarray
+1594        A new, sorted correlation matrix that corresponds to the y data from `yd` when arranged alphabetically by its keys.
+1595
+1596    Example
+1597    -------
+1598    >>> import numpy as np
+1599    >>> import pyerrors as pe
+1600    >>> corr = np.array([[1, 0.2, 0.3], [0.2, 1, 0.4], [0.3, 0.4, 1]])
+1601    >>> kl = ['b', 'a']
+1602    >>> yd = {'a': [1, 2], 'b': [3]}
+1603    >>> sorted_corr = pe.obs.sort_corr(corr, kl, yd)
+1604    >>> print(sorted_corr)
+1605    array([[1. , 0.3, 0.4],
+1606           [0.3, 1. , 0.2],
+1607           [0.4, 0.2, 1. ]])
+1608
+1609    """
+1610    kl_sorted = sorted(kl)
 1611
-1612    """
-1613    kl_sorted = sorted(kl)
-1614
-1615    posd = {}
-1616    ofs = 0
-1617    for ki, k in enumerate(kl):
-1618        posd[k] = [i + ofs for i in range(len(yd[k]))]
-1619        ofs += len(posd[k])
-1620
-1621    mapping = []
-1622    for k in kl_sorted:
-1623        for i in range(len(yd[k])):
-1624            mapping.append(posd[k][i])
-1625
-1626    corr_sorted = np.zeros_like(corr)
-1627    for i in range(corr.shape[0]):
-1628        for j in range(corr.shape[0]):
-1629            corr_sorted[i][j] = corr[mapping[i]][mapping[j]]
-1630
-1631    return corr_sorted
+1612    posd = {}
+1613    ofs = 0
+1614    for ki, k in enumerate(kl):
+1615        posd[k] = [i + ofs for i in range(len(yd[k]))]
+1616        ofs += len(posd[k])
+1617
+1618    mapping = []
+1619    for k in kl_sorted:
+1620        for i in range(len(yd[k])):
+1621            mapping.append(posd[k][i])
+1622
+1623    corr_sorted = np.zeros_like(corr)
+1624    for i in range(corr.shape[0]):
+1625        for j in range(corr.shape[0]):
+1626            corr_sorted[i][j] = corr[mapping[i]][mapping[j]]
+1627
+1628    return corr_sorted
 

 
 
@@ -6069,24 +6063,24 @@ of corr. The lists in the dictionary can be lists of Obs.
- 
1711def import_jackknife(jacks, name, idl=None):
-1712    """Imports jackknife samples and returns an Obs
-1713
-1714    Parameters
-1715    ----------
-1716    jacks : numpy.ndarray
-1717        numpy array containing the mean value as zeroth entry and
-1718        the N jackknife samples as first to Nth entry.
-1719    name : str
-1720        name of the ensemble the samples are defined on.
-1721    """
-1722    length = len(jacks) - 1
-1723    prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
-1724    samples = jacks[1:] @ prj
-1725    mean = np.mean(samples)
-1726    new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
-1727    new_obs._value = jacks[0]
-1728    return new_obs
+            
1708def import_jackknife(jacks, name, idl=None):
+1709    """Imports jackknife samples and returns an Obs
+1710
+1711    Parameters
+1712    ----------
+1713    jacks : numpy.ndarray
+1714        numpy array containing the mean value as zeroth entry and
+1715        the N jackknife samples as first to Nth entry.
+1716    name : str
+1717        name of the ensemble the samples are defined on.
+1718    """
+1719    length = len(jacks) - 1
+1720    prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
+1721    samples = jacks[1:] @ prj
+1722    mean = np.mean(samples)
+1723    new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
+1724    new_obs._value = jacks[0]
+1725    return new_obs
 

 
 
@@ -6116,34 +6110,34 @@ name of the ensemble the samples are defined on.
- 
1731def import_bootstrap(boots, name, random_numbers):
-1732    """Imports bootstrap samples and returns an Obs
-1733
-1734    Parameters
-1735    ----------
-1736    boots : numpy.ndarray
-1737        numpy array containing the mean value as zeroth entry and
-1738        the N bootstrap samples as first to Nth entry.
-1739    name : str
-1740        name of the ensemble the samples are defined on.
-1741    random_numbers : np.ndarray
-1742        Array of shape (samples, length) containing the random numbers to generate the bootstrap samples,
-1743        where samples is the number of bootstrap samples and length is the length of the original Monte Carlo
-1744        chain to be reconstructed.
-1745    """
-1746    samples, length = random_numbers.shape
-1747    if samples != len(boots) - 1:
-1748        raise ValueError("Random numbers do not have the correct shape.")
+            
1728def import_bootstrap(boots, name, random_numbers):
+1729    """Imports bootstrap samples and returns an Obs
+1730
+1731    Parameters
+1732    ----------
+1733    boots : numpy.ndarray
+1734        numpy array containing the mean value as zeroth entry and
+1735        the N bootstrap samples as first to Nth entry.
+1736    name : str
+1737        name of the ensemble the samples are defined on.
+1738    random_numbers : np.ndarray
+1739        Array of shape (samples, length) containing the random numbers to generate the bootstrap samples,
+1740        where samples is the number of bootstrap samples and length is the length of the original Monte Carlo
+1741        chain to be reconstructed.
+1742    """
+1743    samples, length = random_numbers.shape
+1744    if samples != len(boots) - 1:
+1745        raise ValueError("Random numbers do not have the correct shape.")
+1746
+1747    if samples < length:
+1748        raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.")
 1749
-1750    if samples < length:
-1751        raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.")
-1752
-1753    proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
-1754
-1755    samples = scipy.linalg.lstsq(proj, boots[1:])[0]
-1756    ret = Obs([samples], [name])
-1757    ret._value = boots[0]
-1758    return ret
+1750    proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
+1751
+1752    samples = scipy.linalg.lstsq(proj, boots[1:])[0]
+1753    ret = Obs([samples], [name])
+1754    ret._value = boots[0]
+1755    return ret
 

 
 
@@ -6177,34 +6171,34 @@ chain to be reconstructed.
- 
1761def merge_obs(list_of_obs):
-1762    """Combine all observables in list_of_obs into one new observable
-1763
-1764    Parameters
-1765    ----------
-1766    list_of_obs : list
-1767        list of the Obs object to be combined
-1768
-1769    Notes
-1770    -----
-1771    It is not possible to combine obs which are based on the same replicum
-1772    """
-1773    replist = [item for obs in list_of_obs for item in obs.names]
-1774    if (len(replist) == len(set(replist))) is False:
-1775        raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
-1776    if any([len(o.cov_names) for o in list_of_obs]):
-1777        raise Exception('Not possible to merge data that contains covobs!')
-1778    new_dict = {}
-1779    idl_dict = {}
-1780    for o in list_of_obs:
-1781        new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
-1782                        for key in set(o.deltas) | set(o.r_values)})
-1783        idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
-1784
-1785    names = sorted(new_dict.keys())
-1786    o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
-1787    o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
-1788    return o
+            
1758def merge_obs(list_of_obs):
+1759    """Combine all observables in list_of_obs into one new observable
+1760
+1761    Parameters
+1762    ----------
+1763    list_of_obs : list
+1764        list of the Obs object to be combined
+1765
+1766    Notes
+1767    -----
+1768    It is not possible to combine obs which are based on the same replicum
+1769    """
+1770    replist = [item for obs in list_of_obs for item in obs.names]
+1771    if (len(replist) == len(set(replist))) is False:
+1772        raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
+1773    if any([len(o.cov_names) for o in list_of_obs]):
+1774        raise Exception('Not possible to merge data that contains covobs!')
+1775    new_dict = {}
+1776    idl_dict = {}
+1777    for o in list_of_obs:
+1778        new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
+1779                        for key in set(o.deltas) | set(o.r_values)})
+1780        idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
+1781
+1782    names = sorted(new_dict.keys())
+1783    o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
+1784    o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
+1785    return o
 

 
 
@@ -6235,47 +6229,47 @@ list of the Obs object to be combined
- 
1791def cov_Obs(means, cov, name, grad=None):
-1792    """Create an Obs based on mean(s) and a covariance matrix
-1793
-1794    Parameters
-1795    ----------
-1796    mean : list of floats or float
-1797        N mean value(s) of the new Obs
-1798    cov : list or array
-1799        2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
-1800    name : str
-1801        identifier for the covariance matrix
-1802    grad : list or array
-1803        Gradient of the Covobs wrt. the means belonging to cov.
-1804    """
+            
1788def cov_Obs(means, cov, name, grad=None):
+1789    """Create an Obs based on mean(s) and a covariance matrix
+1790
+1791    Parameters
+1792    ----------
+1793    mean : list of floats or float
+1794        N mean value(s) of the new Obs
+1795    cov : list or array
+1796        2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
+1797    name : str
+1798        identifier for the covariance matrix
+1799    grad : list or array
+1800        Gradient of the Covobs wrt. the means belonging to cov.
+1801    """
+1802
+1803    def covobs_to_obs(co):
+1804        """Make an Obs out of a Covobs
 1805
-1806    def covobs_to_obs(co):
-1807        """Make an Obs out of a Covobs
-1808
-1809        Parameters
-1810        ----------
-1811        co : Covobs
-1812            Covobs to be embedded into the Obs
-1813        """
-1814        o = Obs([], [], means=[])
-1815        o._value = co.value
-1816        o.names.append(co.name)
-1817        o._covobs[co.name] = co
-1818        o._dvalue = np.sqrt(co.errsq())
-1819        return o
-1820
-1821    ol = []
-1822    if isinstance(means, (float, int)):
-1823        means = [means]
-1824
-1825    for i in range(len(means)):
-1826        ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
-1827    if ol[0].covobs[name].N != len(means):
-1828        raise Exception('You have to provide %d mean values!' % (ol[0].N))
-1829    if len(ol) == 1:
-1830        return ol[0]
-1831    return ol
+1806        Parameters
+1807        ----------
+1808        co : Covobs
+1809            Covobs to be embedded into the Obs
+1810        """
+1811        o = Obs([], [], means=[])
+1812        o._value = co.value
+1813        o.names.append(co.name)
+1814        o._covobs[co.name] = co
+1815        o._dvalue = np.sqrt(co.errsq())
+1816        return o
+1817
+1818    ol = []
+1819    if isinstance(means, (float, int)):
+1820        means = [means]
+1821
+1822    for i in range(len(means)):
+1823        ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
+1824    if ol[0].covobs[name].N != len(means):
+1825        raise Exception('You have to provide %d mean values!' % (ol[0].N))
+1826    if len(ol) == 1:
+1827        return ol[0]
+1828    return ol