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)])
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])
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)])
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)])
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)])
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])
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])
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])
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])
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)])
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)])
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])
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])
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])
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])
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)"
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
929 @property
-930 def real(self):
-931 return self._real
+ 926 @property
+927 def real(self):
+928 return self._real
933 @property
-934 def imag(self):
-935 return self._imag
+ 930 @property
+931 def imag(self):
+932 return self._imag
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)
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
948 def conjugate(self):
-949 return CObs(self.real, -self.imag)
+ 945 def conjugate(self):
+946 return CObs(self.real, -self.imag)
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)
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)
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
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
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
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
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
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
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
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
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
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