diff --git a/docs/pyerrors/obs.html b/docs/pyerrors/obs.html
index 08587452..1d790a77 100644
--- a/docs/pyerrors/obs.html
+++ b/docs/pyerrors/obs.html
@@ -902,948 +902,961 @@
695 def __str__(self):
696 return _format_uncertainty(self.value, self._dvalue)
697
- 698 def __hash__(self):
- 699 hash_tuple = (np.array([self.value]).astype(np.float32).data.tobytes(),)
- 700 hash_tuple += tuple([o.astype(np.float32).data.tobytes() for o in self.deltas.values()])
- 701 hash_tuple += tuple([np.array([o.errsq()]).astype(np.float32).data.tobytes() for o in self.covobs.values()])
- 702 hash_tuple += tuple([o.encode() for o in self.names])
- 703 m = hashlib.md5()
- 704 [m.update(o) for o in hash_tuple]
- 705 return int(m.hexdigest(), 16) & 0xFFFFFFFF
+ 698 def __format__(self, format_type):
+ 699 my_str = _format_uncertainty(self.value, self._dvalue,
+ 700 significance=int(float(format_type.replace("+", "").replace("-", ""))))
+ 701 for char in ["+", " "]:
+ 702 if format_type.startswith(char):
+ 703 if my_str[0] != "-":
+ 704 my_str = char + my_str
+ 705 return my_str
706
- 707 # Overload comparisons
- 708 def __lt__(self, other):
- 709 return self.value < other
- 710
- 711 def __le__(self, other):
- 712 return self.value <= other
- 713
- 714 def __gt__(self, other):
- 715 return self.value > other
- 716
- 717 def __ge__(self, other):
- 718 return self.value >= other
+ 707 def __hash__(self):
+ 708 hash_tuple = (np.array([self.value]).astype(np.float32).data.tobytes(),)
+ 709 hash_tuple += tuple([o.astype(np.float32).data.tobytes() for o in self.deltas.values()])
+ 710 hash_tuple += tuple([np.array([o.errsq()]).astype(np.float32).data.tobytes() for o in self.covobs.values()])
+ 711 hash_tuple += tuple([o.encode() for o in self.names])
+ 712 m = hashlib.md5()
+ 713 [m.update(o) for o in hash_tuple]
+ 714 return int(m.hexdigest(), 16) & 0xFFFFFFFF
+ 715
+ 716 # Overload comparisons
+ 717 def __lt__(self, other):
+ 718 return self.value < other
719
- 720 def __eq__(self, other):
- 721 return (self - other).is_zero()
+ 720 def __le__(self, other):
+ 721 return self.value <= other
722
- 723 def __ne__(self, other):
- 724 return not (self - other).is_zero()
+ 723 def __gt__(self, other):
+ 724 return self.value > other
725
- 726 # Overload math operations
- 727 def __add__(self, y):
- 728 if isinstance(y, Obs):
- 729 return derived_observable(lambda x, **kwargs: x[0] + x[1], [self, y], man_grad=[1, 1])
- 730 else:
- 731 if isinstance(y, np.ndarray):
- 732 return np.array([self + o for o in y])
- 733 elif y.__class__.__name__ in ['Corr', 'CObs']:
- 734 return NotImplemented
- 735 else:
- 736 return derived_observable(lambda x, **kwargs: x[0] + y, [self], man_grad=[1])
- 737
- 738 def __radd__(self, y):
- 739 return self + y
- 740
- 741 def __mul__(self, y):
- 742 if isinstance(y, Obs):
- 743 return derived_observable(lambda x, **kwargs: x[0] * x[1], [self, y], man_grad=[y.value, self.value])
- 744 else:
- 745 if isinstance(y, np.ndarray):
- 746 return np.array([self * o for o in y])
- 747 elif isinstance(y, complex):
- 748 return CObs(self * y.real, self * y.imag)
- 749 elif y.__class__.__name__ in ['Corr', 'CObs']:
- 750 return NotImplemented
- 751 else:
- 752 return derived_observable(lambda x, **kwargs: x[0] * y, [self], man_grad=[y])
- 753
- 754 def __rmul__(self, y):
- 755 return self * y
- 756
- 757 def __sub__(self, y):
- 758 if isinstance(y, Obs):
- 759 return derived_observable(lambda x, **kwargs: x[0] - x[1], [self, y], man_grad=[1, -1])
- 760 else:
- 761 if isinstance(y, np.ndarray):
- 762 return np.array([self - o for o in y])
- 763 elif y.__class__.__name__ in ['Corr', 'CObs']:
- 764 return NotImplemented
- 765 else:
- 766 return derived_observable(lambda x, **kwargs: x[0] - y, [self], man_grad=[1])
- 767
- 768 def __rsub__(self, y):
- 769 return -1 * (self - y)
- 770
- 771 def __pos__(self):
- 772 return self
- 773
- 774 def __neg__(self):
- 775 return -1 * self
+ 726 def __ge__(self, other):
+ 727 return self.value >= other
+ 728
+ 729 def __eq__(self, other):
+ 730 return (self - other).is_zero()
+ 731
+ 732 def __ne__(self, other):
+ 733 return not (self - other).is_zero()
+ 734
+ 735 # Overload math operations
+ 736 def __add__(self, y):
+ 737 if isinstance(y, Obs):
+ 738 return derived_observable(lambda x, **kwargs: x[0] + x[1], [self, y], man_grad=[1, 1])
+ 739 else:
+ 740 if isinstance(y, np.ndarray):
+ 741 return np.array([self + o for o in y])
+ 742 elif y.__class__.__name__ in ['Corr', 'CObs']:
+ 743 return NotImplemented
+ 744 else:
+ 745 return derived_observable(lambda x, **kwargs: x[0] + y, [self], man_grad=[1])
+ 746
+ 747 def __radd__(self, y):
+ 748 return self + y
+ 749
+ 750 def __mul__(self, y):
+ 751 if isinstance(y, Obs):
+ 752 return derived_observable(lambda x, **kwargs: x[0] * x[1], [self, y], man_grad=[y.value, self.value])
+ 753 else:
+ 754 if isinstance(y, np.ndarray):
+ 755 return np.array([self * o for o in y])
+ 756 elif isinstance(y, complex):
+ 757 return CObs(self * y.real, self * y.imag)
+ 758 elif y.__class__.__name__ in ['Corr', 'CObs']:
+ 759 return NotImplemented
+ 760 else:
+ 761 return derived_observable(lambda x, **kwargs: x[0] * y, [self], man_grad=[y])
+ 762
+ 763 def __rmul__(self, y):
+ 764 return self * y
+ 765
+ 766 def __sub__(self, y):
+ 767 if isinstance(y, Obs):
+ 768 return derived_observable(lambda x, **kwargs: x[0] - x[1], [self, y], man_grad=[1, -1])
+ 769 else:
+ 770 if isinstance(y, np.ndarray):
+ 771 return np.array([self - o for o in y])
+ 772 elif y.__class__.__name__ in ['Corr', 'CObs']:
+ 773 return NotImplemented
+ 774 else:
+ 775 return derived_observable(lambda x, **kwargs: x[0] - y, [self], man_grad=[1])
776
- 777 def __truediv__(self, y):
- 778 if isinstance(y, Obs):
- 779 return derived_observable(lambda x, **kwargs: x[0] / x[1], [self, y], man_grad=[1 / y.value, - self.value / y.value ** 2])
- 780 else:
- 781 if isinstance(y, np.ndarray):
- 782 return np.array([self / o for o in y])
- 783 elif y.__class__.__name__ in ['Corr', 'CObs']:
- 784 return NotImplemented
- 785 else:
- 786 return derived_observable(lambda x, **kwargs: x[0] / y, [self], man_grad=[1 / y])
- 787
- 788 def __rtruediv__(self, y):
- 789 if isinstance(y, Obs):
- 790 return derived_observable(lambda x, **kwargs: x[0] / x[1], [y, self], man_grad=[1 / self.value, - y.value / self.value ** 2])
- 791 else:
- 792 if isinstance(y, np.ndarray):
- 793 return np.array([o / self for o in y])
- 794 elif y.__class__.__name__ in ['Corr', 'CObs']:
- 795 return NotImplemented
- 796 else:
- 797 return derived_observable(lambda x, **kwargs: y / x[0], [self], man_grad=[-y / self.value ** 2])
- 798
- 799 def __pow__(self, y):
- 800 if isinstance(y, Obs):
- 801 return derived_observable(lambda x: x[0] ** x[1], [self, y])
- 802 else:
- 803 return derived_observable(lambda x: x[0] ** y, [self])
- 804
- 805 def __rpow__(self, y):
- 806 if isinstance(y, Obs):
- 807 return derived_observable(lambda x: x[0] ** x[1], [y, self])
- 808 else:
- 809 return derived_observable(lambda x: y ** x[0], [self])
- 810
- 811 def __abs__(self):
- 812 return derived_observable(lambda x: anp.abs(x[0]), [self])
+ 777 def __rsub__(self, y):
+ 778 return -1 * (self - y)
+ 779
+ 780 def __pos__(self):
+ 781 return self
+ 782
+ 783 def __neg__(self):
+ 784 return -1 * self
+ 785
+ 786 def __truediv__(self, y):
+ 787 if isinstance(y, Obs):
+ 788 return derived_observable(lambda x, **kwargs: x[0] / x[1], [self, y], man_grad=[1 / y.value, - self.value / y.value ** 2])
+ 789 else:
+ 790 if isinstance(y, np.ndarray):
+ 791 return np.array([self / o for o in y])
+ 792 elif y.__class__.__name__ in ['Corr', 'CObs']:
+ 793 return NotImplemented
+ 794 else:
+ 795 return derived_observable(lambda x, **kwargs: x[0] / y, [self], man_grad=[1 / y])
+ 796
+ 797 def __rtruediv__(self, y):
+ 798 if isinstance(y, Obs):
+ 799 return derived_observable(lambda x, **kwargs: x[0] / x[1], [y, self], man_grad=[1 / self.value, - y.value / self.value ** 2])
+ 800 else:
+ 801 if isinstance(y, np.ndarray):
+ 802 return np.array([o / self for o in y])
+ 803 elif y.__class__.__name__ in ['Corr', 'CObs']:
+ 804 return NotImplemented
+ 805 else:
+ 806 return derived_observable(lambda x, **kwargs: y / x[0], [self], man_grad=[-y / self.value ** 2])
+ 807
+ 808 def __pow__(self, y):
+ 809 if isinstance(y, Obs):
+ 810 return derived_observable(lambda x: x[0] ** x[1], [self, y])
+ 811 else:
+ 812 return derived_observable(lambda x: x[0] ** y, [self])
813
- 814 # Overload numpy functions
- 815 def sqrt(self):
- 816 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
- 817
- 818 def log(self):
- 819 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
- 820
- 821 def exp(self):
- 822 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
- 823
- 824 def sin(self):
- 825 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
+ 814 def __rpow__(self, y):
+ 815 if isinstance(y, Obs):
+ 816 return derived_observable(lambda x: x[0] ** x[1], [y, self])
+ 817 else:
+ 818 return derived_observable(lambda x: y ** x[0], [self])
+ 819
+ 820 def __abs__(self):
+ 821 return derived_observable(lambda x: anp.abs(x[0]), [self])
+ 822
+ 823 # Overload numpy functions
+ 824 def sqrt(self):
+ 825 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
826
- 827 def cos(self):
- 828 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
+ 827 def log(self):
+ 828 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
829
- 830 def tan(self):
- 831 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
+ 830 def exp(self):
+ 831 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
832
- 833 def arcsin(self):
- 834 return derived_observable(lambda x: anp.arcsin(x[0]), [self])
+ 833 def sin(self):
+ 834 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
835
- 836 def arccos(self):
- 837 return derived_observable(lambda x: anp.arccos(x[0]), [self])
+ 836 def cos(self):
+ 837 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
838
- 839 def arctan(self):
- 840 return derived_observable(lambda x: anp.arctan(x[0]), [self])
+ 839 def tan(self):
+ 840 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
841
- 842 def sinh(self):
- 843 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
+ 842 def arcsin(self):
+ 843 return derived_observable(lambda x: anp.arcsin(x[0]), [self])
844
- 845 def cosh(self):
- 846 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
+ 845 def arccos(self):
+ 846 return derived_observable(lambda x: anp.arccos(x[0]), [self])
847
- 848 def tanh(self):
- 849 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
+ 848 def arctan(self):
+ 849 return derived_observable(lambda x: anp.arctan(x[0]), [self])
850
- 851 def arcsinh(self):
- 852 return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
+ 851 def sinh(self):
+ 852 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
853
- 854 def arccosh(self):
- 855 return derived_observable(lambda x: anp.arccosh(x[0]), [self])
+ 854 def cosh(self):
+ 855 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
856
- 857 def arctanh(self):
- 858 return derived_observable(lambda x: anp.arctanh(x[0]), [self])
+ 857 def tanh(self):
+ 858 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
859
- 860
- 861class CObs:
- 862 """Class for a complex valued observable."""
- 863 __slots__ = ['_real', '_imag', 'tag']
- 864
- 865 def __init__(self, real, imag=0.0):
- 866 self._real = real
- 867 self._imag = imag
- 868 self.tag = None
+ 860 def arcsinh(self):
+ 861 return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
+ 862
+ 863 def arccosh(self):
+ 864 return derived_observable(lambda x: anp.arccosh(x[0]), [self])
+ 865
+ 866 def arctanh(self):
+ 867 return derived_observable(lambda x: anp.arctanh(x[0]), [self])
+ 868
869
- 870 @property
- 871 def real(self):
- 872 return self._real
+ 870class CObs:
+ 871 """Class for a complex valued observable."""
+ 872 __slots__ = ['_real', '_imag', 'tag']
873
- 874 @property
- 875 def imag(self):
- 876 return self._imag
- 877
- 878 def gamma_method(self, **kwargs):
- 879 """Executes the gamma_method for the real and the imaginary part."""
- 880 if isinstance(self.real, Obs):
- 881 self.real.gamma_method(**kwargs)
- 882 if isinstance(self.imag, Obs):
- 883 self.imag.gamma_method(**kwargs)
- 884
- 885 def is_zero(self):
- 886 """Checks whether both real and imaginary part are zero within machine precision."""
- 887 return self.real == 0.0 and self.imag == 0.0
- 888
- 889 def conjugate(self):
- 890 return CObs(self.real, -self.imag)
- 891
- 892 def __add__(self, other):
- 893 if isinstance(other, np.ndarray):
- 894 return other + self
- 895 elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 896 return CObs(self.real + other.real,
- 897 self.imag + other.imag)
- 898 else:
- 899 return CObs(self.real + other, self.imag)
+ 874 def __init__(self, real, imag=0.0):
+ 875 self._real = real
+ 876 self._imag = imag
+ 877 self.tag = None
+ 878
+ 879 @property
+ 880 def real(self):
+ 881 return self._real
+ 882
+ 883 @property
+ 884 def imag(self):
+ 885 return self._imag
+ 886
+ 887 def gamma_method(self, **kwargs):
+ 888 """Executes the gamma_method for the real and the imaginary part."""
+ 889 if isinstance(self.real, Obs):
+ 890 self.real.gamma_method(**kwargs)
+ 891 if isinstance(self.imag, Obs):
+ 892 self.imag.gamma_method(**kwargs)
+ 893
+ 894 def is_zero(self):
+ 895 """Checks whether both real and imaginary part are zero within machine precision."""
+ 896 return self.real == 0.0 and self.imag == 0.0
+ 897
+ 898 def conjugate(self):
+ 899 return CObs(self.real, -self.imag)
900
- 901 def __radd__(self, y):
- 902 return self + y
- 903
- 904 def __sub__(self, other):
- 905 if isinstance(other, np.ndarray):
- 906 return -1 * (other - self)
- 907 elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 908 return CObs(self.real - other.real, self.imag - other.imag)
- 909 else:
- 910 return CObs(self.real - other, self.imag)
- 911
- 912 def __rsub__(self, other):
- 913 return -1 * (self - other)
- 914
- 915 def __mul__(self, other):
- 916 if isinstance(other, np.ndarray):
- 917 return other * self
- 918 elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 919 if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
- 920 return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
- 921 [self.real, other.real, self.imag, other.imag],
- 922 man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
- 923 derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
- 924 [self.real, other.real, self.imag, other.imag],
- 925 man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
- 926 elif getattr(other, 'imag', 0) != 0:
- 927 return CObs(self.real * other.real - self.imag * other.imag,
- 928 self.imag * other.real + self.real * other.imag)
- 929 else:
- 930 return CObs(self.real * other.real, self.imag * other.real)
- 931 else:
- 932 return CObs(self.real * other, self.imag * other)
- 933
- 934 def __rmul__(self, other):
- 935 return self * other
- 936
- 937 def __truediv__(self, other):
- 938 if isinstance(other, np.ndarray):
- 939 return 1 / (other / self)
- 940 elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 941 r = other.real ** 2 + other.imag ** 2
- 942 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
- 943 else:
- 944 return CObs(self.real / other, self.imag / other)
+ 901 def __add__(self, other):
+ 902 if isinstance(other, np.ndarray):
+ 903 return other + self
+ 904 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 905 return CObs(self.real + other.real,
+ 906 self.imag + other.imag)
+ 907 else:
+ 908 return CObs(self.real + other, self.imag)
+ 909
+ 910 def __radd__(self, y):
+ 911 return self + y
+ 912
+ 913 def __sub__(self, other):
+ 914 if isinstance(other, np.ndarray):
+ 915 return -1 * (other - self)
+ 916 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 917 return CObs(self.real - other.real, self.imag - other.imag)
+ 918 else:
+ 919 return CObs(self.real - other, self.imag)
+ 920
+ 921 def __rsub__(self, other):
+ 922 return -1 * (self - other)
+ 923
+ 924 def __mul__(self, other):
+ 925 if isinstance(other, np.ndarray):
+ 926 return other * self
+ 927 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 928 if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
+ 929 return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
+ 930 [self.real, other.real, self.imag, other.imag],
+ 931 man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
+ 932 derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
+ 933 [self.real, other.real, self.imag, other.imag],
+ 934 man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
+ 935 elif getattr(other, 'imag', 0) != 0:
+ 936 return CObs(self.real * other.real - self.imag * other.imag,
+ 937 self.imag * other.real + self.real * other.imag)
+ 938 else:
+ 939 return CObs(self.real * other.real, self.imag * other.real)
+ 940 else:
+ 941 return CObs(self.real * other, self.imag * other)
+ 942
+ 943 def __rmul__(self, other):
+ 944 return self * other
945
- 946 def __rtruediv__(self, other):
- 947 r = self.real ** 2 + self.imag ** 2
- 948 if hasattr(other, 'real') and hasattr(other, 'imag'):
- 949 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
- 950 else:
- 951 return CObs(self.real * other / r, -self.imag * other / r)
- 952
- 953 def __abs__(self):
- 954 return np.sqrt(self.real**2 + self.imag**2)
- 955
- 956 def __pos__(self):
- 957 return self
- 958
- 959 def __neg__(self):
- 960 return -1 * self
+ 946 def __truediv__(self, other):
+ 947 if isinstance(other, np.ndarray):
+ 948 return 1 / (other / self)
+ 949 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 950 r = other.real ** 2 + other.imag ** 2
+ 951 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
+ 952 else:
+ 953 return CObs(self.real / other, self.imag / other)
+ 954
+ 955 def __rtruediv__(self, other):
+ 956 r = self.real ** 2 + self.imag ** 2
+ 957 if hasattr(other, 'real') and hasattr(other, 'imag'):
+ 958 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
+ 959 else:
+ 960 return CObs(self.real * other / r, -self.imag * other / r)
961
- 962 def __eq__(self, other):
- 963 return self.real == other.real and self.imag == other.imag
+ 962 def __abs__(self):
+ 963 return np.sqrt(self.real**2 + self.imag**2)
964
- 965 def __str__(self):
- 966 return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
+ 965 def __pos__(self):
+ 966 return self
967
- 968 def __repr__(self):
- 969 return 'CObs[' + str(self) + ']'
+ 968 def __neg__(self):
+ 969 return -1 * self
970
- 971
- 972def _format_uncertainty(value, dvalue):
- 973 """Creates a string of a value and its error in paranthesis notation, e.g., 13.02(45)"""
- 974 if dvalue == 0.0:
- 975 return str(value)
- 976 fexp = np.floor(np.log10(dvalue))
- 977 if fexp < 0.0:
- 978 return '{:{form}}({:2.0f})'.format(value, dvalue * 10 ** (-fexp + 1), form='.' + str(-int(fexp) + 1) + 'f')
- 979 elif fexp == 0.0:
- 980 return '{:.1f}({:1.1f})'.format(value, dvalue)
- 981 else:
- 982 return '{:.0f}({:2.0f})'.format(value, dvalue)
- 983
- 984
- 985def _expand_deltas(deltas, idx, shape, gapsize):
- 986 """Expand deltas defined on idx to a regular range with spacing gapsize between two
- 987 configurations and where holes are filled by 0.
- 988 If idx is of type range, the deltas are not changed if the idx.step == gapsize.
- 989
- 990 Parameters
- 991 ----------
- 992 deltas : list
- 993 List of fluctuations
- 994 idx : list
- 995 List or range of configs on which the deltas are defined, has to be sorted in ascending order.
- 996 shape : int
- 997 Number of configs in idx.
- 998 gapsize : int
- 999 The target distance between two configurations. If longer distances
-1000 are found in idx, the data is expanded.
-1001 """
-1002 if isinstance(idx, range):
-1003 if (idx.step == gapsize):
-1004 return deltas
-1005 ret = np.zeros((idx[-1] - idx[0] + gapsize) // gapsize)
-1006 for i in range(shape):
-1007 ret[(idx[i] - idx[0]) // gapsize] = deltas[i]
-1008 return ret
-1009
-1010
-1011def _merge_idx(idl):
-1012 """Returns the union of all lists in idl as sorted list
-1013
-1014 Parameters
-1015 ----------
-1016 idl : list
-1017 List of lists or ranges.
-1018 """
-1019
-1020 # Use groupby to efficiently check whether all elements of idl are identical
-1021 try:
-1022 g = groupby(idl)
-1023 if next(g, True) and not next(g, False):
-1024 return idl[0]
-1025 except Exception:
-1026 pass
-1027
-1028 if np.all([type(idx) is range for idx in idl]):
-1029 if len(set([idx[0] for idx in idl])) == 1:
-1030 idstart = min([idx.start for idx in idl])
-1031 idstop = max([idx.stop for idx in idl])
-1032 idstep = min([idx.step for idx in idl])
-1033 return range(idstart, idstop, idstep)
-1034
-1035 return sorted(set().union(*idl))
-1036
-1037
-1038def _intersection_idx(idl):
-1039 """Returns the intersection of all lists in idl as sorted list
+ 971 def __eq__(self, other):
+ 972 return self.real == other.real and self.imag == other.imag
+ 973
+ 974 def __str__(self):
+ 975 return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
+ 976
+ 977 def __repr__(self):
+ 978 return 'CObs[' + str(self) + ']'
+ 979
+ 980
+ 981def _format_uncertainty(value, dvalue, significance=2):
+ 982 """Creates a string of a value and its error in paranthesis notation, e.g., 13.02(45)"""
+ 983 if dvalue == 0.0:
+ 984 return str(value)
+ 985 if not isinstance(significance, int):
+ 986 raise TypeError("significance needs to be an integer.")
+ 987 if significance < 1:
+ 988 raise ValueError("significance needs to be larger than zero.")
+ 989 fexp = np.floor(np.log10(dvalue))
+ 990 if fexp < 0.0:
+ 991 return '{:{form}}({:1.0f})'.format(value, dvalue * 10 ** (-fexp + significance - 1), form='.' + str(-int(fexp) + significance - 1) + 'f')
+ 992 elif fexp == 0.0:
+ 993 return f"{value:.{significance - 1}f}({dvalue:1.{significance - 1}f})"
+ 994 else:
+ 995 return f"{value:.{max(0, int(significance - fexp - 1))}f}({dvalue:2.{max(0, int(significance - fexp - 1))}f})"
+ 996
+ 997
+ 998def _expand_deltas(deltas, idx, shape, gapsize):
+ 999 """Expand deltas defined on idx to a regular range with spacing gapsize between two
+1000 configurations and where holes are filled by 0.
+1001 If idx is of type range, the deltas are not changed if the idx.step == gapsize.
+1002
+1003 Parameters
+1004 ----------
+1005 deltas : list
+1006 List of fluctuations
+1007 idx : list
+1008 List or range of configs on which the deltas are defined, has to be sorted in ascending order.
+1009 shape : int
+1010 Number of configs in idx.
+1011 gapsize : int
+1012 The target distance between two configurations. If longer distances
+1013 are found in idx, the data is expanded.
+1014 """
+1015 if isinstance(idx, range):
+1016 if (idx.step == gapsize):
+1017 return deltas
+1018 ret = np.zeros((idx[-1] - idx[0] + gapsize) // gapsize)
+1019 for i in range(shape):
+1020 ret[(idx[i] - idx[0]) // gapsize] = deltas[i]
+1021 return ret
+1022
+1023
+1024def _merge_idx(idl):
+1025 """Returns the union of all lists in idl as sorted list
+1026
+1027 Parameters
+1028 ----------
+1029 idl : list
+1030 List of lists or ranges.
+1031 """
+1032
+1033 # Use groupby to efficiently check whether all elements of idl are identical
+1034 try:
+1035 g = groupby(idl)
+1036 if next(g, True) and not next(g, False):
+1037 return idl[0]
+1038 except Exception:
+1039 pass
1040
-1041 Parameters
-1042 ----------
-1043 idl : list
-1044 List of lists or ranges.
-1045 """
-1046
-1047 def _lcm(*args):
-1048 """Returns the lowest common multiple of args.
+1041 if np.all([type(idx) is range for idx in idl]):
+1042 if len(set([idx[0] for idx in idl])) == 1:
+1043 idstart = min([idx.start for idx in idl])
+1044 idstop = max([idx.stop for idx in idl])
+1045 idstep = min([idx.step for idx in idl])
+1046 return range(idstart, idstop, idstep)
+1047
+1048 return sorted(set().union(*idl))
1049
-1050 From python 3.9 onwards the math library contains an lcm function."""
-1051 return reduce(lambda a, b: a * b // gcd(a, b), args)
-1052
-1053 # Use groupby to efficiently check whether all elements of idl are identical
-1054 try:
-1055 g = groupby(idl)
-1056 if next(g, True) and not next(g, False):
-1057 return idl[0]
-1058 except Exception:
-1059 pass
-1060
-1061 if np.all([type(idx) is range for idx in idl]):
-1062 if len(set([idx[0] for idx in idl])) == 1:
-1063 idstart = max([idx.start for idx in idl])
-1064 idstop = min([idx.stop for idx in idl])
-1065 idstep = _lcm(*[idx.step for idx in idl])
-1066 return range(idstart, idstop, idstep)
-1067
-1068 return sorted(set.intersection(*[set(o) for o in idl]))
-1069
-1070
-1071def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
-1072 """Expand deltas defined on idx to the list of configs that is defined by new_idx.
-1073 New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest
-1074 common divisor of the step sizes is used as new step size.
-1075
-1076 Parameters
-1077 ----------
-1078 deltas : list
-1079 List of fluctuations
-1080 idx : list
-1081 List or range of configs on which the deltas are defined.
-1082 Has to be a subset of new_idx and has to be sorted in ascending order.
-1083 shape : list
-1084 Number of configs in idx.
-1085 new_idx : list
-1086 List of configs that defines the new range, has to be sorted in ascending order.
-1087 """
+1050
+1051def _intersection_idx(idl):
+1052 """Returns the intersection of all lists in idl as sorted list
+1053
+1054 Parameters
+1055 ----------
+1056 idl : list
+1057 List of lists or ranges.
+1058 """
+1059
+1060 def _lcm(*args):
+1061 """Returns the lowest common multiple of args.
+1062
+1063 From python 3.9 onwards the math library contains an lcm function."""
+1064 return reduce(lambda a, b: a * b // gcd(a, b), args)
+1065
+1066 # Use groupby to efficiently check whether all elements of idl are identical
+1067 try:
+1068 g = groupby(idl)
+1069 if next(g, True) and not next(g, False):
+1070 return idl[0]
+1071 except Exception:
+1072 pass
+1073
+1074 if np.all([type(idx) is range for idx in idl]):
+1075 if len(set([idx[0] for idx in idl])) == 1:
+1076 idstart = max([idx.start for idx in idl])
+1077 idstop = min([idx.stop for idx in idl])
+1078 idstep = _lcm(*[idx.step for idx in idl])
+1079 return range(idstart, idstop, idstep)
+1080
+1081 return sorted(set.intersection(*[set(o) for o in idl]))
+1082
+1083
+1084def _expand_deltas_for_merge(deltas, idx, shape, new_idx):
+1085 """Expand deltas defined on idx to the list of configs that is defined by new_idx.
+1086 New, empty entries are filled by 0. If idx and new_idx are of type range, the smallest
+1087 common divisor of the step sizes is used as new step size.
1088
-1089 if type(idx) is range and type(new_idx) is range:
-1090 if idx == new_idx:
-1091 return deltas
-1092 ret = np.zeros(new_idx[-1] - new_idx[0] + 1)
-1093 for i in range(shape):
-1094 ret[idx[i] - new_idx[0]] = deltas[i]
-1095 return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx)
-1096
-1097
-1098def derived_observable(func, data, array_mode=False, **kwargs):
-1099 """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
-1100
-1101 Parameters
-1102 ----------
-1103 func : object
-1104 arbitrary function of the form func(data, **kwargs). For the
-1105 automatic differentiation to work, all numpy functions have to have
-1106 the autograd wrapper (use 'import autograd.numpy as anp').
-1107 data : list
-1108 list of Obs, e.g. [obs1, obs2, obs3].
-1109 num_grad : bool
-1110 if True, numerical derivatives are used instead of autograd
-1111 (default False). To control the numerical differentiation the
-1112 kwargs of numdifftools.step_generators.MaxStepGenerator
-1113 can be used.
-1114 man_grad : list
-1115 manually supply a list or an array which contains the jacobian
-1116 of func. Use cautiously, supplying the wrong derivative will
-1117 not be intercepted.
-1118
-1119 Notes
-1120 -----
-1121 For simple mathematical operations it can be practical to use anonymous
-1122 functions. For the ratio of two observables one can e.g. use
-1123
-1124 new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
-1125 """
-1126
-1127 data = np.asarray(data)
-1128 raveled_data = data.ravel()
-1129
-1130 # Workaround for matrix operations containing non Obs data
-1131 if not all(isinstance(x, Obs) for x in raveled_data):
-1132 for i in range(len(raveled_data)):
-1133 if isinstance(raveled_data[i], (int, float)):
-1134 raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
-1135
-1136 allcov = {}
-1137 for o in raveled_data:
-1138 for name in o.cov_names:
-1139 if name in allcov:
-1140 if not np.allclose(allcov[name], o.covobs[name].cov):
-1141 raise Exception('Inconsistent covariance matrices for %s!' % (name))
-1142 else:
-1143 allcov[name] = o.covobs[name].cov
-1144
-1145 n_obs = len(raveled_data)
-1146 new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
-1147 new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
-1148 new_sample_names = sorted(set(new_names) - set(new_cov_names))
-1149
-1150 reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
-1151
-1152 if data.ndim == 1:
-1153 values = np.array([o.value for o in data])
-1154 else:
-1155 values = np.vectorize(lambda x: x.value)(data)
-1156
-1157 new_values = func(values, **kwargs)
-1158
-1159 multi = int(isinstance(new_values, np.ndarray))
-1160
-1161 new_r_values = {}
-1162 new_idl_d = {}
-1163 for name in new_sample_names:
-1164 idl = []
-1165 tmp_values = np.zeros(n_obs)
-1166 for i, item in enumerate(raveled_data):
-1167 tmp_values[i] = item.r_values.get(name, item.value)
-1168 tmp_idl = item.idl.get(name)
-1169 if tmp_idl is not None:
-1170 idl.append(tmp_idl)
-1171 if multi > 0:
-1172 tmp_values = np.array(tmp_values).reshape(data.shape)
-1173 new_r_values[name] = func(tmp_values, **kwargs)
-1174 new_idl_d[name] = _merge_idx(idl)
-1175
-1176 if 'man_grad' in kwargs:
-1177 deriv = np.asarray(kwargs.get('man_grad'))
-1178 if new_values.shape + data.shape != deriv.shape:
-1179 raise Exception('Manual derivative does not have correct shape.')
-1180 elif kwargs.get('num_grad') is True:
-1181 if multi > 0:
-1182 raise Exception('Multi mode currently not supported for numerical derivative')
-1183 options = {
-1184 'base_step': 0.1,
-1185 'step_ratio': 2.5}
-1186 for key in options.keys():
-1187 kwarg = kwargs.get(key)
-1188 if kwarg is not None:
-1189 options[key] = kwarg
-1190 tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
-1191 if tmp_df.size == 1:
-1192 deriv = np.array([tmp_df.real])
-1193 else:
-1194 deriv = tmp_df.real
-1195 else:
-1196 deriv = jacobian(func)(values, **kwargs)
-1197
-1198 final_result = np.zeros(new_values.shape, dtype=object)
-1199
-1200 if array_mode is True:
-1201
-1202 class _Zero_grad():
-1203 def __init__(self, N):
-1204 self.grad = np.zeros((N, 1))
-1205
-1206 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]))
-1207 d_extracted = {}
-1208 g_extracted = {}
-1209 for name in new_sample_names:
-1210 d_extracted[name] = []
-1211 ens_length = len(new_idl_d[name])
-1212 for i_dat, dat in enumerate(data):
-1213 d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name]) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
-1214 for name in new_cov_names:
-1215 g_extracted[name] = []
-1216 zero_grad = _Zero_grad(new_covobs_lengths[name])
-1217 for i_dat, dat in enumerate(data):
-1218 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)))
-1219
-1220 for i_val, new_val in np.ndenumerate(new_values):
-1221 new_deltas = {}
-1222 new_grad = {}
-1223 if array_mode is True:
-1224 for name in new_sample_names:
-1225 ens_length = d_extracted[name][0].shape[-1]
-1226 new_deltas[name] = np.zeros(ens_length)
-1227 for i_dat, dat in enumerate(d_extracted[name]):
-1228 new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1229 for name in new_cov_names:
-1230 new_grad[name] = 0
-1231 for i_dat, dat in enumerate(g_extracted[name]):
-1232 new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1233 else:
-1234 for j_obs, obs in np.ndenumerate(data):
-1235 for name in obs.names:
-1236 if name in obs.cov_names:
-1237 new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
-1238 else:
-1239 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])
-1240
-1241 new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
-1242
-1243 if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
-1244 raise Exception('The same name has been used for deltas and covobs!')
-1245 new_samples = []
-1246 new_means = []
-1247 new_idl = []
-1248 new_names_obs = []
-1249 for name in new_names:
-1250 if name not in new_covobs:
-1251 new_samples.append(new_deltas[name])
-1252 new_idl.append(new_idl_d[name])
-1253 new_means.append(new_r_values[name][i_val])
-1254 new_names_obs.append(name)
-1255 final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
-1256 for name in new_covobs:
-1257 final_result[i_val].names.append(name)
-1258 final_result[i_val]._covobs = new_covobs
-1259 final_result[i_val]._value = new_val
-1260 final_result[i_val].reweighted = reweighted
-1261
-1262 if multi == 0:
-1263 final_result = final_result.item()
-1264
-1265 return final_result
-1266
-1267
-1268def _reduce_deltas(deltas, idx_old, idx_new):
-1269 """Extract deltas defined on idx_old on all configs of idx_new.
-1270
-1271 Assumes, that idx_old and idx_new are correctly defined idl, i.e., they
-1272 are ordered in an ascending order.
-1273
-1274 Parameters
-1275 ----------
-1276 deltas : list
-1277 List of fluctuations
-1278 idx_old : list
-1279 List or range of configs on which the deltas are defined
-1280 idx_new : list
-1281 List of configs for which we want to extract the deltas.
-1282 Has to be a subset of idx_old.
-1283 """
-1284 if not len(deltas) == len(idx_old):
-1285 raise Exception('Length of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old)))
-1286 if type(idx_old) is range and type(idx_new) is range:
-1287 if idx_old == idx_new:
-1288 return deltas
-1289 # Use groupby to efficiently check whether all elements of idx_old and idx_new are identical
-1290 try:
-1291 g = groupby([idx_old, idx_new])
-1292 if next(g, True) and not next(g, False):
-1293 return deltas
-1294 except Exception:
-1295 pass
-1296 indices = np.intersect1d(idx_old, idx_new, assume_unique=True, return_indices=True)[1]
-1297 if len(indices) < len(idx_new):
-1298 raise Exception('Error in _reduce_deltas: Config of idx_new not in idx_old')
-1299 return np.array(deltas)[indices]
-1300
-1301
-1302def reweight(weight, obs, **kwargs):
-1303 """Reweight a list of observables.
-1304
-1305 Parameters
-1306 ----------
-1307 weight : Obs
-1308 Reweighting factor. An Observable that has to be defined on a superset of the
-1309 configurations in obs[i].idl for all i.
-1310 obs : list
-1311 list of Obs, e.g. [obs1, obs2, obs3].
-1312 all_configs : bool
-1313 if True, the reweighted observables are normalized by the average of
-1314 the reweighting factor on all configurations in weight.idl and not
-1315 on the configurations in obs[i].idl. Default False.
-1316 """
-1317 result = []
-1318 for i in range(len(obs)):
-1319 if len(obs[i].cov_names):
-1320 raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
-1321 if not set(obs[i].names).issubset(weight.names):
-1322 raise Exception('Error: Ensembles do not fit')
-1323 for name in obs[i].names:
-1324 if not set(obs[i].idl[name]).issubset(weight.idl[name]):
-1325 raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
-1326 new_samples = []
-1327 w_deltas = {}
-1328 for name in sorted(obs[i].names):
-1329 w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
-1330 new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
-1331 tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
-1332
-1333 if kwargs.get('all_configs'):
-1334 new_weight = weight
-1335 else:
-1336 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)])
-1337
-1338 result.append(tmp_obs / new_weight)
-1339 result[-1].reweighted = True
-1340
-1341 return result
-1342
-1343
-1344def correlate(obs_a, obs_b):
-1345 """Correlate two observables.
-1346
-1347 Parameters
-1348 ----------
-1349 obs_a : Obs
-1350 First observable
-1351 obs_b : Obs
-1352 Second observable
+1089 Parameters
+1090 ----------
+1091 deltas : list
+1092 List of fluctuations
+1093 idx : list
+1094 List or range of configs on which the deltas are defined.
+1095 Has to be a subset of new_idx and has to be sorted in ascending order.
+1096 shape : list
+1097 Number of configs in idx.
+1098 new_idx : list
+1099 List of configs that defines the new range, has to be sorted in ascending order.
+1100 """
+1101
+1102 if type(idx) is range and type(new_idx) is range:
+1103 if idx == new_idx:
+1104 return deltas
+1105 ret = np.zeros(new_idx[-1] - new_idx[0] + 1)
+1106 for i in range(shape):
+1107 ret[idx[i] - new_idx[0]] = deltas[i]
+1108 return np.array([ret[new_idx[i] - new_idx[0]] for i in range(len(new_idx))]) * len(new_idx) / len(idx)
+1109
+1110
+1111def derived_observable(func, data, array_mode=False, **kwargs):
+1112 """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
+1113
+1114 Parameters
+1115 ----------
+1116 func : object
+1117 arbitrary function of the form func(data, **kwargs). For the
+1118 automatic differentiation to work, all numpy functions have to have
+1119 the autograd wrapper (use 'import autograd.numpy as anp').
+1120 data : list
+1121 list of Obs, e.g. [obs1, obs2, obs3].
+1122 num_grad : bool
+1123 if True, numerical derivatives are used instead of autograd
+1124 (default False). To control the numerical differentiation the
+1125 kwargs of numdifftools.step_generators.MaxStepGenerator
+1126 can be used.
+1127 man_grad : list
+1128 manually supply a list or an array which contains the jacobian
+1129 of func. Use cautiously, supplying the wrong derivative will
+1130 not be intercepted.
+1131
+1132 Notes
+1133 -----
+1134 For simple mathematical operations it can be practical to use anonymous
+1135 functions. For the ratio of two observables one can e.g. use
+1136
+1137 new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
+1138 """
+1139
+1140 data = np.asarray(data)
+1141 raveled_data = data.ravel()
+1142
+1143 # Workaround for matrix operations containing non Obs data
+1144 if not all(isinstance(x, Obs) for x in raveled_data):
+1145 for i in range(len(raveled_data)):
+1146 if isinstance(raveled_data[i], (int, float)):
+1147 raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
+1148
+1149 allcov = {}
+1150 for o in raveled_data:
+1151 for name in o.cov_names:
+1152 if name in allcov:
+1153 if not np.allclose(allcov[name], o.covobs[name].cov):
+1154 raise Exception('Inconsistent covariance matrices for %s!' % (name))
+1155 else:
+1156 allcov[name] = o.covobs[name].cov
+1157
+1158 n_obs = len(raveled_data)
+1159 new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
+1160 new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
+1161 new_sample_names = sorted(set(new_names) - set(new_cov_names))
+1162
+1163 reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
+1164
+1165 if data.ndim == 1:
+1166 values = np.array([o.value for o in data])
+1167 else:
+1168 values = np.vectorize(lambda x: x.value)(data)
+1169
+1170 new_values = func(values, **kwargs)
+1171
+1172 multi = int(isinstance(new_values, np.ndarray))
+1173
+1174 new_r_values = {}
+1175 new_idl_d = {}
+1176 for name in new_sample_names:
+1177 idl = []
+1178 tmp_values = np.zeros(n_obs)
+1179 for i, item in enumerate(raveled_data):
+1180 tmp_values[i] = item.r_values.get(name, item.value)
+1181 tmp_idl = item.idl.get(name)
+1182 if tmp_idl is not None:
+1183 idl.append(tmp_idl)
+1184 if multi > 0:
+1185 tmp_values = np.array(tmp_values).reshape(data.shape)
+1186 new_r_values[name] = func(tmp_values, **kwargs)
+1187 new_idl_d[name] = _merge_idx(idl)
+1188
+1189 if 'man_grad' in kwargs:
+1190 deriv = np.asarray(kwargs.get('man_grad'))
+1191 if new_values.shape + data.shape != deriv.shape:
+1192 raise Exception('Manual derivative does not have correct shape.')
+1193 elif kwargs.get('num_grad') is True:
+1194 if multi > 0:
+1195 raise Exception('Multi mode currently not supported for numerical derivative')
+1196 options = {
+1197 'base_step': 0.1,
+1198 'step_ratio': 2.5}
+1199 for key in options.keys():
+1200 kwarg = kwargs.get(key)
+1201 if kwarg is not None:
+1202 options[key] = kwarg
+1203 tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
+1204 if tmp_df.size == 1:
+1205 deriv = np.array([tmp_df.real])
+1206 else:
+1207 deriv = tmp_df.real
+1208 else:
+1209 deriv = jacobian(func)(values, **kwargs)
+1210
+1211 final_result = np.zeros(new_values.shape, dtype=object)
+1212
+1213 if array_mode is True:
+1214
+1215 class _Zero_grad():
+1216 def __init__(self, N):
+1217 self.grad = np.zeros((N, 1))
+1218
+1219 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]))
+1220 d_extracted = {}
+1221 g_extracted = {}
+1222 for name in new_sample_names:
+1223 d_extracted[name] = []
+1224 ens_length = len(new_idl_d[name])
+1225 for i_dat, dat in enumerate(data):
+1226 d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name]) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
+1227 for name in new_cov_names:
+1228 g_extracted[name] = []
+1229 zero_grad = _Zero_grad(new_covobs_lengths[name])
+1230 for i_dat, dat in enumerate(data):
+1231 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)))
+1232
+1233 for i_val, new_val in np.ndenumerate(new_values):
+1234 new_deltas = {}
+1235 new_grad = {}
+1236 if array_mode is True:
+1237 for name in new_sample_names:
+1238 ens_length = d_extracted[name][0].shape[-1]
+1239 new_deltas[name] = np.zeros(ens_length)
+1240 for i_dat, dat in enumerate(d_extracted[name]):
+1241 new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1242 for name in new_cov_names:
+1243 new_grad[name] = 0
+1244 for i_dat, dat in enumerate(g_extracted[name]):
+1245 new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1246 else:
+1247 for j_obs, obs in np.ndenumerate(data):
+1248 for name in obs.names:
+1249 if name in obs.cov_names:
+1250 new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
+1251 else:
+1252 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])
+1253
+1254 new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
+1255
+1256 if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
+1257 raise Exception('The same name has been used for deltas and covobs!')
+1258 new_samples = []
+1259 new_means = []
+1260 new_idl = []
+1261 new_names_obs = []
+1262 for name in new_names:
+1263 if name not in new_covobs:
+1264 new_samples.append(new_deltas[name])
+1265 new_idl.append(new_idl_d[name])
+1266 new_means.append(new_r_values[name][i_val])
+1267 new_names_obs.append(name)
+1268 final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
+1269 for name in new_covobs:
+1270 final_result[i_val].names.append(name)
+1271 final_result[i_val]._covobs = new_covobs
+1272 final_result[i_val]._value = new_val
+1273 final_result[i_val].reweighted = reweighted
+1274
+1275 if multi == 0:
+1276 final_result = final_result.item()
+1277
+1278 return final_result
+1279
+1280
+1281def _reduce_deltas(deltas, idx_old, idx_new):
+1282 """Extract deltas defined on idx_old on all configs of idx_new.
+1283
+1284 Assumes, that idx_old and idx_new are correctly defined idl, i.e., they
+1285 are ordered in an ascending order.
+1286
+1287 Parameters
+1288 ----------
+1289 deltas : list
+1290 List of fluctuations
+1291 idx_old : list
+1292 List or range of configs on which the deltas are defined
+1293 idx_new : list
+1294 List of configs for which we want to extract the deltas.
+1295 Has to be a subset of idx_old.
+1296 """
+1297 if not len(deltas) == len(idx_old):
+1298 raise Exception('Length of deltas and idx_old have to be the same: %d != %d' % (len(deltas), len(idx_old)))
+1299 if type(idx_old) is range and type(idx_new) is range:
+1300 if idx_old == idx_new:
+1301 return deltas
+1302 # Use groupby to efficiently check whether all elements of idx_old and idx_new are identical
+1303 try:
+1304 g = groupby([idx_old, idx_new])
+1305 if next(g, True) and not next(g, False):
+1306 return deltas
+1307 except Exception:
+1308 pass
+1309 indices = np.intersect1d(idx_old, idx_new, assume_unique=True, return_indices=True)[1]
+1310 if len(indices) < len(idx_new):
+1311 raise Exception('Error in _reduce_deltas: Config of idx_new not in idx_old')
+1312 return np.array(deltas)[indices]
+1313
+1314
+1315def reweight(weight, obs, **kwargs):
+1316 """Reweight a list of observables.
+1317
+1318 Parameters
+1319 ----------
+1320 weight : Obs
+1321 Reweighting factor. An Observable that has to be defined on a superset of the
+1322 configurations in obs[i].idl for all i.
+1323 obs : list
+1324 list of Obs, e.g. [obs1, obs2, obs3].
+1325 all_configs : bool
+1326 if True, the reweighted observables are normalized by the average of
+1327 the reweighting factor on all configurations in weight.idl and not
+1328 on the configurations in obs[i].idl. Default False.
+1329 """
+1330 result = []
+1331 for i in range(len(obs)):
+1332 if len(obs[i].cov_names):
+1333 raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
+1334 if not set(obs[i].names).issubset(weight.names):
+1335 raise Exception('Error: Ensembles do not fit')
+1336 for name in obs[i].names:
+1337 if not set(obs[i].idl[name]).issubset(weight.idl[name]):
+1338 raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
+1339 new_samples = []
+1340 w_deltas = {}
+1341 for name in sorted(obs[i].names):
+1342 w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
+1343 new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
+1344 tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
+1345
+1346 if kwargs.get('all_configs'):
+1347 new_weight = weight
+1348 else:
+1349 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)])
+1350
+1351 result.append(tmp_obs / new_weight)
+1352 result[-1].reweighted = True
1353
-1354 Notes
-1355 -----
-1356 Keep in mind to only correlate primary observables which have not been reweighted
-1357 yet. The reweighting has to be applied after correlating the observables.
-1358 Currently only works if ensembles are identical (this is not strictly necessary).
-1359 """
-1360
-1361 if sorted(obs_a.names) != sorted(obs_b.names):
-1362 raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
-1363 if len(obs_a.cov_names) or len(obs_b.cov_names):
-1364 raise Exception('Error: Not possible to correlate Obs that contain covobs!')
-1365 for name in obs_a.names:
-1366 if obs_a.shape[name] != obs_b.shape[name]:
-1367 raise Exception('Shapes of ensemble', name, 'do not fit')
-1368 if obs_a.idl[name] != obs_b.idl[name]:
-1369 raise Exception('idl of ensemble', name, 'do not fit')
-1370
-1371 if obs_a.reweighted is True:
-1372 warnings.warn("The first observable is already reweighted.", RuntimeWarning)
-1373 if obs_b.reweighted is True:
-1374 warnings.warn("The second observable is already reweighted.", RuntimeWarning)
-1375
-1376 new_samples = []
-1377 new_idl = []
-1378 for name in sorted(obs_a.names):
-1379 new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
-1380 new_idl.append(obs_a.idl[name])
-1381
-1382 o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
-1383 o.reweighted = obs_a.reweighted or obs_b.reweighted
-1384 return o
-1385
-1386
-1387def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
-1388 r'''Calculates the error covariance matrix of a set of observables.
-1389
-1390 WARNING: This function should be used with care, especially for observables with support on multiple
-1391 ensembles with differing autocorrelations. See the notes below for details.
-1392
-1393 The gamma method has to be applied first to all observables.
+1354 return result
+1355
+1356
+1357def correlate(obs_a, obs_b):
+1358 """Correlate two observables.
+1359
+1360 Parameters
+1361 ----------
+1362 obs_a : Obs
+1363 First observable
+1364 obs_b : Obs
+1365 Second observable
+1366
+1367 Notes
+1368 -----
+1369 Keep in mind to only correlate primary observables which have not been reweighted
+1370 yet. The reweighting has to be applied after correlating the observables.
+1371 Currently only works if ensembles are identical (this is not strictly necessary).
+1372 """
+1373
+1374 if sorted(obs_a.names) != sorted(obs_b.names):
+1375 raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
+1376 if len(obs_a.cov_names) or len(obs_b.cov_names):
+1377 raise Exception('Error: Not possible to correlate Obs that contain covobs!')
+1378 for name in obs_a.names:
+1379 if obs_a.shape[name] != obs_b.shape[name]:
+1380 raise Exception('Shapes of ensemble', name, 'do not fit')
+1381 if obs_a.idl[name] != obs_b.idl[name]:
+1382 raise Exception('idl of ensemble', name, 'do not fit')
+1383
+1384 if obs_a.reweighted is True:
+1385 warnings.warn("The first observable is already reweighted.", RuntimeWarning)
+1386 if obs_b.reweighted is True:
+1387 warnings.warn("The second observable is already reweighted.", RuntimeWarning)
+1388
+1389 new_samples = []
+1390 new_idl = []
+1391 for name in sorted(obs_a.names):
+1392 new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
+1393 new_idl.append(obs_a.idl[name])
1394
-1395 Parameters
-1396 ----------
-1397 obs : list or numpy.ndarray
-1398 List or one dimensional array of Obs
-1399 visualize : bool
-1400 If True plots the corresponding normalized correlation matrix (default False).
-1401 correlation : bool
-1402 If True the correlation matrix instead of the error covariance matrix is returned (default False).
-1403 smooth : None or int
-1404 If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
-1405 smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
-1406 largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
-1407 small ones.
-1408
-1409 Notes
-1410 -----
-1411 The error covariance is defined such that it agrees with the squared standard error for two identical observables
-1412 $$\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$$
-1413 in the absence of autocorrelation.
-1414 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
-1415 $$\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.
-1416 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.
-1417 $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
-1418 This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
-1419 '''
-1420
-1421 length = len(obs)
-1422
-1423 max_samples = np.max([o.N for o in obs])
-1424 if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
-1425 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)
-1426
-1427 cov = np.zeros((length, length))
-1428 for i in range(length):
-1429 for j in range(i, length):
-1430 cov[i, j] = _covariance_element(obs[i], obs[j])
-1431 cov = cov + cov.T - np.diag(np.diag(cov))
-1432
-1433 corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
-1434
-1435 if isinstance(smooth, int):
-1436 corr = _smooth_eigenvalues(corr, smooth)
-1437
-1438 if visualize:
-1439 plt.matshow(corr, vmin=-1, vmax=1)
-1440 plt.set_cmap('RdBu')
-1441 plt.colorbar()
-1442 plt.draw()
-1443
-1444 if correlation is True:
-1445 return corr
-1446
-1447 errors = [o.dvalue for o in obs]
-1448 cov = np.diag(errors) @ corr @ np.diag(errors)
-1449
-1450 eigenvalues = np.linalg.eigh(cov)[0]
-1451 if not np.all(eigenvalues >= 0):
-1452 warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
-1453
-1454 return cov
-1455
+1395 o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
+1396 o.reweighted = obs_a.reweighted or obs_b.reweighted
+1397 return o
+1398
+1399
+1400def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
+1401 r'''Calculates the error covariance matrix of a set of observables.
+1402
+1403 WARNING: This function should be used with care, especially for observables with support on multiple
+1404 ensembles with differing autocorrelations. See the notes below for details.
+1405
+1406 The gamma method has to be applied first to all observables.
+1407
+1408 Parameters
+1409 ----------
+1410 obs : list or numpy.ndarray
+1411 List or one dimensional array of Obs
+1412 visualize : bool
+1413 If True plots the corresponding normalized correlation matrix (default False).
+1414 correlation : bool
+1415 If True the correlation matrix instead of the error covariance matrix is returned (default False).
+1416 smooth : None or int
+1417 If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
+1418 smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
+1419 largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
+1420 small ones.
+1421
+1422 Notes
+1423 -----
+1424 The error covariance is defined such that it agrees with the squared standard error for two identical observables
+1425 $$\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$$
+1426 in the absence of autocorrelation.
+1427 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
+1428 $$\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.
+1429 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.
+1430 $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
+1431 This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
+1432 '''
+1433
+1434 length = len(obs)
+1435
+1436 max_samples = np.max([o.N for o in obs])
+1437 if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
+1438 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)
+1439
+1440 cov = np.zeros((length, length))
+1441 for i in range(length):
+1442 for j in range(i, length):
+1443 cov[i, j] = _covariance_element(obs[i], obs[j])
+1444 cov = cov + cov.T - np.diag(np.diag(cov))
+1445
+1446 corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+1447
+1448 if isinstance(smooth, int):
+1449 corr = _smooth_eigenvalues(corr, smooth)
+1450
+1451 if visualize:
+1452 plt.matshow(corr, vmin=-1, vmax=1)
+1453 plt.set_cmap('RdBu')
+1454 plt.colorbar()
+1455 plt.draw()
1456
-1457def _smooth_eigenvalues(corr, E):
-1458 """Eigenvalue smoothing as described in hep-lat/9412087
+1457 if correlation is True:
+1458 return corr
1459
-1460 corr : np.ndarray
-1461 correlation matrix
-1462 E : integer
-1463 Number of eigenvalues to be left substantially unchanged
-1464 """
-1465 if not (2 < E < corr.shape[0] - 1):
-1466 raise Exception(f"'E' has to be between 2 and the dimension of the correlation matrix minus 1 ({corr.shape[0] - 1}).")
-1467 vals, vec = np.linalg.eigh(corr)
-1468 lambda_min = np.mean(vals[:-E])
-1469 vals[vals < lambda_min] = lambda_min
-1470 vals /= np.mean(vals)
-1471 return vec @ np.diag(vals) @ vec.T
+1460 errors = [o.dvalue for o in obs]
+1461 cov = np.diag(errors) @ corr @ np.diag(errors)
+1462
+1463 eigenvalues = np.linalg.eigh(cov)[0]
+1464 if not np.all(eigenvalues >= 0):
+1465 warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
+1466
+1467 return cov
+1468
+1469
+1470def _smooth_eigenvalues(corr, E):
+1471 """Eigenvalue smoothing as described in hep-lat/9412087
1472
-1473
-1474def _covariance_element(obs1, obs2):
-1475 """Estimates the covariance of two Obs objects, neglecting autocorrelations."""
-1476
-1477 def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx):
-1478 deltas1 = _reduce_deltas(deltas1, idx1, new_idx)
-1479 deltas2 = _reduce_deltas(deltas2, idx2, new_idx)
-1480 return np.sum(deltas1 * deltas2)
-1481
-1482 if set(obs1.names).isdisjoint(set(obs2.names)):
-1483 return 0.0
-1484
-1485 if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'):
-1486 raise Exception('The gamma method has to be applied to both Obs first.')
-1487
-1488 dvalue = 0.0
+1473 corr : np.ndarray
+1474 correlation matrix
+1475 E : integer
+1476 Number of eigenvalues to be left substantially unchanged
+1477 """
+1478 if not (2 < E < corr.shape[0] - 1):
+1479 raise Exception(f"'E' has to be between 2 and the dimension of the correlation matrix minus 1 ({corr.shape[0] - 1}).")
+1480 vals, vec = np.linalg.eigh(corr)
+1481 lambda_min = np.mean(vals[:-E])
+1482 vals[vals < lambda_min] = lambda_min
+1483 vals /= np.mean(vals)
+1484 return vec @ np.diag(vals) @ vec.T
+1485
+1486
+1487def _covariance_element(obs1, obs2):
+1488 """Estimates the covariance of two Obs objects, neglecting autocorrelations."""
1489
-1490 for e_name in obs1.mc_names:
-1491
-1492 if e_name not in obs2.mc_names:
-1493 continue
+1490 def calc_gamma(deltas1, deltas2, idx1, idx2, new_idx):
+1491 deltas1 = _reduce_deltas(deltas1, idx1, new_idx)
+1492 deltas2 = _reduce_deltas(deltas2, idx2, new_idx)
+1493 return np.sum(deltas1 * deltas2)
1494
-1495 idl_d = {}
-1496 for r_name in obs1.e_content[e_name]:
-1497 if r_name not in obs2.e_content[e_name]:
-1498 continue
-1499 idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]])
+1495 if set(obs1.names).isdisjoint(set(obs2.names)):
+1496 return 0.0
+1497
+1498 if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'):
+1499 raise Exception('The gamma method has to be applied to both Obs first.')
1500
-1501 gamma = 0.0
+1501 dvalue = 0.0
1502
-1503 for r_name in obs1.e_content[e_name]:
-1504 if r_name not in obs2.e_content[e_name]:
-1505 continue
-1506 if len(idl_d[r_name]) == 0:
-1507 continue
-1508 gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name])
-1509
-1510 if gamma == 0.0:
-1511 continue
-1512
-1513 gamma_div = 0.0
-1514 for r_name in obs1.e_content[e_name]:
-1515 if r_name not in obs2.e_content[e_name]:
-1516 continue
-1517 if len(idl_d[r_name]) == 0:
+1503 for e_name in obs1.mc_names:
+1504
+1505 if e_name not in obs2.mc_names:
+1506 continue
+1507
+1508 idl_d = {}
+1509 for r_name in obs1.e_content[e_name]:
+1510 if r_name not in obs2.e_content[e_name]:
+1511 continue
+1512 idl_d[r_name] = _intersection_idx([obs1.idl[r_name], obs2.idl[r_name]])
+1513
+1514 gamma = 0.0
+1515
+1516 for r_name in obs1.e_content[e_name]:
+1517 if r_name not in obs2.e_content[e_name]:
1518 continue
-1519 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]))
-1520 gamma /= gamma_div
-1521
-1522 dvalue += gamma
-1523
-1524 for e_name in obs1.cov_names:
+1519 if len(idl_d[r_name]) == 0:
+1520 continue
+1521 gamma += calc_gamma(obs1.deltas[r_name], obs2.deltas[r_name], obs1.idl[r_name], obs2.idl[r_name], idl_d[r_name])
+1522
+1523 if gamma == 0.0:
+1524 continue
1525
-1526 if e_name not in obs2.cov_names:
-1527 continue
-1528
-1529 dvalue += float(np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad)))
-1530
-1531 return dvalue
-1532
-1533
-1534def import_jackknife(jacks, name, idl=None):
-1535 """Imports jackknife samples and returns an Obs
+1526 gamma_div = 0.0
+1527 for r_name in obs1.e_content[e_name]:
+1528 if r_name not in obs2.e_content[e_name]:
+1529 continue
+1530 if len(idl_d[r_name]) == 0:
+1531 continue
+1532 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]))
+1533 gamma /= gamma_div
+1534
+1535 dvalue += gamma
1536
-1537 Parameters
-1538 ----------
-1539 jacks : numpy.ndarray
-1540 numpy array containing the mean value as zeroth entry and
-1541 the N jackknife samples as first to Nth entry.
-1542 name : str
-1543 name of the ensemble the samples are defined on.
-1544 """
-1545 length = len(jacks) - 1
-1546 prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
-1547 samples = jacks[1:] @ prj
-1548 mean = np.mean(samples)
-1549 new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
-1550 new_obs._value = jacks[0]
-1551 return new_obs
-1552
-1553
-1554def merge_obs(list_of_obs):
-1555 """Combine all observables in list_of_obs into one new observable
-1556
-1557 Parameters
-1558 ----------
-1559 list_of_obs : list
-1560 list of the Obs object to be combined
-1561
-1562 Notes
-1563 -----
-1564 It is not possible to combine obs which are based on the same replicum
-1565 """
-1566 replist = [item for obs in list_of_obs for item in obs.names]
-1567 if (len(replist) == len(set(replist))) is False:
-1568 raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
-1569 if any([len(o.cov_names) for o in list_of_obs]):
-1570 raise Exception('Not possible to merge data that contains covobs!')
-1571 new_dict = {}
-1572 idl_dict = {}
-1573 for o in list_of_obs:
-1574 new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
-1575 for key in set(o.deltas) | set(o.r_values)})
-1576 idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
-1577
-1578 names = sorted(new_dict.keys())
-1579 o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
-1580 o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
-1581 return o
-1582
-1583
-1584def cov_Obs(means, cov, name, grad=None):
-1585 """Create an Obs based on mean(s) and a covariance matrix
-1586
-1587 Parameters
-1588 ----------
-1589 mean : list of floats or float
-1590 N mean value(s) of the new Obs
-1591 cov : list or array
-1592 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
-1593 name : str
-1594 identifier for the covariance matrix
-1595 grad : list or array
-1596 Gradient of the Covobs wrt. the means belonging to cov.
-1597 """
-1598
-1599 def covobs_to_obs(co):
-1600 """Make an Obs out of a Covobs
-1601
-1602 Parameters
-1603 ----------
-1604 co : Covobs
-1605 Covobs to be embedded into the Obs
-1606 """
-1607 o = Obs([], [], means=[])
-1608 o._value = co.value
-1609 o.names.append(co.name)
-1610 o._covobs[co.name] = co
-1611 o._dvalue = np.sqrt(co.errsq())
-1612 return o
-1613
-1614 ol = []
-1615 if isinstance(means, (float, int)):
-1616 means = [means]
-1617
-1618 for i in range(len(means)):
-1619 ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
-1620 if ol[0].covobs[name].N != len(means):
-1621 raise Exception('You have to provide %d mean values!' % (ol[0].N))
-1622 if len(ol) == 1:
-1623 return ol[0]
-1624 return ol
-1625
+1537 for e_name in obs1.cov_names:
+1538
+1539 if e_name not in obs2.cov_names:
+1540 continue
+1541
+1542 dvalue += float(np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad)))
+1543
+1544 return dvalue
+1545
+1546
+1547def import_jackknife(jacks, name, idl=None):
+1548 """Imports jackknife samples and returns an Obs
+1549
+1550 Parameters
+1551 ----------
+1552 jacks : numpy.ndarray
+1553 numpy array containing the mean value as zeroth entry and
+1554 the N jackknife samples as first to Nth entry.
+1555 name : str
+1556 name of the ensemble the samples are defined on.
+1557 """
+1558 length = len(jacks) - 1
+1559 prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
+1560 samples = jacks[1:] @ prj
+1561 mean = np.mean(samples)
+1562 new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
+1563 new_obs._value = jacks[0]
+1564 return new_obs
+1565
+1566
+1567def merge_obs(list_of_obs):
+1568 """Combine all observables in list_of_obs into one new observable
+1569
+1570 Parameters
+1571 ----------
+1572 list_of_obs : list
+1573 list of the Obs object to be combined
+1574
+1575 Notes
+1576 -----
+1577 It is not possible to combine obs which are based on the same replicum
+1578 """
+1579 replist = [item for obs in list_of_obs for item in obs.names]
+1580 if (len(replist) == len(set(replist))) is False:
+1581 raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
+1582 if any([len(o.cov_names) for o in list_of_obs]):
+1583 raise Exception('Not possible to merge data that contains covobs!')
+1584 new_dict = {}
+1585 idl_dict = {}
+1586 for o in list_of_obs:
+1587 new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
+1588 for key in set(o.deltas) | set(o.r_values)})
+1589 idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
+1590
+1591 names = sorted(new_dict.keys())
+1592 o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
+1593 o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
+1594 return o
+1595
+1596
+1597def cov_Obs(means, cov, name, grad=None):
+1598 """Create an Obs based on mean(s) and a covariance matrix
+1599
+1600 Parameters
+1601 ----------
+1602 mean : list of floats or float
+1603 N mean value(s) of the new Obs
+1604 cov : list or array
+1605 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
+1606 name : str
+1607 identifier for the covariance matrix
+1608 grad : list or array
+1609 Gradient of the Covobs wrt. the means belonging to cov.
+1610 """
+1611
+1612 def covobs_to_obs(co):
+1613 """Make an Obs out of a Covobs
+1614
+1615 Parameters
+1616 ----------
+1617 co : Covobs
+1618 Covobs to be embedded into the Obs
+1619 """
+1620 o = Obs([], [], means=[])
+1621 o._value = co.value
+1622 o.names.append(co.name)
+1623 o._covobs[co.name] = co
+1624 o._dvalue = np.sqrt(co.errsq())
+1625 return o
1626
-1627def _determine_gap(o, e_content, e_name):
-1628 gaps = []
-1629 for r_name in e_content[e_name]:
-1630 if isinstance(o.idl[r_name], range):
-1631 gaps.append(o.idl[r_name].step)
-1632 else:
-1633 gaps.append(np.min(np.diff(o.idl[r_name])))
-1634
-1635 gap = min(gaps)
-1636 if not np.all([gi % gap == 0 for gi in gaps]):
-1637 raise Exception(f"Replica for ensemble {e_name} do not have a common spacing.", gaps)
+1627 ol = []
+1628 if isinstance(means, (float, int)):
+1629 means = [means]
+1630
+1631 for i in range(len(means)):
+1632 ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
+1633 if ol[0].covobs[name].N != len(means):
+1634 raise Exception('You have to provide %d mean values!' % (ol[0].N))
+1635 if len(ol) == 1:
+1636 return ol[0]
+1637 return ol
1638
-1639 return gap
+1639
+1640def _determine_gap(o, e_content, e_name):
+1641 gaps = []
+1642 for r_name in e_content[e_name]:
+1643 if isinstance(o.idl[r_name], range):
+1644 gaps.append(o.idl[r_name].step)
+1645 else:
+1646 gaps.append(np.min(np.diff(o.idl[r_name])))
+1647
+1648 gap = min(gaps)
+1649 if not np.all([gi % gap == 0 for gi in gaps]):
+1650 raise Exception(f"Replica for ensemble {e_name} do not have a common spacing.", gaps)
+1651
+1652 return gap
@@ -2538,167 +2551,176 @@
696 def __str__(self):
697 return _format_uncertainty(self.value, self._dvalue)
698
-699 def __hash__(self):
-700 hash_tuple = (np.array([self.value]).astype(np.float32).data.tobytes(),)
-701 hash_tuple += tuple([o.astype(np.float32).data.tobytes() for o in self.deltas.values()])
-702 hash_tuple += tuple([np.array([o.errsq()]).astype(np.float32).data.tobytes() for o in self.covobs.values()])
-703 hash_tuple += tuple([o.encode() for o in self.names])
-704 m = hashlib.md5()
-705 [m.update(o) for o in hash_tuple]
-706 return int(m.hexdigest(), 16) & 0xFFFFFFFF
+699 def __format__(self, format_type):
+700 my_str = _format_uncertainty(self.value, self._dvalue,
+701 significance=int(float(format_type.replace("+", "").replace("-", ""))))
+702 for char in ["+", " "]:
+703 if format_type.startswith(char):
+704 if my_str[0] != "-":
+705 my_str = char + my_str
+706 return my_str
707
-708 # Overload comparisons
-709 def __lt__(self, other):
-710 return self.value < other
-711
-712 def __le__(self, other):
-713 return self.value <= other
-714
-715 def __gt__(self, other):
-716 return self.value > other
-717
-718 def __ge__(self, other):
-719 return self.value >= other
+708 def __hash__(self):
+709 hash_tuple = (np.array([self.value]).astype(np.float32).data.tobytes(),)
+710 hash_tuple += tuple([o.astype(np.float32).data.tobytes() for o in self.deltas.values()])
+711 hash_tuple += tuple([np.array([o.errsq()]).astype(np.float32).data.tobytes() for o in self.covobs.values()])
+712 hash_tuple += tuple([o.encode() for o in self.names])
+713 m = hashlib.md5()
+714 [m.update(o) for o in hash_tuple]
+715 return int(m.hexdigest(), 16) & 0xFFFFFFFF
+716
+717 # Overload comparisons
+718 def __lt__(self, other):
+719 return self.value < other
720
-721 def __eq__(self, other):
-722 return (self - other).is_zero()
+721 def __le__(self, other):
+722 return self.value <= other
723
-724 def __ne__(self, other):
-725 return not (self - other).is_zero()
+724 def __gt__(self, other):
+725 return self.value > other
726
-727 # Overload math operations
-728 def __add__(self, y):
-729 if isinstance(y, Obs):
-730 return derived_observable(lambda x, **kwargs: x[0] + x[1], [self, y], man_grad=[1, 1])
-731 else:
-732 if isinstance(y, np.ndarray):
-733 return np.array([self + o for o in y])
-734 elif y.__class__.__name__ in ['Corr', 'CObs']:
-735 return NotImplemented
-736 else:
-737 return derived_observable(lambda x, **kwargs: x[0] + y, [self], man_grad=[1])
-738
-739 def __radd__(self, y):
-740 return self + y
-741
-742 def __mul__(self, y):
-743 if isinstance(y, Obs):
-744 return derived_observable(lambda x, **kwargs: x[0] * x[1], [self, y], man_grad=[y.value, self.value])
-745 else:
-746 if isinstance(y, np.ndarray):
-747 return np.array([self * o for o in y])
-748 elif isinstance(y, complex):
-749 return CObs(self * y.real, self * y.imag)
-750 elif y.__class__.__name__ in ['Corr', 'CObs']:
-751 return NotImplemented
-752 else:
-753 return derived_observable(lambda x, **kwargs: x[0] * y, [self], man_grad=[y])
-754
-755 def __rmul__(self, y):
-756 return self * y
-757
-758 def __sub__(self, y):
-759 if isinstance(y, Obs):
-760 return derived_observable(lambda x, **kwargs: x[0] - x[1], [self, y], man_grad=[1, -1])
-761 else:
-762 if isinstance(y, np.ndarray):
-763 return np.array([self - o for o in y])
-764 elif y.__class__.__name__ in ['Corr', 'CObs']:
-765 return NotImplemented
-766 else:
-767 return derived_observable(lambda x, **kwargs: x[0] - y, [self], man_grad=[1])
-768
-769 def __rsub__(self, y):
-770 return -1 * (self - y)
-771
-772 def __pos__(self):
-773 return self
-774
-775 def __neg__(self):
-776 return -1 * self
+727 def __ge__(self, other):
+728 return self.value >= other
+729
+730 def __eq__(self, other):
+731 return (self - other).is_zero()
+732
+733 def __ne__(self, other):
+734 return not (self - other).is_zero()
+735
+736 # Overload math operations
+737 def __add__(self, y):
+738 if isinstance(y, Obs):
+739 return derived_observable(lambda x, **kwargs: x[0] + x[1], [self, y], man_grad=[1, 1])
+740 else:
+741 if isinstance(y, np.ndarray):
+742 return np.array([self + o for o in y])
+743 elif y.__class__.__name__ in ['Corr', 'CObs']:
+744 return NotImplemented
+745 else:
+746 return derived_observable(lambda x, **kwargs: x[0] + y, [self], man_grad=[1])
+747
+748 def __radd__(self, y):
+749 return self + y
+750
+751 def __mul__(self, y):
+752 if isinstance(y, Obs):
+753 return derived_observable(lambda x, **kwargs: x[0] * x[1], [self, y], man_grad=[y.value, self.value])
+754 else:
+755 if isinstance(y, np.ndarray):
+756 return np.array([self * o for o in y])
+757 elif isinstance(y, complex):
+758 return CObs(self * y.real, self * y.imag)
+759 elif y.__class__.__name__ in ['Corr', 'CObs']:
+760 return NotImplemented
+761 else:
+762 return derived_observable(lambda x, **kwargs: x[0] * y, [self], man_grad=[y])
+763
+764 def __rmul__(self, y):
+765 return self * y
+766
+767 def __sub__(self, y):
+768 if isinstance(y, Obs):
+769 return derived_observable(lambda x, **kwargs: x[0] - x[1], [self, y], man_grad=[1, -1])
+770 else:
+771 if isinstance(y, np.ndarray):
+772 return np.array([self - o for o in y])
+773 elif y.__class__.__name__ in ['Corr', 'CObs']:
+774 return NotImplemented
+775 else:
+776 return derived_observable(lambda x, **kwargs: x[0] - y, [self], man_grad=[1])
777
-778 def __truediv__(self, y):
-779 if isinstance(y, Obs):
-780 return derived_observable(lambda x, **kwargs: x[0] / x[1], [self, y], man_grad=[1 / y.value, - self.value / y.value ** 2])
-781 else:
-782 if isinstance(y, np.ndarray):
-783 return np.array([self / o for o in y])
-784 elif y.__class__.__name__ in ['Corr', 'CObs']:
-785 return NotImplemented
-786 else:
-787 return derived_observable(lambda x, **kwargs: x[0] / y, [self], man_grad=[1 / y])
-788
-789 def __rtruediv__(self, y):
-790 if isinstance(y, Obs):
-791 return derived_observable(lambda x, **kwargs: x[0] / x[1], [y, self], man_grad=[1 / self.value, - y.value / self.value ** 2])
-792 else:
-793 if isinstance(y, np.ndarray):
-794 return np.array([o / self for o in y])
-795 elif y.__class__.__name__ in ['Corr', 'CObs']:
-796 return NotImplemented
-797 else:
-798 return derived_observable(lambda x, **kwargs: y / x[0], [self], man_grad=[-y / self.value ** 2])
-799
-800 def __pow__(self, y):
-801 if isinstance(y, Obs):
-802 return derived_observable(lambda x: x[0] ** x[1], [self, y])
-803 else:
-804 return derived_observable(lambda x: x[0] ** y, [self])
-805
-806 def __rpow__(self, y):
-807 if isinstance(y, Obs):
-808 return derived_observable(lambda x: x[0] ** x[1], [y, self])
-809 else:
-810 return derived_observable(lambda x: y ** x[0], [self])
-811
-812 def __abs__(self):
-813 return derived_observable(lambda x: anp.abs(x[0]), [self])
+778 def __rsub__(self, y):
+779 return -1 * (self - y)
+780
+781 def __pos__(self):
+782 return self
+783
+784 def __neg__(self):
+785 return -1 * self
+786
+787 def __truediv__(self, y):
+788 if isinstance(y, Obs):
+789 return derived_observable(lambda x, **kwargs: x[0] / x[1], [self, y], man_grad=[1 / y.value, - self.value / y.value ** 2])
+790 else:
+791 if isinstance(y, np.ndarray):
+792 return np.array([self / o for o in y])
+793 elif y.__class__.__name__ in ['Corr', 'CObs']:
+794 return NotImplemented
+795 else:
+796 return derived_observable(lambda x, **kwargs: x[0] / y, [self], man_grad=[1 / y])
+797
+798 def __rtruediv__(self, y):
+799 if isinstance(y, Obs):
+800 return derived_observable(lambda x, **kwargs: x[0] / x[1], [y, self], man_grad=[1 / self.value, - y.value / self.value ** 2])
+801 else:
+802 if isinstance(y, np.ndarray):
+803 return np.array([o / self for o in y])
+804 elif y.__class__.__name__ in ['Corr', 'CObs']:
+805 return NotImplemented
+806 else:
+807 return derived_observable(lambda x, **kwargs: y / x[0], [self], man_grad=[-y / self.value ** 2])
+808
+809 def __pow__(self, y):
+810 if isinstance(y, Obs):
+811 return derived_observable(lambda x: x[0] ** x[1], [self, y])
+812 else:
+813 return derived_observable(lambda x: x[0] ** y, [self])
814
-815 # Overload numpy functions
-816 def sqrt(self):
-817 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
-818
-819 def log(self):
-820 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
-821
-822 def exp(self):
-823 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
-824
-825 def sin(self):
-826 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
+815 def __rpow__(self, y):
+816 if isinstance(y, Obs):
+817 return derived_observable(lambda x: x[0] ** x[1], [y, self])
+818 else:
+819 return derived_observable(lambda x: y ** x[0], [self])
+820
+821 def __abs__(self):
+822 return derived_observable(lambda x: anp.abs(x[0]), [self])
+823
+824 # Overload numpy functions
+825 def sqrt(self):
+826 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
827
-828 def cos(self):
-829 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
+828 def log(self):
+829 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
830
-831 def tan(self):
-832 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
+831 def exp(self):
+832 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
833
-834 def arcsin(self):
-835 return derived_observable(lambda x: anp.arcsin(x[0]), [self])
+834 def sin(self):
+835 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
836
-837 def arccos(self):
-838 return derived_observable(lambda x: anp.arccos(x[0]), [self])
+837 def cos(self):
+838 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
839
-840 def arctan(self):
-841 return derived_observable(lambda x: anp.arctan(x[0]), [self])
+840 def tan(self):
+841 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
842
-843 def sinh(self):
-844 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
+843 def arcsin(self):
+844 return derived_observable(lambda x: anp.arcsin(x[0]), [self])
845
-846 def cosh(self):
-847 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
+846 def arccos(self):
+847 return derived_observable(lambda x: anp.arccos(x[0]), [self])
848
-849 def tanh(self):
-850 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
+849 def arctan(self):
+850 return derived_observable(lambda x: anp.arctan(x[0]), [self])
851
-852 def arcsinh(self):
-853 return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
+852 def sinh(self):
+853 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
854
-855 def arccosh(self):
-856 return derived_observable(lambda x: anp.arccosh(x[0]), [self])
+855 def cosh(self):
+856 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
857
-858 def arctanh(self):
-859 return derived_observable(lambda x: anp.arctanh(x[0]), [self])
+858 def tanh(self):
+859 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
+860
+861 def arcsinh(self):
+862 return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
+863
+864 def arccosh(self):
+865 return derived_observable(lambda x: anp.arccosh(x[0]), [self])
+866
+867 def arctanh(self):
+868 return derived_observable(lambda x: anp.arctanh(x[0]), [self])
@@ -3856,8 +3878,8 @@ should agree with samples from a full jackknife analysis up to O(1/N).
- 816 def sqrt(self):
-817 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
+ 825 def sqrt(self):
+826 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
819 def log(self):
-820 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
+ 828 def log(self):
+829 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
822 def exp(self):
-823 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
+ 831 def exp(self):
+832 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
825 def sin(self):
-826 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
+ 834 def sin(self):
+835 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
828 def cos(self):
-829 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
+ 837 def cos(self):
+838 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
831 def tan(self):
-832 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
+ 840 def tan(self):
+841 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
834 def arcsin(self):
-835 return derived_observable(lambda x: anp.arcsin(x[0]), [self])
+ 843 def arcsin(self):
+844 return derived_observable(lambda x: anp.arcsin(x[0]), [self])
837 def arccos(self):
-838 return derived_observable(lambda x: anp.arccos(x[0]), [self])
+ 846 def arccos(self):
+847 return derived_observable(lambda x: anp.arccos(x[0]), [self])
840 def arctan(self):
-841 return derived_observable(lambda x: anp.arctan(x[0]), [self])
+ 849 def arctan(self):
+850 return derived_observable(lambda x: anp.arctan(x[0]), [self])
843 def sinh(self):
-844 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
+ 852 def sinh(self):
+853 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
846 def cosh(self):
-847 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
+ 855 def cosh(self):
+856 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
849 def tanh(self):
-850 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
+ 858 def tanh(self):
+859 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
852 def arcsinh(self):
-853 return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
+ 861 def arcsinh(self):
+862 return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
855 def arccosh(self):
-856 return derived_observable(lambda x: anp.arccosh(x[0]), [self])
+ 864 def arccosh(self):
+865 return derived_observable(lambda x: anp.arccosh(x[0]), [self])
858 def arctanh(self):
-859 return derived_observable(lambda x: anp.arctanh(x[0]), [self])
+ 867 def arctanh(self):
+868 return derived_observable(lambda x: anp.arctanh(x[0]), [self])
862class CObs:
-863 """Class for a complex valued observable."""
-864 __slots__ = ['_real', '_imag', 'tag']
-865
-866 def __init__(self, real, imag=0.0):
-867 self._real = real
-868 self._imag = imag
-869 self.tag = None
-870
-871 @property
-872 def real(self):
-873 return self._real
+ 871class CObs:
+872 """Class for a complex valued observable."""
+873 __slots__ = ['_real', '_imag', 'tag']
874
-875 @property
-876 def imag(self):
-877 return self._imag
-878
-879 def gamma_method(self, **kwargs):
-880 """Executes the gamma_method for the real and the imaginary part."""
-881 if isinstance(self.real, Obs):
-882 self.real.gamma_method(**kwargs)
-883 if isinstance(self.imag, Obs):
-884 self.imag.gamma_method(**kwargs)
-885
-886 def is_zero(self):
-887 """Checks whether both real and imaginary part are zero within machine precision."""
-888 return self.real == 0.0 and self.imag == 0.0
-889
-890 def conjugate(self):
-891 return CObs(self.real, -self.imag)
-892
-893 def __add__(self, other):
-894 if isinstance(other, np.ndarray):
-895 return other + self
-896 elif hasattr(other, 'real') and hasattr(other, 'imag'):
-897 return CObs(self.real + other.real,
-898 self.imag + other.imag)
-899 else:
-900 return CObs(self.real + other, self.imag)
+875 def __init__(self, real, imag=0.0):
+876 self._real = real
+877 self._imag = imag
+878 self.tag = None
+879
+880 @property
+881 def real(self):
+882 return self._real
+883
+884 @property
+885 def imag(self):
+886 return self._imag
+887
+888 def gamma_method(self, **kwargs):
+889 """Executes the gamma_method for the real and the imaginary part."""
+890 if isinstance(self.real, Obs):
+891 self.real.gamma_method(**kwargs)
+892 if isinstance(self.imag, Obs):
+893 self.imag.gamma_method(**kwargs)
+894
+895 def is_zero(self):
+896 """Checks whether both real and imaginary part are zero within machine precision."""
+897 return self.real == 0.0 and self.imag == 0.0
+898
+899 def conjugate(self):
+900 return CObs(self.real, -self.imag)
901
-902 def __radd__(self, y):
-903 return self + y
-904
-905 def __sub__(self, other):
-906 if isinstance(other, np.ndarray):
-907 return -1 * (other - self)
-908 elif hasattr(other, 'real') and hasattr(other, 'imag'):
-909 return CObs(self.real - other.real, self.imag - other.imag)
-910 else:
-911 return CObs(self.real - other, self.imag)
-912
-913 def __rsub__(self, other):
-914 return -1 * (self - other)
-915
-916 def __mul__(self, other):
-917 if isinstance(other, np.ndarray):
-918 return other * self
-919 elif hasattr(other, 'real') and hasattr(other, 'imag'):
-920 if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
-921 return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
-922 [self.real, other.real, self.imag, other.imag],
-923 man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
-924 derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
-925 [self.real, other.real, self.imag, other.imag],
-926 man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
-927 elif getattr(other, 'imag', 0) != 0:
-928 return CObs(self.real * other.real - self.imag * other.imag,
-929 self.imag * other.real + self.real * other.imag)
-930 else:
-931 return CObs(self.real * other.real, self.imag * other.real)
-932 else:
-933 return CObs(self.real * other, self.imag * other)
-934
-935 def __rmul__(self, other):
-936 return self * other
-937
-938 def __truediv__(self, other):
-939 if isinstance(other, np.ndarray):
-940 return 1 / (other / self)
-941 elif hasattr(other, 'real') and hasattr(other, 'imag'):
-942 r = other.real ** 2 + other.imag ** 2
-943 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
-944 else:
-945 return CObs(self.real / other, self.imag / other)
+902 def __add__(self, other):
+903 if isinstance(other, np.ndarray):
+904 return other + self
+905 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+906 return CObs(self.real + other.real,
+907 self.imag + other.imag)
+908 else:
+909 return CObs(self.real + other, self.imag)
+910
+911 def __radd__(self, y):
+912 return self + y
+913
+914 def __sub__(self, other):
+915 if isinstance(other, np.ndarray):
+916 return -1 * (other - self)
+917 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+918 return CObs(self.real - other.real, self.imag - other.imag)
+919 else:
+920 return CObs(self.real - other, self.imag)
+921
+922 def __rsub__(self, other):
+923 return -1 * (self - other)
+924
+925 def __mul__(self, other):
+926 if isinstance(other, np.ndarray):
+927 return other * self
+928 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+929 if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
+930 return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
+931 [self.real, other.real, self.imag, other.imag],
+932 man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
+933 derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
+934 [self.real, other.real, self.imag, other.imag],
+935 man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
+936 elif getattr(other, 'imag', 0) != 0:
+937 return CObs(self.real * other.real - self.imag * other.imag,
+938 self.imag * other.real + self.real * other.imag)
+939 else:
+940 return CObs(self.real * other.real, self.imag * other.real)
+941 else:
+942 return CObs(self.real * other, self.imag * other)
+943
+944 def __rmul__(self, other):
+945 return self * other
946
-947 def __rtruediv__(self, other):
-948 r = self.real ** 2 + self.imag ** 2
-949 if hasattr(other, 'real') and hasattr(other, 'imag'):
-950 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
-951 else:
-952 return CObs(self.real * other / r, -self.imag * other / r)
-953
-954 def __abs__(self):
-955 return np.sqrt(self.real**2 + self.imag**2)
-956
-957 def __pos__(self):
-958 return self
-959
-960 def __neg__(self):
-961 return -1 * self
+947 def __truediv__(self, other):
+948 if isinstance(other, np.ndarray):
+949 return 1 / (other / self)
+950 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+951 r = other.real ** 2 + other.imag ** 2
+952 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
+953 else:
+954 return CObs(self.real / other, self.imag / other)
+955
+956 def __rtruediv__(self, other):
+957 r = self.real ** 2 + self.imag ** 2
+958 if hasattr(other, 'real') and hasattr(other, 'imag'):
+959 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
+960 else:
+961 return CObs(self.real * other / r, -self.imag * other / r)
962
-963 def __eq__(self, other):
-964 return self.real == other.real and self.imag == other.imag
+963 def __abs__(self):
+964 return np.sqrt(self.real**2 + self.imag**2)
965
-966 def __str__(self):
-967 return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
+966 def __pos__(self):
+967 return self
968
-969 def __repr__(self):
-970 return 'CObs[' + str(self) + ']'
+969 def __neg__(self):
+970 return -1 * self
+971
+972 def __eq__(self, other):
+973 return self.real == other.real and self.imag == other.imag
+974
+975 def __str__(self):
+976 return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
+977
+978 def __repr__(self):
+979 return 'CObs[' + str(self) + ']'
866 def __init__(self, real, imag=0.0):
-867 self._real = real
-868 self._imag = imag
-869 self.tag = None
+ 875 def __init__(self, real, imag=0.0):
+876 self._real = real
+877 self._imag = imag
+878 self.tag = None
879 def gamma_method(self, **kwargs):
-880 """Executes the gamma_method for the real and the imaginary part."""
-881 if isinstance(self.real, Obs):
-882 self.real.gamma_method(**kwargs)
-883 if isinstance(self.imag, Obs):
-884 self.imag.gamma_method(**kwargs)
+ 888 def gamma_method(self, **kwargs):
+889 """Executes the gamma_method for the real and the imaginary part."""
+890 if isinstance(self.real, Obs):
+891 self.real.gamma_method(**kwargs)
+892 if isinstance(self.imag, Obs):
+893 self.imag.gamma_method(**kwargs)
886 def is_zero(self):
-887 """Checks whether both real and imaginary part are zero within machine precision."""
-888 return self.real == 0.0 and self.imag == 0.0
+ 895 def is_zero(self):
+896 """Checks whether both real and imaginary part are zero within machine precision."""
+897 return self.real == 0.0 and self.imag == 0.0
890 def conjugate(self):
-891 return CObs(self.real, -self.imag)
+ 899 def conjugate(self):
+900 return CObs(self.real, -self.imag)
1099def derived_observable(func, data, array_mode=False, **kwargs):
-1100 """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
-1101
-1102 Parameters
-1103 ----------
-1104 func : object
-1105 arbitrary function of the form func(data, **kwargs). For the
-1106 automatic differentiation to work, all numpy functions have to have
-1107 the autograd wrapper (use 'import autograd.numpy as anp').
-1108 data : list
-1109 list of Obs, e.g. [obs1, obs2, obs3].
-1110 num_grad : bool
-1111 if True, numerical derivatives are used instead of autograd
-1112 (default False). To control the numerical differentiation the
-1113 kwargs of numdifftools.step_generators.MaxStepGenerator
-1114 can be used.
-1115 man_grad : list
-1116 manually supply a list or an array which contains the jacobian
-1117 of func. Use cautiously, supplying the wrong derivative will
-1118 not be intercepted.
-1119
-1120 Notes
-1121 -----
-1122 For simple mathematical operations it can be practical to use anonymous
-1123 functions. For the ratio of two observables one can e.g. use
-1124
-1125 new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
-1126 """
-1127
-1128 data = np.asarray(data)
-1129 raveled_data = data.ravel()
-1130
-1131 # Workaround for matrix operations containing non Obs data
-1132 if not all(isinstance(x, Obs) for x in raveled_data):
-1133 for i in range(len(raveled_data)):
-1134 if isinstance(raveled_data[i], (int, float)):
-1135 raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
-1136
-1137 allcov = {}
-1138 for o in raveled_data:
-1139 for name in o.cov_names:
-1140 if name in allcov:
-1141 if not np.allclose(allcov[name], o.covobs[name].cov):
-1142 raise Exception('Inconsistent covariance matrices for %s!' % (name))
-1143 else:
-1144 allcov[name] = o.covobs[name].cov
-1145
-1146 n_obs = len(raveled_data)
-1147 new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
-1148 new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
-1149 new_sample_names = sorted(set(new_names) - set(new_cov_names))
-1150
-1151 reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
-1152
-1153 if data.ndim == 1:
-1154 values = np.array([o.value for o in data])
-1155 else:
-1156 values = np.vectorize(lambda x: x.value)(data)
-1157
-1158 new_values = func(values, **kwargs)
-1159
-1160 multi = int(isinstance(new_values, np.ndarray))
-1161
-1162 new_r_values = {}
-1163 new_idl_d = {}
-1164 for name in new_sample_names:
-1165 idl = []
-1166 tmp_values = np.zeros(n_obs)
-1167 for i, item in enumerate(raveled_data):
-1168 tmp_values[i] = item.r_values.get(name, item.value)
-1169 tmp_idl = item.idl.get(name)
-1170 if tmp_idl is not None:
-1171 idl.append(tmp_idl)
-1172 if multi > 0:
-1173 tmp_values = np.array(tmp_values).reshape(data.shape)
-1174 new_r_values[name] = func(tmp_values, **kwargs)
-1175 new_idl_d[name] = _merge_idx(idl)
-1176
-1177 if 'man_grad' in kwargs:
-1178 deriv = np.asarray(kwargs.get('man_grad'))
-1179 if new_values.shape + data.shape != deriv.shape:
-1180 raise Exception('Manual derivative does not have correct shape.')
-1181 elif kwargs.get('num_grad') is True:
-1182 if multi > 0:
-1183 raise Exception('Multi mode currently not supported for numerical derivative')
-1184 options = {
-1185 'base_step': 0.1,
-1186 'step_ratio': 2.5}
-1187 for key in options.keys():
-1188 kwarg = kwargs.get(key)
-1189 if kwarg is not None:
-1190 options[key] = kwarg
-1191 tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
-1192 if tmp_df.size == 1:
-1193 deriv = np.array([tmp_df.real])
-1194 else:
-1195 deriv = tmp_df.real
-1196 else:
-1197 deriv = jacobian(func)(values, **kwargs)
-1198
-1199 final_result = np.zeros(new_values.shape, dtype=object)
-1200
-1201 if array_mode is True:
-1202
-1203 class _Zero_grad():
-1204 def __init__(self, N):
-1205 self.grad = np.zeros((N, 1))
-1206
-1207 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]))
-1208 d_extracted = {}
-1209 g_extracted = {}
-1210 for name in new_sample_names:
-1211 d_extracted[name] = []
-1212 ens_length = len(new_idl_d[name])
-1213 for i_dat, dat in enumerate(data):
-1214 d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name]) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
-1215 for name in new_cov_names:
-1216 g_extracted[name] = []
-1217 zero_grad = _Zero_grad(new_covobs_lengths[name])
-1218 for i_dat, dat in enumerate(data):
-1219 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)))
-1220
-1221 for i_val, new_val in np.ndenumerate(new_values):
-1222 new_deltas = {}
-1223 new_grad = {}
-1224 if array_mode is True:
-1225 for name in new_sample_names:
-1226 ens_length = d_extracted[name][0].shape[-1]
-1227 new_deltas[name] = np.zeros(ens_length)
-1228 for i_dat, dat in enumerate(d_extracted[name]):
-1229 new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1230 for name in new_cov_names:
-1231 new_grad[name] = 0
-1232 for i_dat, dat in enumerate(g_extracted[name]):
-1233 new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
-1234 else:
-1235 for j_obs, obs in np.ndenumerate(data):
-1236 for name in obs.names:
-1237 if name in obs.cov_names:
-1238 new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
-1239 else:
-1240 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])
-1241
-1242 new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
-1243
-1244 if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
-1245 raise Exception('The same name has been used for deltas and covobs!')
-1246 new_samples = []
-1247 new_means = []
-1248 new_idl = []
-1249 new_names_obs = []
-1250 for name in new_names:
-1251 if name not in new_covobs:
-1252 new_samples.append(new_deltas[name])
-1253 new_idl.append(new_idl_d[name])
-1254 new_means.append(new_r_values[name][i_val])
-1255 new_names_obs.append(name)
-1256 final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
-1257 for name in new_covobs:
-1258 final_result[i_val].names.append(name)
-1259 final_result[i_val]._covobs = new_covobs
-1260 final_result[i_val]._value = new_val
-1261 final_result[i_val].reweighted = reweighted
-1262
-1263 if multi == 0:
-1264 final_result = final_result.item()
-1265
-1266 return final_result
+ 1112def derived_observable(func, data, array_mode=False, **kwargs):
+1113 """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
+1114
+1115 Parameters
+1116 ----------
+1117 func : object
+1118 arbitrary function of the form func(data, **kwargs). For the
+1119 automatic differentiation to work, all numpy functions have to have
+1120 the autograd wrapper (use 'import autograd.numpy as anp').
+1121 data : list
+1122 list of Obs, e.g. [obs1, obs2, obs3].
+1123 num_grad : bool
+1124 if True, numerical derivatives are used instead of autograd
+1125 (default False). To control the numerical differentiation the
+1126 kwargs of numdifftools.step_generators.MaxStepGenerator
+1127 can be used.
+1128 man_grad : list
+1129 manually supply a list or an array which contains the jacobian
+1130 of func. Use cautiously, supplying the wrong derivative will
+1131 not be intercepted.
+1132
+1133 Notes
+1134 -----
+1135 For simple mathematical operations it can be practical to use anonymous
+1136 functions. For the ratio of two observables one can e.g. use
+1137
+1138 new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
+1139 """
+1140
+1141 data = np.asarray(data)
+1142 raveled_data = data.ravel()
+1143
+1144 # Workaround for matrix operations containing non Obs data
+1145 if not all(isinstance(x, Obs) for x in raveled_data):
+1146 for i in range(len(raveled_data)):
+1147 if isinstance(raveled_data[i], (int, float)):
+1148 raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
+1149
+1150 allcov = {}
+1151 for o in raveled_data:
+1152 for name in o.cov_names:
+1153 if name in allcov:
+1154 if not np.allclose(allcov[name], o.covobs[name].cov):
+1155 raise Exception('Inconsistent covariance matrices for %s!' % (name))
+1156 else:
+1157 allcov[name] = o.covobs[name].cov
+1158
+1159 n_obs = len(raveled_data)
+1160 new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
+1161 new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
+1162 new_sample_names = sorted(set(new_names) - set(new_cov_names))
+1163
+1164 reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
+1165
+1166 if data.ndim == 1:
+1167 values = np.array([o.value for o in data])
+1168 else:
+1169 values = np.vectorize(lambda x: x.value)(data)
+1170
+1171 new_values = func(values, **kwargs)
+1172
+1173 multi = int(isinstance(new_values, np.ndarray))
+1174
+1175 new_r_values = {}
+1176 new_idl_d = {}
+1177 for name in new_sample_names:
+1178 idl = []
+1179 tmp_values = np.zeros(n_obs)
+1180 for i, item in enumerate(raveled_data):
+1181 tmp_values[i] = item.r_values.get(name, item.value)
+1182 tmp_idl = item.idl.get(name)
+1183 if tmp_idl is not None:
+1184 idl.append(tmp_idl)
+1185 if multi > 0:
+1186 tmp_values = np.array(tmp_values).reshape(data.shape)
+1187 new_r_values[name] = func(tmp_values, **kwargs)
+1188 new_idl_d[name] = _merge_idx(idl)
+1189
+1190 if 'man_grad' in kwargs:
+1191 deriv = np.asarray(kwargs.get('man_grad'))
+1192 if new_values.shape + data.shape != deriv.shape:
+1193 raise Exception('Manual derivative does not have correct shape.')
+1194 elif kwargs.get('num_grad') is True:
+1195 if multi > 0:
+1196 raise Exception('Multi mode currently not supported for numerical derivative')
+1197 options = {
+1198 'base_step': 0.1,
+1199 'step_ratio': 2.5}
+1200 for key in options.keys():
+1201 kwarg = kwargs.get(key)
+1202 if kwarg is not None:
+1203 options[key] = kwarg
+1204 tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
+1205 if tmp_df.size == 1:
+1206 deriv = np.array([tmp_df.real])
+1207 else:
+1208 deriv = tmp_df.real
+1209 else:
+1210 deriv = jacobian(func)(values, **kwargs)
+1211
+1212 final_result = np.zeros(new_values.shape, dtype=object)
+1213
+1214 if array_mode is True:
+1215
+1216 class _Zero_grad():
+1217 def __init__(self, N):
+1218 self.grad = np.zeros((N, 1))
+1219
+1220 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]))
+1221 d_extracted = {}
+1222 g_extracted = {}
+1223 for name in new_sample_names:
+1224 d_extracted[name] = []
+1225 ens_length = len(new_idl_d[name])
+1226 for i_dat, dat in enumerate(data):
+1227 d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name]) for o in dat.reshape(np.prod(dat.shape))]).reshape(dat.shape + (ens_length, )))
+1228 for name in new_cov_names:
+1229 g_extracted[name] = []
+1230 zero_grad = _Zero_grad(new_covobs_lengths[name])
+1231 for i_dat, dat in enumerate(data):
+1232 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)))
+1233
+1234 for i_val, new_val in np.ndenumerate(new_values):
+1235 new_deltas = {}
+1236 new_grad = {}
+1237 if array_mode is True:
+1238 for name in new_sample_names:
+1239 ens_length = d_extracted[name][0].shape[-1]
+1240 new_deltas[name] = np.zeros(ens_length)
+1241 for i_dat, dat in enumerate(d_extracted[name]):
+1242 new_deltas[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1243 for name in new_cov_names:
+1244 new_grad[name] = 0
+1245 for i_dat, dat in enumerate(g_extracted[name]):
+1246 new_grad[name] += np.tensordot(deriv[i_val + (i_dat, )], dat)
+1247 else:
+1248 for j_obs, obs in np.ndenumerate(data):
+1249 for name in obs.names:
+1250 if name in obs.cov_names:
+1251 new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
+1252 else:
+1253 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])
+1254
+1255 new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
+1256
+1257 if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
+1258 raise Exception('The same name has been used for deltas and covobs!')
+1259 new_samples = []
+1260 new_means = []
+1261 new_idl = []
+1262 new_names_obs = []
+1263 for name in new_names:
+1264 if name not in new_covobs:
+1265 new_samples.append(new_deltas[name])
+1266 new_idl.append(new_idl_d[name])
+1267 new_means.append(new_r_values[name][i_val])
+1268 new_names_obs.append(name)
+1269 final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
+1270 for name in new_covobs:
+1271 final_result[i_val].names.append(name)
+1272 final_result[i_val]._covobs = new_covobs
+1273 final_result[i_val]._value = new_val
+1274 final_result[i_val].reweighted = reweighted
+1275
+1276 if multi == 0:
+1277 final_result = final_result.item()
+1278
+1279 return final_result
1303def reweight(weight, obs, **kwargs):
-1304 """Reweight a list of observables.
-1305
-1306 Parameters
-1307 ----------
-1308 weight : Obs
-1309 Reweighting factor. An Observable that has to be defined on a superset of the
-1310 configurations in obs[i].idl for all i.
-1311 obs : list
-1312 list of Obs, e.g. [obs1, obs2, obs3].
-1313 all_configs : bool
-1314 if True, the reweighted observables are normalized by the average of
-1315 the reweighting factor on all configurations in weight.idl and not
-1316 on the configurations in obs[i].idl. Default False.
-1317 """
-1318 result = []
-1319 for i in range(len(obs)):
-1320 if len(obs[i].cov_names):
-1321 raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
-1322 if not set(obs[i].names).issubset(weight.names):
-1323 raise Exception('Error: Ensembles do not fit')
-1324 for name in obs[i].names:
-1325 if not set(obs[i].idl[name]).issubset(weight.idl[name]):
-1326 raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
-1327 new_samples = []
-1328 w_deltas = {}
-1329 for name in sorted(obs[i].names):
-1330 w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
-1331 new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
-1332 tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
-1333
-1334 if kwargs.get('all_configs'):
-1335 new_weight = weight
-1336 else:
-1337 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)])
-1338
-1339 result.append(tmp_obs / new_weight)
-1340 result[-1].reweighted = True
-1341
-1342 return result
+ 1316def reweight(weight, obs, **kwargs):
+1317 """Reweight a list of observables.
+1318
+1319 Parameters
+1320 ----------
+1321 weight : Obs
+1322 Reweighting factor. An Observable that has to be defined on a superset of the
+1323 configurations in obs[i].idl for all i.
+1324 obs : list
+1325 list of Obs, e.g. [obs1, obs2, obs3].
+1326 all_configs : bool
+1327 if True, the reweighted observables are normalized by the average of
+1328 the reweighting factor on all configurations in weight.idl and not
+1329 on the configurations in obs[i].idl. Default False.
+1330 """
+1331 result = []
+1332 for i in range(len(obs)):
+1333 if len(obs[i].cov_names):
+1334 raise Exception('Error: Not possible to reweight an Obs that contains covobs!')
+1335 if not set(obs[i].names).issubset(weight.names):
+1336 raise Exception('Error: Ensembles do not fit')
+1337 for name in obs[i].names:
+1338 if not set(obs[i].idl[name]).issubset(weight.idl[name]):
+1339 raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name))
+1340 new_samples = []
+1341 w_deltas = {}
+1342 for name in sorted(obs[i].names):
+1343 w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
+1344 new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
+1345 tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
+1346
+1347 if kwargs.get('all_configs'):
+1348 new_weight = weight
+1349 else:
+1350 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)])
+1351
+1352 result.append(tmp_obs / new_weight)
+1353 result[-1].reweighted = True
+1354
+1355 return result
1345def correlate(obs_a, obs_b):
-1346 """Correlate two observables.
-1347
-1348 Parameters
-1349 ----------
-1350 obs_a : Obs
-1351 First observable
-1352 obs_b : Obs
-1353 Second observable
-1354
-1355 Notes
-1356 -----
-1357 Keep in mind to only correlate primary observables which have not been reweighted
-1358 yet. The reweighting has to be applied after correlating the observables.
-1359 Currently only works if ensembles are identical (this is not strictly necessary).
-1360 """
-1361
-1362 if sorted(obs_a.names) != sorted(obs_b.names):
-1363 raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
-1364 if len(obs_a.cov_names) or len(obs_b.cov_names):
-1365 raise Exception('Error: Not possible to correlate Obs that contain covobs!')
-1366 for name in obs_a.names:
-1367 if obs_a.shape[name] != obs_b.shape[name]:
-1368 raise Exception('Shapes of ensemble', name, 'do not fit')
-1369 if obs_a.idl[name] != obs_b.idl[name]:
-1370 raise Exception('idl of ensemble', name, 'do not fit')
-1371
-1372 if obs_a.reweighted is True:
-1373 warnings.warn("The first observable is already reweighted.", RuntimeWarning)
-1374 if obs_b.reweighted is True:
-1375 warnings.warn("The second observable is already reweighted.", RuntimeWarning)
-1376
-1377 new_samples = []
-1378 new_idl = []
-1379 for name in sorted(obs_a.names):
-1380 new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
-1381 new_idl.append(obs_a.idl[name])
-1382
-1383 o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
-1384 o.reweighted = obs_a.reweighted or obs_b.reweighted
-1385 return o
+ 1358def correlate(obs_a, obs_b):
+1359 """Correlate two observables.
+1360
+1361 Parameters
+1362 ----------
+1363 obs_a : Obs
+1364 First observable
+1365 obs_b : Obs
+1366 Second observable
+1367
+1368 Notes
+1369 -----
+1370 Keep in mind to only correlate primary observables which have not been reweighted
+1371 yet. The reweighting has to be applied after correlating the observables.
+1372 Currently only works if ensembles are identical (this is not strictly necessary).
+1373 """
+1374
+1375 if sorted(obs_a.names) != sorted(obs_b.names):
+1376 raise Exception(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
+1377 if len(obs_a.cov_names) or len(obs_b.cov_names):
+1378 raise Exception('Error: Not possible to correlate Obs that contain covobs!')
+1379 for name in obs_a.names:
+1380 if obs_a.shape[name] != obs_b.shape[name]:
+1381 raise Exception('Shapes of ensemble', name, 'do not fit')
+1382 if obs_a.idl[name] != obs_b.idl[name]:
+1383 raise Exception('idl of ensemble', name, 'do not fit')
+1384
+1385 if obs_a.reweighted is True:
+1386 warnings.warn("The first observable is already reweighted.", RuntimeWarning)
+1387 if obs_b.reweighted is True:
+1388 warnings.warn("The second observable is already reweighted.", RuntimeWarning)
+1389
+1390 new_samples = []
+1391 new_idl = []
+1392 for name in sorted(obs_a.names):
+1393 new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
+1394 new_idl.append(obs_a.idl[name])
+1395
+1396 o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
+1397 o.reweighted = obs_a.reweighted or obs_b.reweighted
+1398 return o
1388def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
-1389 r'''Calculates the error covariance matrix of a set of observables.
-1390
-1391 WARNING: This function should be used with care, especially for observables with support on multiple
-1392 ensembles with differing autocorrelations. See the notes below for details.
-1393
-1394 The gamma method has to be applied first to all observables.
-1395
-1396 Parameters
-1397 ----------
-1398 obs : list or numpy.ndarray
-1399 List or one dimensional array of Obs
-1400 visualize : bool
-1401 If True plots the corresponding normalized correlation matrix (default False).
-1402 correlation : bool
-1403 If True the correlation matrix instead of the error covariance matrix is returned (default False).
-1404 smooth : None or int
-1405 If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
-1406 smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
-1407 largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
-1408 small ones.
-1409
-1410 Notes
-1411 -----
-1412 The error covariance is defined such that it agrees with the squared standard error for two identical observables
-1413 $$\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$$
-1414 in the absence of autocorrelation.
-1415 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
-1416 $$\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.
-1417 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.
-1418 $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
-1419 This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
-1420 '''
-1421
-1422 length = len(obs)
-1423
-1424 max_samples = np.max([o.N for o in obs])
-1425 if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
-1426 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)
-1427
-1428 cov = np.zeros((length, length))
-1429 for i in range(length):
-1430 for j in range(i, length):
-1431 cov[i, j] = _covariance_element(obs[i], obs[j])
-1432 cov = cov + cov.T - np.diag(np.diag(cov))
-1433
-1434 corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
-1435
-1436 if isinstance(smooth, int):
-1437 corr = _smooth_eigenvalues(corr, smooth)
-1438
-1439 if visualize:
-1440 plt.matshow(corr, vmin=-1, vmax=1)
-1441 plt.set_cmap('RdBu')
-1442 plt.colorbar()
-1443 plt.draw()
-1444
-1445 if correlation is True:
-1446 return corr
-1447
-1448 errors = [o.dvalue for o in obs]
-1449 cov = np.diag(errors) @ corr @ np.diag(errors)
-1450
-1451 eigenvalues = np.linalg.eigh(cov)[0]
-1452 if not np.all(eigenvalues >= 0):
-1453 warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
-1454
-1455 return cov
+ 1401def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
+1402 r'''Calculates the error covariance matrix of a set of observables.
+1403
+1404 WARNING: This function should be used with care, especially for observables with support on multiple
+1405 ensembles with differing autocorrelations. See the notes below for details.
+1406
+1407 The gamma method has to be applied first to all observables.
+1408
+1409 Parameters
+1410 ----------
+1411 obs : list or numpy.ndarray
+1412 List or one dimensional array of Obs
+1413 visualize : bool
+1414 If True plots the corresponding normalized correlation matrix (default False).
+1415 correlation : bool
+1416 If True the correlation matrix instead of the error covariance matrix is returned (default False).
+1417 smooth : None or int
+1418 If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
+1419 smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
+1420 largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
+1421 small ones.
+1422
+1423 Notes
+1424 -----
+1425 The error covariance is defined such that it agrees with the squared standard error for two identical observables
+1426 $$\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$$
+1427 in the absence of autocorrelation.
+1428 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
+1429 $$\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.
+1430 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.
+1431 $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
+1432 This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
+1433 '''
+1434
+1435 length = len(obs)
+1436
+1437 max_samples = np.max([o.N for o in obs])
+1438 if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
+1439 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)
+1440
+1441 cov = np.zeros((length, length))
+1442 for i in range(length):
+1443 for j in range(i, length):
+1444 cov[i, j] = _covariance_element(obs[i], obs[j])
+1445 cov = cov + cov.T - np.diag(np.diag(cov))
+1446
+1447 corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+1448
+1449 if isinstance(smooth, int):
+1450 corr = _smooth_eigenvalues(corr, smooth)
+1451
+1452 if visualize:
+1453 plt.matshow(corr, vmin=-1, vmax=1)
+1454 plt.set_cmap('RdBu')
+1455 plt.colorbar()
+1456 plt.draw()
+1457
+1458 if correlation is True:
+1459 return corr
+1460
+1461 errors = [o.dvalue for o in obs]
+1462 cov = np.diag(errors) @ corr @ np.diag(errors)
+1463
+1464 eigenvalues = np.linalg.eigh(cov)[0]
+1465 if not np.all(eigenvalues >= 0):
+1466 warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning)
+1467
+1468 return cov
1535def import_jackknife(jacks, name, idl=None):
-1536 """Imports jackknife samples and returns an Obs
-1537
-1538 Parameters
-1539 ----------
-1540 jacks : numpy.ndarray
-1541 numpy array containing the mean value as zeroth entry and
-1542 the N jackknife samples as first to Nth entry.
-1543 name : str
-1544 name of the ensemble the samples are defined on.
-1545 """
-1546 length = len(jacks) - 1
-1547 prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
-1548 samples = jacks[1:] @ prj
-1549 mean = np.mean(samples)
-1550 new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
-1551 new_obs._value = jacks[0]
-1552 return new_obs
+ 1548def import_jackknife(jacks, name, idl=None):
+1549 """Imports jackknife samples and returns an Obs
+1550
+1551 Parameters
+1552 ----------
+1553 jacks : numpy.ndarray
+1554 numpy array containing the mean value as zeroth entry and
+1555 the N jackknife samples as first to Nth entry.
+1556 name : str
+1557 name of the ensemble the samples are defined on.
+1558 """
+1559 length = len(jacks) - 1
+1560 prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
+1561 samples = jacks[1:] @ prj
+1562 mean = np.mean(samples)
+1563 new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
+1564 new_obs._value = jacks[0]
+1565 return new_obs
1555def merge_obs(list_of_obs):
-1556 """Combine all observables in list_of_obs into one new observable
-1557
-1558 Parameters
-1559 ----------
-1560 list_of_obs : list
-1561 list of the Obs object to be combined
-1562
-1563 Notes
-1564 -----
-1565 It is not possible to combine obs which are based on the same replicum
-1566 """
-1567 replist = [item for obs in list_of_obs for item in obs.names]
-1568 if (len(replist) == len(set(replist))) is False:
-1569 raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
-1570 if any([len(o.cov_names) for o in list_of_obs]):
-1571 raise Exception('Not possible to merge data that contains covobs!')
-1572 new_dict = {}
-1573 idl_dict = {}
-1574 for o in list_of_obs:
-1575 new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
-1576 for key in set(o.deltas) | set(o.r_values)})
-1577 idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
-1578
-1579 names = sorted(new_dict.keys())
-1580 o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
-1581 o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
-1582 return o
+ 1568def merge_obs(list_of_obs):
+1569 """Combine all observables in list_of_obs into one new observable
+1570
+1571 Parameters
+1572 ----------
+1573 list_of_obs : list
+1574 list of the Obs object to be combined
+1575
+1576 Notes
+1577 -----
+1578 It is not possible to combine obs which are based on the same replicum
+1579 """
+1580 replist = [item for obs in list_of_obs for item in obs.names]
+1581 if (len(replist) == len(set(replist))) is False:
+1582 raise Exception('list_of_obs contains duplicate replica: %s' % (str(replist)))
+1583 if any([len(o.cov_names) for o in list_of_obs]):
+1584 raise Exception('Not possible to merge data that contains covobs!')
+1585 new_dict = {}
+1586 idl_dict = {}
+1587 for o in list_of_obs:
+1588 new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
+1589 for key in set(o.deltas) | set(o.r_values)})
+1590 idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
+1591
+1592 names = sorted(new_dict.keys())
+1593 o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
+1594 o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
+1595 return o
1585def cov_Obs(means, cov, name, grad=None):
-1586 """Create an Obs based on mean(s) and a covariance matrix
-1587
-1588 Parameters
-1589 ----------
-1590 mean : list of floats or float
-1591 N mean value(s) of the new Obs
-1592 cov : list or array
-1593 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
-1594 name : str
-1595 identifier for the covariance matrix
-1596 grad : list or array
-1597 Gradient of the Covobs wrt. the means belonging to cov.
-1598 """
-1599
-1600 def covobs_to_obs(co):
-1601 """Make an Obs out of a Covobs
-1602
-1603 Parameters
-1604 ----------
-1605 co : Covobs
-1606 Covobs to be embedded into the Obs
-1607 """
-1608 o = Obs([], [], means=[])
-1609 o._value = co.value
-1610 o.names.append(co.name)
-1611 o._covobs[co.name] = co
-1612 o._dvalue = np.sqrt(co.errsq())
-1613 return o
-1614
-1615 ol = []
-1616 if isinstance(means, (float, int)):
-1617 means = [means]
-1618
-1619 for i in range(len(means)):
-1620 ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
-1621 if ol[0].covobs[name].N != len(means):
-1622 raise Exception('You have to provide %d mean values!' % (ol[0].N))
-1623 if len(ol) == 1:
-1624 return ol[0]
-1625 return ol
+ 1598def cov_Obs(means, cov, name, grad=None):
+1599 """Create an Obs based on mean(s) and a covariance matrix
+1600
+1601 Parameters
+1602 ----------
+1603 mean : list of floats or float
+1604 N mean value(s) of the new Obs
+1605 cov : list or array
+1606 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
+1607 name : str
+1608 identifier for the covariance matrix
+1609 grad : list or array
+1610 Gradient of the Covobs wrt. the means belonging to cov.
+1611 """
+1612
+1613 def covobs_to_obs(co):
+1614 """Make an Obs out of a Covobs
+1615
+1616 Parameters
+1617 ----------
+1618 co : Covobs
+1619 Covobs to be embedded into the Obs
+1620 """
+1621 o = Obs([], [], means=[])
+1622 o._value = co.value
+1623 o.names.append(co.name)
+1624 o._covobs[co.name] = co
+1625 o._dvalue = np.sqrt(co.errsq())
+1626 return o
+1627
+1628 ol = []
+1629 if isinstance(means, (float, int)):
+1630 means = [means]
+1631
+1632 for i in range(len(means)):
+1633 ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
+1634 if ol[0].covobs[name].N != len(means):
+1635 raise Exception('You have to provide %d mean values!' % (ol[0].N))
+1636 if len(ol) == 1:
+1637 return ol[0]
+1638 return ol