diff --git a/docs/pyerrors/correlators.html b/docs/pyerrors/correlators.html
index aa01d708..700e56d6 100644
--- a/docs/pyerrors/correlators.html
+++ b/docs/pyerrors/correlators.html
@@ -970,800 +970,795 @@
727 if variant == 'log':
728 newcontent = []
729 for t in range(self.T - 1):
- 730 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
+ 730 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0) or self.content[t][0].value / self.content[t + 1][0].value < 0:
731 newcontent.append(None)
- 732 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
- 733 newcontent.append(None)
- 734 else:
- 735 newcontent.append(self.content[t] / self.content[t + 1])
- 736 if (all([x is None for x in newcontent])):
- 737 raise ValueError('m_eff is undefined at all timeslices')
+ 732 else:
+ 733 newcontent.append(self.content[t] / self.content[t + 1])
+ 734 if (all([x is None for x in newcontent])):
+ 735 raise ValueError('m_eff is undefined at all timeslices')
+ 736
+ 737 return np.log(Corr(newcontent, padding=[0, 1]))
738
- 739 return np.log(Corr(newcontent, padding=[0, 1]))
- 740
- 741 elif variant == 'logsym':
- 742 newcontent = []
- 743 for t in range(1, self.T - 1):
- 744 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
- 745 newcontent.append(None)
- 746 elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
- 747 newcontent.append(None)
- 748 else:
- 749 newcontent.append(self.content[t - 1] / self.content[t + 1])
- 750 if (all([x is None for x in newcontent])):
- 751 raise ValueError('m_eff is undefined at all timeslices')
- 752
- 753 return np.log(Corr(newcontent, padding=[1, 1])) / 2
- 754
- 755 elif variant in ['periodic', 'cosh', 'sinh']:
- 756 if variant in ['periodic', 'cosh']:
- 757 func = anp.cosh
- 758 else:
- 759 func = anp.sinh
- 760
- 761 def root_function(x, d):
- 762 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
- 763
- 764 newcontent = []
- 765 for t in range(self.T - 1):
- 766 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
- 767 newcontent.append(None)
- 768 # Fill the two timeslices in the middle of the lattice with their predecessors
- 769 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
- 770 newcontent.append(newcontent[-1])
- 771 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
- 772 newcontent.append(None)
- 773 else:
- 774 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
- 775 if (all([x is None for x in newcontent])):
- 776 raise ValueError('m_eff is undefined at all timeslices')
- 777
- 778 return Corr(newcontent, padding=[0, 1])
- 779
- 780 elif variant == 'arccosh':
- 781 newcontent = []
- 782 for t in range(1, self.T - 1):
- 783 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0):
- 784 newcontent.append(None)
- 785 else:
- 786 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
- 787 if (all([x is None for x in newcontent])):
- 788 raise ValueError("m_eff is undefined at all timeslices")
- 789 return np.arccosh(Corr(newcontent, padding=[1, 1]))
- 790
- 791 else:
- 792 raise ValueError('Unknown variant.')
- 793
- 794 def fit(self, function, fitrange=None, silent=False, **kwargs):
- 795 r'''Fits function to the data
- 796
- 797 Parameters
- 798 ----------
- 799 function : obj
- 800 function to fit to the data. See fits.least_squares for details.
- 801 fitrange : list
- 802 Two element list containing the timeslices on which the fit is supposed to start and stop.
- 803 Caution: This range is inclusive as opposed to standard python indexing.
- 804 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
- 805 If not specified, self.prange or all timeslices are used.
- 806 silent : bool
- 807 Decides whether output is printed to the standard output.
- 808 '''
- 809 if self.N != 1:
- 810 raise ValueError("Correlator must be projected before fitting")
- 811
- 812 if fitrange is None:
- 813 if self.prange:
- 814 fitrange = self.prange
- 815 else:
- 816 fitrange = [0, self.T - 1]
- 817 else:
- 818 if not isinstance(fitrange, list):
- 819 raise TypeError("fitrange has to be a list with two elements")
- 820 if len(fitrange) != 2:
- 821 raise ValueError("fitrange has to have exactly two elements [fit_start, fit_stop]")
- 822
- 823 xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
- 824 ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
- 825 result = least_squares(xs, ys, function, silent=silent, **kwargs)
- 826 return result
- 827
- 828 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
- 829 """ Extract a plateau value from a Corr object
- 830
- 831 Parameters
- 832 ----------
- 833 plateau_range : list
- 834 list with two entries, indicating the first and the last timeslice
- 835 of the plateau region.
- 836 method : str
- 837 method to extract the plateau.
- 838 'fit' fits a constant to the plateau region
- 839 'avg', 'average' or 'mean' just average over the given timeslices.
- 840 auto_gamma : bool
- 841 apply gamma_method with default parameters to the Corr. Defaults to None
- 842 """
- 843 if not plateau_range:
- 844 if self.prange:
- 845 plateau_range = self.prange
- 846 else:
- 847 raise ValueError("no plateau range provided")
- 848 if self.N != 1:
- 849 raise ValueError("Correlator must be projected before getting a plateau.")
- 850 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
- 851 raise ValueError("plateau is undefined at all timeslices in plateaurange.")
- 852 if auto_gamma:
- 853 self.gamma_method()
- 854 if method == "fit":
- 855 def const_func(a, t):
- 856 return a[0]
- 857 return self.fit(const_func, plateau_range)[0]
- 858 elif method in ["avg", "average", "mean"]:
- 859 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
- 860 return returnvalue
- 861
- 862 else:
- 863 raise ValueError("Unsupported plateau method: " + method)
- 864
- 865 def set_prange(self, prange):
- 866 """Sets the attribute prange of the Corr object."""
- 867 if not len(prange) == 2:
- 868 raise ValueError("prange must be a list or array with two values")
- 869 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
- 870 raise TypeError("Start and end point must be integers")
- 871 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] <= prange[1]):
- 872 raise ValueError("Start and end point must define a range in the interval 0,T")
- 873
- 874 self.prange = prange
- 875 return
- 876
- 877 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
- 878 """Plots the correlator using the tag of the correlator as label if available.
- 879
- 880 Parameters
- 881 ----------
- 882 x_range : list
- 883 list of two values, determining the range of the x-axis e.g. [4, 8].
- 884 comp : Corr or list of Corr
- 885 Correlator or list of correlators which are plotted for comparison.
- 886 The tags of these correlators are used as labels if available.
- 887 logscale : bool
- 888 Sets y-axis to logscale.
- 889 plateau : Obs
- 890 Plateau value to be visualized in the figure.
- 891 fit_res : Fit_result
- 892 Fit_result object to be visualized.
- 893 fit_key : str
- 894 Key for the fit function in Fit_result.fit_function (for combined fits).
- 895 ylabel : str
- 896 Label for the y-axis.
- 897 save : str
- 898 path to file in which the figure should be saved.
- 899 auto_gamma : bool
- 900 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
- 901 hide_sigma : float
- 902 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
- 903 references : list
- 904 List of floating point values that are displayed as horizontal lines for reference.
- 905 title : string
- 906 Optional title of the figure.
- 907 """
- 908 if self.N != 1:
- 909 raise ValueError("Correlator must be projected before plotting")
- 910
- 911 if auto_gamma:
- 912 self.gamma_method()
- 913
- 914 if x_range is None:
- 915 x_range = [0, self.T - 1]
- 916
- 917 fig = plt.figure()
- 918 ax1 = fig.add_subplot(111)
- 919
- 920 x, y, y_err = self.plottable()
- 921 if hide_sigma:
- 922 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+ 739 elif variant == 'logsym':
+ 740 newcontent = []
+ 741 for t in range(1, self.T - 1):
+ 742 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0) or self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
+ 743 newcontent.append(None)
+ 744 else:
+ 745 newcontent.append(self.content[t - 1] / self.content[t + 1])
+ 746 if (all([x is None for x in newcontent])):
+ 747 raise ValueError('m_eff is undefined at all timeslices')
+ 748
+ 749 return np.log(Corr(newcontent, padding=[1, 1])) / 2
+ 750
+ 751 elif variant in ['periodic', 'cosh', 'sinh']:
+ 752 if variant in ['periodic', 'cosh']:
+ 753 func = anp.cosh
+ 754 else:
+ 755 func = anp.sinh
+ 756
+ 757 def root_function(x, d):
+ 758 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
+ 759
+ 760 newcontent = []
+ 761 for t in range(self.T - 1):
+ 762 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
+ 763 newcontent.append(None)
+ 764 # Fill the two timeslices in the middle of the lattice with their predecessors
+ 765 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+ 766 newcontent.append(newcontent[-1])
+ 767 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
+ 768 newcontent.append(None)
+ 769 else:
+ 770 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+ 771 if (all([x is None for x in newcontent])):
+ 772 raise ValueError('m_eff is undefined at all timeslices')
+ 773
+ 774 return Corr(newcontent, padding=[0, 1])
+ 775
+ 776 elif variant == 'arccosh':
+ 777 newcontent = []
+ 778 for t in range(1, self.T - 1):
+ 779 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0):
+ 780 newcontent.append(None)
+ 781 else:
+ 782 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+ 783 if (all([x is None for x in newcontent])):
+ 784 raise ValueError("m_eff is undefined at all timeslices")
+ 785 return np.arccosh(Corr(newcontent, padding=[1, 1]))
+ 786
+ 787 else:
+ 788 raise ValueError('Unknown variant.')
+ 789
+ 790 def fit(self, function, fitrange=None, silent=False, **kwargs):
+ 791 r'''Fits function to the data
+ 792
+ 793 Parameters
+ 794 ----------
+ 795 function : obj
+ 796 function to fit to the data. See fits.least_squares for details.
+ 797 fitrange : list
+ 798 Two element list containing the timeslices on which the fit is supposed to start and stop.
+ 799 Caution: This range is inclusive as opposed to standard python indexing.
+ 800 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+ 801 If not specified, self.prange or all timeslices are used.
+ 802 silent : bool
+ 803 Decides whether output is printed to the standard output.
+ 804 '''
+ 805 if self.N != 1:
+ 806 raise ValueError("Correlator must be projected before fitting")
+ 807
+ 808 if fitrange is None:
+ 809 if self.prange:
+ 810 fitrange = self.prange
+ 811 else:
+ 812 fitrange = [0, self.T - 1]
+ 813 else:
+ 814 if not isinstance(fitrange, list):
+ 815 raise TypeError("fitrange has to be a list with two elements")
+ 816 if len(fitrange) != 2:
+ 817 raise ValueError("fitrange has to have exactly two elements [fit_start, fit_stop]")
+ 818
+ 819 xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
+ 820 ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
+ 821 result = least_squares(xs, ys, function, silent=silent, **kwargs)
+ 822 return result
+ 823
+ 824 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+ 825 """ Extract a plateau value from a Corr object
+ 826
+ 827 Parameters
+ 828 ----------
+ 829 plateau_range : list
+ 830 list with two entries, indicating the first and the last timeslice
+ 831 of the plateau region.
+ 832 method : str
+ 833 method to extract the plateau.
+ 834 'fit' fits a constant to the plateau region
+ 835 'avg', 'average' or 'mean' just average over the given timeslices.
+ 836 auto_gamma : bool
+ 837 apply gamma_method with default parameters to the Corr. Defaults to None
+ 838 """
+ 839 if not plateau_range:
+ 840 if self.prange:
+ 841 plateau_range = self.prange
+ 842 else:
+ 843 raise ValueError("no plateau range provided")
+ 844 if self.N != 1:
+ 845 raise ValueError("Correlator must be projected before getting a plateau.")
+ 846 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+ 847 raise ValueError("plateau is undefined at all timeslices in plateaurange.")
+ 848 if auto_gamma:
+ 849 self.gamma_method()
+ 850 if method == "fit":
+ 851 def const_func(a, t):
+ 852 return a[0]
+ 853 return self.fit(const_func, plateau_range)[0]
+ 854 elif method in ["avg", "average", "mean"]:
+ 855 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+ 856 return returnvalue
+ 857
+ 858 else:
+ 859 raise ValueError("Unsupported plateau method: " + method)
+ 860
+ 861 def set_prange(self, prange):
+ 862 """Sets the attribute prange of the Corr object."""
+ 863 if not len(prange) == 2:
+ 864 raise ValueError("prange must be a list or array with two values")
+ 865 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+ 866 raise TypeError("Start and end point must be integers")
+ 867 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] <= prange[1]):
+ 868 raise ValueError("Start and end point must define a range in the interval 0,T")
+ 869
+ 870 self.prange = prange
+ 871
+ 872 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
+ 873 """Plots the correlator using the tag of the correlator as label if available.
+ 874
+ 875 Parameters
+ 876 ----------
+ 877 x_range : list
+ 878 list of two values, determining the range of the x-axis e.g. [4, 8].
+ 879 comp : Corr or list of Corr
+ 880 Correlator or list of correlators which are plotted for comparison.
+ 881 The tags of these correlators are used as labels if available.
+ 882 logscale : bool
+ 883 Sets y-axis to logscale.
+ 884 plateau : Obs
+ 885 Plateau value to be visualized in the figure.
+ 886 fit_res : Fit_result
+ 887 Fit_result object to be visualized.
+ 888 fit_key : str
+ 889 Key for the fit function in Fit_result.fit_function (for combined fits).
+ 890 ylabel : str
+ 891 Label for the y-axis.
+ 892 save : str
+ 893 path to file in which the figure should be saved.
+ 894 auto_gamma : bool
+ 895 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+ 896 hide_sigma : float
+ 897 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+ 898 references : list
+ 899 List of floating point values that are displayed as horizontal lines for reference.
+ 900 title : string
+ 901 Optional title of the figure.
+ 902 """
+ 903 if self.N != 1:
+ 904 raise ValueError("Correlator must be projected before plotting")
+ 905
+ 906 if auto_gamma:
+ 907 self.gamma_method()
+ 908
+ 909 if x_range is None:
+ 910 x_range = [0, self.T - 1]
+ 911
+ 912 fig = plt.figure()
+ 913 ax1 = fig.add_subplot(111)
+ 914
+ 915 x, y, y_err = self.plottable()
+ 916 if hide_sigma:
+ 917 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+ 918 else:
+ 919 hide_from = None
+ 920 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+ 921 if logscale:
+ 922 ax1.set_yscale('log')
923 else:
- 924 hide_from = None
- 925 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
- 926 if logscale:
- 927 ax1.set_yscale('log')
- 928 else:
- 929 if y_range is None:
- 930 try:
- 931 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
- 932 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
- 933 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
- 934 except Exception:
- 935 pass
- 936 else:
- 937 ax1.set_ylim(y_range)
- 938 if comp:
- 939 if isinstance(comp, (Corr, list)):
- 940 for corr in comp if isinstance(comp, list) else [comp]:
- 941 if auto_gamma:
- 942 corr.gamma_method()
- 943 x, y, y_err = corr.plottable()
- 944 if hide_sigma:
- 945 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
- 946 else:
- 947 hide_from = None
- 948 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
- 949 else:
- 950 raise TypeError("'comp' must be a correlator or a list of correlators.")
- 951
- 952 if plateau:
- 953 if isinstance(plateau, Obs):
- 954 if auto_gamma:
- 955 plateau.gamma_method()
- 956 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
- 957 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
- 958 else:
- 959 raise TypeError("'plateau' must be an Obs")
- 960
- 961 if references:
- 962 if isinstance(references, list):
- 963 for ref in references:
- 964 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
- 965 else:
- 966 raise TypeError("'references' must be a list of floating pint values.")
- 967
- 968 if self.prange:
- 969 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
- 970 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
- 971
- 972 if fit_res:
- 973 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
- 974 if isinstance(fit_res.fit_function, dict):
- 975 if fit_key:
- 976 ax1.plot(x_samples, fit_res.fit_function[fit_key]([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
- 977 else:
- 978 raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
- 979 else:
- 980 ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
+ 924 if y_range is None:
+ 925 try:
+ 926 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+ 927 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+ 928 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+ 929 except Exception:
+ 930 pass
+ 931 else:
+ 932 ax1.set_ylim(y_range)
+ 933 if comp:
+ 934 if isinstance(comp, (Corr, list)):
+ 935 for corr in comp if isinstance(comp, list) else [comp]:
+ 936 if auto_gamma:
+ 937 corr.gamma_method()
+ 938 x, y, y_err = corr.plottable()
+ 939 if hide_sigma:
+ 940 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+ 941 else:
+ 942 hide_from = None
+ 943 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+ 944 else:
+ 945 raise TypeError("'comp' must be a correlator or a list of correlators.")
+ 946
+ 947 if plateau:
+ 948 if isinstance(plateau, Obs):
+ 949 if auto_gamma:
+ 950 plateau.gamma_method()
+ 951 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+ 952 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+ 953 else:
+ 954 raise TypeError("'plateau' must be an Obs")
+ 955
+ 956 if references:
+ 957 if isinstance(references, list):
+ 958 for ref in references:
+ 959 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+ 960 else:
+ 961 raise TypeError("'references' must be a list of floating pint values.")
+ 962
+ 963 if self.prange:
+ 964 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
+ 965 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
+ 966
+ 967 if fit_res:
+ 968 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+ 969 if isinstance(fit_res.fit_function, dict):
+ 970 if fit_key:
+ 971 ax1.plot(x_samples, fit_res.fit_function[fit_key]([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
+ 972 else:
+ 973 raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
+ 974 else:
+ 975 ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
+ 976
+ 977 ax1.set_xlabel(r'$x_0 / a$')
+ 978 if ylabel:
+ 979 ax1.set_ylabel(ylabel)
+ 980 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
981
- 982 ax1.set_xlabel(r'$x_0 / a$')
- 983 if ylabel:
- 984 ax1.set_ylabel(ylabel)
- 985 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
- 986
- 987 _handles, labels = ax1.get_legend_handles_labels()
- 988 if labels:
- 989 ax1.legend()
+ 982 _handles, labels = ax1.get_legend_handles_labels()
+ 983 if labels:
+ 984 ax1.legend()
+ 985
+ 986 if title:
+ 987 plt.title(title)
+ 988
+ 989 plt.draw()
990
- 991 if title:
- 992 plt.title(title)
- 993
- 994 plt.draw()
- 995
- 996 if save:
- 997 if isinstance(save, str):
- 998 fig.savefig(save, bbox_inches='tight')
- 999 else:
-1000 raise TypeError("'save' has to be a string.")
-1001
-1002 def spaghetti_plot(self, logscale=True):
-1003 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-1004
-1005 Parameters
-1006 ----------
-1007 logscale : bool
-1008 Determines whether the scale of the y-axis is logarithmic or standard.
-1009 """
-1010 if self.N != 1:
-1011 raise ValueError("Correlator needs to be projected first.")
-1012
-1013 mc_names = list(set([item for sublist in [list(itertools.chain.from_iterable(map(o[0].e_content.get, o[0].mc_names))) for o in self.content if o is not None] for item in sublist]))
-1014 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content, strict=True) if o is not None]
-1015
-1016 for name in mc_names:
-1017 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+ 991 if save:
+ 992 if isinstance(save, str):
+ 993 fig.savefig(save, bbox_inches='tight')
+ 994 else:
+ 995 raise TypeError("'save' has to be a string.")
+ 996
+ 997 def spaghetti_plot(self, logscale=True):
+ 998 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+ 999
+1000 Parameters
+1001 ----------
+1002 logscale : bool
+1003 Determines whether the scale of the y-axis is logarithmic or standard.
+1004 """
+1005 if self.N != 1:
+1006 raise ValueError("Correlator needs to be projected first.")
+1007
+1008 mc_names = list(set([item for sublist in [list(itertools.chain.from_iterable(map(o[0].e_content.get, o[0].mc_names))) for o in self.content if o is not None] for item in sublist]))
+1009 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content, strict=True) if o is not None]
+1010
+1011 for name in mc_names:
+1012 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+1013
+1014 fig = plt.figure()
+1015 ax = fig.add_subplot(111)
+1016 for dat in data:
+1017 ax.plot(x0_vals, dat, ls='-', marker='')
1018
-1019 fig = plt.figure()
-1020 ax = fig.add_subplot(111)
-1021 for dat in data:
-1022 ax.plot(x0_vals, dat, ls='-', marker='')
-1023
-1024 if logscale is True:
-1025 ax.set_yscale('log')
-1026
-1027 ax.set_xlabel(r'$x_0 / a$')
-1028 plt.title(name)
-1029 plt.draw()
-1030
-1031 def dump(self, filename, datatype="json.gz", **kwargs):
-1032 """Dumps the Corr into a file of chosen type
-1033 Parameters
-1034 ----------
-1035 filename : str
-1036 Name of the file to be saved.
-1037 datatype : str
-1038 Format of the exported file. Supported formats include
-1039 "json.gz" and "pickle"
-1040 path : str
-1041 specifies a custom path for the file (default '.')
-1042 """
-1043 if datatype == "json.gz":
-1044 from .input.json import dump_to_json
-1045 if 'path' in kwargs:
-1046 file_name = kwargs.get('path') + '/' + filename
-1047 else:
-1048 file_name = filename
-1049 dump_to_json(self, file_name)
-1050 elif datatype == "pickle":
-1051 dump_object(self, filename, **kwargs)
-1052 else:
-1053 raise ValueError("Unknown datatype " + str(datatype))
-1054
-1055 def print(self, print_range=None):
-1056 print(self.__repr__(print_range))
-1057
-1058 def __repr__(self, print_range=None):
-1059 if print_range is None:
-1060 print_range = [0, None]
-1061
-1062 content_string = ""
-1063 content_string += "Corr T=" + str(self.T) + " N=" + str(self.N) + "\n" # +" filled with"+ str(type(self.content[0][0])) there should be a good solution here
+1019 if logscale is True:
+1020 ax.set_yscale('log')
+1021
+1022 ax.set_xlabel(r'$x_0 / a$')
+1023 plt.title(name)
+1024 plt.draw()
+1025
+1026 def dump(self, filename, datatype="json.gz", **kwargs):
+1027 """Dumps the Corr into a file of chosen type
+1028 Parameters
+1029 ----------
+1030 filename : str
+1031 Name of the file to be saved.
+1032 datatype : str
+1033 Format of the exported file. Supported formats include
+1034 "json.gz" and "pickle"
+1035 path : str
+1036 specifies a custom path for the file (default '.')
+1037 """
+1038 if datatype == "json.gz":
+1039 from .input.json import dump_to_json
+1040 if 'path' in kwargs:
+1041 file_name = kwargs.get('path') + '/' + filename
+1042 else:
+1043 file_name = filename
+1044 dump_to_json(self, file_name)
+1045 elif datatype == "pickle":
+1046 dump_object(self, filename, **kwargs)
+1047 else:
+1048 raise ValueError("Unknown datatype " + str(datatype))
+1049
+1050 def print(self, print_range=None):
+1051 print(self.__repr__(print_range))
+1052
+1053 def __repr__(self, print_range=None):
+1054 if print_range is None:
+1055 print_range = [0, None]
+1056
+1057 content_string = ""
+1058 content_string += "Corr T=" + str(self.T) + " N=" + str(self.N) + "\n" # +" filled with"+ str(type(self.content[0][0])) there should be a good solution here
+1059
+1060 if self.tag is not None:
+1061 content_string += "Description: " + self.tag + "\n"
+1062 if self.N != 1:
+1063 return content_string
1064
-1065 if self.tag is not None:
-1066 content_string += "Description: " + self.tag + "\n"
-1067 if self.N != 1:
-1068 return content_string
-1069
-1070 if print_range[1]:
-1071 print_range[1] += 1
-1072 content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
-1073 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]):
-1074 if sub_corr is None:
-1075 content_string += str(i + print_range[0]) + '\n'
-1076 else:
-1077 content_string += str(i + print_range[0])
-1078 for element in sub_corr:
-1079 content_string += f"\t{element:+2}"
-1080 content_string += '\n'
-1081 return content_string
-1082
-1083 def __str__(self):
-1084 return self.__repr__()
+1065 if print_range[1]:
+1066 print_range[1] += 1
+1067 content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
+1068 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]):
+1069 if sub_corr is None:
+1070 content_string += str(i + print_range[0]) + '\n'
+1071 else:
+1072 content_string += str(i + print_range[0])
+1073 for element in sub_corr:
+1074 content_string += f"\t{element:+2}"
+1075 content_string += '\n'
+1076 return content_string
+1077
+1078 def __str__(self):
+1079 return self.__repr__()
+1080
+1081 # We define the basic operations, that can be performed with correlators.
+1082 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
+1083 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
+1084 # One could try and tell Obs to check if the y in __mul__ is a Corr and
1085
-1086 # We define the basic operations, that can be performed with correlators.
-1087 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
-1088 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
-1089 # One could try and tell Obs to check if the y in __mul__ is a Corr and
-1090
-1091 __array_priority__ = 10000
-1092
-1093 def __eq__(self, y):
-1094 if isinstance(y, Corr):
-1095 comp = np.asarray(y.content, dtype=object)
-1096 else:
-1097 comp = np.asarray(y)
-1098 return np.asarray(self.content, dtype=object) == comp
-1099
-1100 __hash__ = None
-1101
-1102 def __add__(self, y):
-1103 if isinstance(y, Corr):
-1104 if ((self.N != y.N) or (self.T != y.T)):
-1105 raise ValueError("Addition of Corrs with different shape")
-1106 newcontent = []
-1107 for t in range(self.T):
-1108 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
-1109 newcontent.append(None)
-1110 else:
-1111 newcontent.append(self.content[t] + y.content[t])
-1112 return Corr(newcontent)
-1113
-1114 elif isinstance(y, (Obs, int, float, CObs, complex)):
-1115 newcontent = []
-1116 for t in range(self.T):
-1117 if _check_for_none(self, self.content[t]):
-1118 newcontent.append(None)
-1119 else:
-1120 newcontent.append(self.content[t] + y)
-1121 return Corr(newcontent, prange=self.prange)
-1122 elif isinstance(y, np.ndarray):
-1123 if y.shape == (self.T,):
-1124 return Corr(list((np.array(self.content).T + y).T))
-1125 else:
-1126 raise ValueError("operands could not be broadcast together")
-1127 else:
-1128 raise TypeError("Corr + wrong type")
-1129
-1130 def __mul__(self, y):
-1131 if isinstance(y, Corr):
-1132 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
-1133 raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
-1134 newcontent = []
-1135 for t in range(self.T):
-1136 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
-1137 newcontent.append(None)
-1138 else:
-1139 newcontent.append(self.content[t] * y.content[t])
-1140 return Corr(newcontent)
-1141
-1142 elif isinstance(y, (Obs, int, float, CObs, complex)):
-1143 newcontent = []
-1144 for t in range(self.T):
-1145 if _check_for_none(self, self.content[t]):
-1146 newcontent.append(None)
-1147 else:
-1148 newcontent.append(self.content[t] * y)
-1149 return Corr(newcontent, prange=self.prange)
-1150 elif isinstance(y, np.ndarray):
-1151 if y.shape == (self.T,):
-1152 return Corr(list((np.array(self.content).T * y).T))
-1153 else:
-1154 raise ValueError("operands could not be broadcast together")
-1155 else:
-1156 raise TypeError("Corr * wrong type")
-1157
-1158 def __matmul__(self, y):
-1159 if isinstance(y, np.ndarray):
-1160 if y.ndim != 2 or y.shape[0] != y.shape[1]:
-1161 raise ValueError("Can only multiply correlators by square matrices.")
-1162 if not self.N == y.shape[0]:
-1163 raise ValueError("matmul: mismatch of matrix dimensions")
-1164 newcontent = []
-1165 for t in range(self.T):
-1166 if _check_for_none(self, self.content[t]):
-1167 newcontent.append(None)
-1168 else:
-1169 newcontent.append(self.content[t] @ y)
-1170 return Corr(newcontent)
-1171 elif isinstance(y, Corr):
-1172 if not self.N == y.N:
-1173 raise ValueError("matmul: mismatch of matrix dimensions")
-1174 newcontent = []
-1175 for t in range(self.T):
-1176 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
-1177 newcontent.append(None)
-1178 else:
-1179 newcontent.append(self.content[t] @ y.content[t])
-1180 return Corr(newcontent)
-1181
-1182 else:
-1183 return NotImplemented
-1184
-1185 def __rmatmul__(self, y):
-1186 if isinstance(y, np.ndarray):
-1187 if y.ndim != 2 or y.shape[0] != y.shape[1]:
-1188 raise ValueError("Can only multiply correlators by square matrices.")
-1189 if not self.N == y.shape[0]:
-1190 raise ValueError("matmul: mismatch of matrix dimensions")
-1191 newcontent = []
-1192 for t in range(self.T):
-1193 if _check_for_none(self, self.content[t]):
-1194 newcontent.append(None)
-1195 else:
-1196 newcontent.append(y @ self.content[t])
-1197 return Corr(newcontent)
-1198 else:
-1199 return NotImplemented
-1200
-1201 def __truediv__(self, y):
-1202 if isinstance(y, Corr):
-1203 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
-1204 raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
-1205 newcontent = []
+1086 __array_priority__ = 10000
+1087
+1088 def __eq__(self, y):
+1089 if isinstance(y, Corr):
+1090 comp = np.asarray(y.content, dtype=object)
+1091 else:
+1092 comp = np.asarray(y)
+1093 return np.asarray(self.content, dtype=object) == comp
+1094
+1095 __hash__ = None
+1096
+1097 def __add__(self, y):
+1098 if isinstance(y, Corr):
+1099 if ((self.N != y.N) or (self.T != y.T)):
+1100 raise ValueError("Addition of Corrs with different shape")
+1101 newcontent = []
+1102 for t in range(self.T):
+1103 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
+1104 newcontent.append(None)
+1105 else:
+1106 newcontent.append(self.content[t] + y.content[t])
+1107 return Corr(newcontent)
+1108
+1109 elif isinstance(y, (Obs, int, float, CObs, complex)):
+1110 newcontent = []
+1111 for t in range(self.T):
+1112 if _check_for_none(self, self.content[t]):
+1113 newcontent.append(None)
+1114 else:
+1115 newcontent.append(self.content[t] + y)
+1116 return Corr(newcontent, prange=self.prange)
+1117 elif isinstance(y, np.ndarray):
+1118 if y.shape == (self.T,):
+1119 return Corr(list((np.array(self.content).T + y).T))
+1120 else:
+1121 raise ValueError("operands could not be broadcast together")
+1122 else:
+1123 raise TypeError("Corr + wrong type")
+1124
+1125 def __mul__(self, y):
+1126 if isinstance(y, Corr):
+1127 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
+1128 raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
+1129 newcontent = []
+1130 for t in range(self.T):
+1131 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
+1132 newcontent.append(None)
+1133 else:
+1134 newcontent.append(self.content[t] * y.content[t])
+1135 return Corr(newcontent)
+1136
+1137 elif isinstance(y, (Obs, int, float, CObs, complex)):
+1138 newcontent = []
+1139 for t in range(self.T):
+1140 if _check_for_none(self, self.content[t]):
+1141 newcontent.append(None)
+1142 else:
+1143 newcontent.append(self.content[t] * y)
+1144 return Corr(newcontent, prange=self.prange)
+1145 elif isinstance(y, np.ndarray):
+1146 if y.shape == (self.T,):
+1147 return Corr(list((np.array(self.content).T * y).T))
+1148 else:
+1149 raise ValueError("operands could not be broadcast together")
+1150 else:
+1151 raise TypeError("Corr * wrong type")
+1152
+1153 def __matmul__(self, y):
+1154 if isinstance(y, np.ndarray):
+1155 if y.ndim != 2 or y.shape[0] != y.shape[1]:
+1156 raise ValueError("Can only multiply correlators by square matrices.")
+1157 if not self.N == y.shape[0]:
+1158 raise ValueError("matmul: mismatch of matrix dimensions")
+1159 newcontent = []
+1160 for t in range(self.T):
+1161 if _check_for_none(self, self.content[t]):
+1162 newcontent.append(None)
+1163 else:
+1164 newcontent.append(self.content[t] @ y)
+1165 return Corr(newcontent)
+1166 elif isinstance(y, Corr):
+1167 if not self.N == y.N:
+1168 raise ValueError("matmul: mismatch of matrix dimensions")
+1169 newcontent = []
+1170 for t in range(self.T):
+1171 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
+1172 newcontent.append(None)
+1173 else:
+1174 newcontent.append(self.content[t] @ y.content[t])
+1175 return Corr(newcontent)
+1176
+1177 else:
+1178 return NotImplemented
+1179
+1180 def __rmatmul__(self, y):
+1181 if isinstance(y, np.ndarray):
+1182 if y.ndim != 2 or y.shape[0] != y.shape[1]:
+1183 raise ValueError("Can only multiply correlators by square matrices.")
+1184 if not self.N == y.shape[0]:
+1185 raise ValueError("matmul: mismatch of matrix dimensions")
+1186 newcontent = []
+1187 for t in range(self.T):
+1188 if _check_for_none(self, self.content[t]):
+1189 newcontent.append(None)
+1190 else:
+1191 newcontent.append(y @ self.content[t])
+1192 return Corr(newcontent)
+1193 else:
+1194 return NotImplemented
+1195
+1196 def __truediv__(self, y):
+1197 if isinstance(y, Corr):
+1198 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
+1199 raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
+1200 newcontent = []
+1201 for t in range(self.T):
+1202 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
+1203 newcontent.append(None)
+1204 else:
+1205 newcontent.append(self.content[t] / y.content[t])
1206 for t in range(self.T):
-1207 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
-1208 newcontent.append(None)
-1209 else:
-1210 newcontent.append(self.content[t] / y.content[t])
-1211 for t in range(self.T):
-1212 if _check_for_none(self, newcontent[t]):
-1213 continue
-1214 if np.isnan(np.sum(newcontent[t]).value):
-1215 newcontent[t] = None
-1216
-1217 if all([item is None for item in newcontent]):
-1218 raise ValueError("Division returns completely undefined correlator")
-1219 return Corr(newcontent)
-1220
-1221 elif isinstance(y, (Obs, CObs)):
-1222 if isinstance(y, Obs):
-1223 if y.value == 0:
-1224 raise ValueError('Division by zero will return undefined correlator')
-1225 if isinstance(y, CObs):
-1226 if y.is_zero():
-1227 raise ValueError('Division by zero will return undefined correlator')
-1228
-1229 newcontent = []
-1230 for t in range(self.T):
-1231 if _check_for_none(self, self.content[t]):
-1232 newcontent.append(None)
-1233 else:
-1234 newcontent.append(self.content[t] / y)
-1235 return Corr(newcontent, prange=self.prange)
-1236
-1237 elif isinstance(y, (int, float)):
-1238 if y == 0:
-1239 raise ValueError('Division by zero will return undefined correlator')
-1240 newcontent = []
-1241 for t in range(self.T):
-1242 if _check_for_none(self, self.content[t]):
-1243 newcontent.append(None)
-1244 else:
-1245 newcontent.append(self.content[t] / y)
-1246 return Corr(newcontent, prange=self.prange)
-1247 elif isinstance(y, np.ndarray):
-1248 if y.shape == (self.T,):
-1249 return Corr(list((np.array(self.content).T / y).T))
-1250 else:
-1251 raise ValueError("operands could not be broadcast together")
-1252 else:
-1253 raise TypeError('Corr / wrong type')
-1254
-1255 def __neg__(self):
-1256 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
-1257 return Corr(newcontent, prange=self.prange)
-1258
-1259 def __sub__(self, y):
-1260 return self + (-y)
-1261
-1262 def __pow__(self, y):
-1263 if isinstance(y, (Obs, int, float, CObs)):
-1264 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
-1265 return Corr(newcontent, prange=self.prange)
-1266 else:
-1267 raise TypeError('Type of exponent not supported')
-1268
-1269 def __abs__(self):
-1270 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
-1271 return Corr(newcontent, prange=self.prange)
-1272
-1273 # The numpy functions:
-1274 def sqrt(self):
-1275 return self ** 0.5
-1276
-1277 def log(self):
-1278 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
-1279 return Corr(newcontent, prange=self.prange)
-1280
-1281 def exp(self):
-1282 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
-1283 return Corr(newcontent, prange=self.prange)
-1284
-1285 def _apply_func_to_corr(self, func):
-1286 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
-1287 for t in range(self.T):
-1288 if _check_for_none(self, newcontent[t]):
-1289 continue
-1290 tmp_sum = np.sum(newcontent[t])
-1291 if hasattr(tmp_sum, "value"):
-1292 if np.isnan(tmp_sum.value):
-1293 newcontent[t] = None
-1294 if all([item is None for item in newcontent]):
-1295 raise ValueError('Operation returns undefined correlator')
-1296 return Corr(newcontent)
-1297
-1298 def sin(self):
-1299 return self._apply_func_to_corr(np.sin)
-1300
-1301 def cos(self):
-1302 return self._apply_func_to_corr(np.cos)
-1303
-1304 def tan(self):
-1305 return self._apply_func_to_corr(np.tan)
-1306
-1307 def sinh(self):
-1308 return self._apply_func_to_corr(np.sinh)
-1309
-1310 def cosh(self):
-1311 return self._apply_func_to_corr(np.cosh)
-1312
-1313 def tanh(self):
-1314 return self._apply_func_to_corr(np.tanh)
-1315
-1316 def arcsin(self):
-1317 return self._apply_func_to_corr(np.arcsin)
-1318
-1319 def arccos(self):
-1320 return self._apply_func_to_corr(np.arccos)
-1321
-1322 def arctan(self):
-1323 return self._apply_func_to_corr(np.arctan)
-1324
-1325 def arcsinh(self):
-1326 return self._apply_func_to_corr(np.arcsinh)
-1327
-1328 def arccosh(self):
-1329 return self._apply_func_to_corr(np.arccosh)
-1330
-1331 def arctanh(self):
-1332 return self._apply_func_to_corr(np.arctanh)
-1333
-1334 # Right hand side operations (require tweak in main module to work)
-1335 def __radd__(self, y):
-1336 return self + y
-1337
-1338 def __rsub__(self, y):
-1339 return -self + y
-1340
-1341 def __rmul__(self, y):
-1342 return self * y
-1343
-1344 def __rtruediv__(self, y):
-1345 return (self / y) ** (-1)
-1346
-1347 @property
-1348 def real(self):
-1349 def return_real(obs_OR_cobs):
-1350 if isinstance(obs_OR_cobs.flatten()[0], CObs):
-1351 return np.vectorize(lambda x: x.real)(obs_OR_cobs)
-1352 else:
-1353 return obs_OR_cobs
-1354
-1355 return self._apply_func_to_corr(return_real)
-1356
-1357 @property
-1358 def imag(self):
-1359 def return_imag(obs_OR_cobs):
-1360 if isinstance(obs_OR_cobs.flatten()[0], CObs):
-1361 return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
-1362 else:
-1363 return obs_OR_cobs * 0 # So it stays the right type
+1207 if _check_for_none(self, newcontent[t]):
+1208 continue
+1209 if np.isnan(np.sum(newcontent[t]).value):
+1210 newcontent[t] = None
+1211
+1212 if all([item is None for item in newcontent]):
+1213 raise ValueError("Division returns completely undefined correlator")
+1214 return Corr(newcontent)
+1215
+1216 elif isinstance(y, (Obs, CObs)):
+1217 if isinstance(y, Obs):
+1218 if y.value == 0:
+1219 raise ValueError('Division by zero will return undefined correlator')
+1220 if isinstance(y, CObs):
+1221 if y.is_zero():
+1222 raise ValueError('Division by zero will return undefined correlator')
+1223
+1224 newcontent = []
+1225 for t in range(self.T):
+1226 if _check_for_none(self, self.content[t]):
+1227 newcontent.append(None)
+1228 else:
+1229 newcontent.append(self.content[t] / y)
+1230 return Corr(newcontent, prange=self.prange)
+1231
+1232 elif isinstance(y, (int, float)):
+1233 if y == 0:
+1234 raise ValueError('Division by zero will return undefined correlator')
+1235 newcontent = []
+1236 for t in range(self.T):
+1237 if _check_for_none(self, self.content[t]):
+1238 newcontent.append(None)
+1239 else:
+1240 newcontent.append(self.content[t] / y)
+1241 return Corr(newcontent, prange=self.prange)
+1242 elif isinstance(y, np.ndarray):
+1243 if y.shape == (self.T,):
+1244 return Corr(list((np.array(self.content).T / y).T))
+1245 else:
+1246 raise ValueError("operands could not be broadcast together")
+1247 else:
+1248 raise TypeError('Corr / wrong type')
+1249
+1250 def __neg__(self):
+1251 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
+1252 return Corr(newcontent, prange=self.prange)
+1253
+1254 def __sub__(self, y):
+1255 return self + (-y)
+1256
+1257 def __pow__(self, y):
+1258 if isinstance(y, (Obs, int, float, CObs)):
+1259 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
+1260 return Corr(newcontent, prange=self.prange)
+1261 else:
+1262 raise TypeError('Type of exponent not supported')
+1263
+1264 def __abs__(self):
+1265 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
+1266 return Corr(newcontent, prange=self.prange)
+1267
+1268 # The numpy functions:
+1269 def sqrt(self):
+1270 return self ** 0.5
+1271
+1272 def log(self):
+1273 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
+1274 return Corr(newcontent, prange=self.prange)
+1275
+1276 def exp(self):
+1277 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
+1278 return Corr(newcontent, prange=self.prange)
+1279
+1280 def _apply_func_to_corr(self, func):
+1281 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
+1282 for t in range(self.T):
+1283 if _check_for_none(self, newcontent[t]):
+1284 continue
+1285 tmp_sum = np.sum(newcontent[t])
+1286 if hasattr(tmp_sum, "value"):
+1287 if np.isnan(tmp_sum.value):
+1288 newcontent[t] = None
+1289 if all([item is None for item in newcontent]):
+1290 raise ValueError('Operation returns undefined correlator')
+1291 return Corr(newcontent)
+1292
+1293 def sin(self):
+1294 return self._apply_func_to_corr(np.sin)
+1295
+1296 def cos(self):
+1297 return self._apply_func_to_corr(np.cos)
+1298
+1299 def tan(self):
+1300 return self._apply_func_to_corr(np.tan)
+1301
+1302 def sinh(self):
+1303 return self._apply_func_to_corr(np.sinh)
+1304
+1305 def cosh(self):
+1306 return self._apply_func_to_corr(np.cosh)
+1307
+1308 def tanh(self):
+1309 return self._apply_func_to_corr(np.tanh)
+1310
+1311 def arcsin(self):
+1312 return self._apply_func_to_corr(np.arcsin)
+1313
+1314 def arccos(self):
+1315 return self._apply_func_to_corr(np.arccos)
+1316
+1317 def arctan(self):
+1318 return self._apply_func_to_corr(np.arctan)
+1319
+1320 def arcsinh(self):
+1321 return self._apply_func_to_corr(np.arcsinh)
+1322
+1323 def arccosh(self):
+1324 return self._apply_func_to_corr(np.arccosh)
+1325
+1326 def arctanh(self):
+1327 return self._apply_func_to_corr(np.arctanh)
+1328
+1329 # Right hand side operations (require tweak in main module to work)
+1330 def __radd__(self, y):
+1331 return self + y
+1332
+1333 def __rsub__(self, y):
+1334 return -self + y
+1335
+1336 def __rmul__(self, y):
+1337 return self * y
+1338
+1339 def __rtruediv__(self, y):
+1340 return (self / y) ** (-1)
+1341
+1342 @property
+1343 def real(self):
+1344 def return_real(obs_OR_cobs):
+1345 if isinstance(obs_OR_cobs.flatten()[0], CObs):
+1346 return np.vectorize(lambda x: x.real)(obs_OR_cobs)
+1347 else:
+1348 return obs_OR_cobs
+1349
+1350 return self._apply_func_to_corr(return_real)
+1351
+1352 @property
+1353 def imag(self):
+1354 def return_imag(obs_OR_cobs):
+1355 if isinstance(obs_OR_cobs.flatten()[0], CObs):
+1356 return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
+1357 else:
+1358 return obs_OR_cobs * 0 # So it stays the right type
+1359
+1360 return self._apply_func_to_corr(return_imag)
+1361
+1362 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1363 r''' Project large correlation matrix to lowest states
1364
-1365 return self._apply_func_to_corr(return_imag)
-1366
-1367 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1368 r''' Project large correlation matrix to lowest states
-1369
-1370 This method can be used to reduce the size of an (N x N) correlation matrix
-1371 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1372 is still small.
-1373
-1374 Parameters
-1375 ----------
-1376 Ntrunc: int
-1377 Rank of the target matrix.
-1378 tproj: int
-1379 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1380 The default value is 3.
-1381 t0proj: int
-1382 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1383 discouraged for O(a) improved theories, since the correctness of the procedure
-1384 cannot be granted in this case. The default value is 2.
-1385 basematrix : Corr
-1386 Correlation matrix that is used to determine the eigenvectors of the
-1387 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1388 is is not specified.
-1389
-1390 Notes
-1391 -----
-1392 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1393 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
-1394 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1395 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1396 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1397 correlation matrix and to remove some noise that is added by irrelevant operators.
-1398 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1399 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1400 '''
-1401
-1402 if self.N == 1:
-1403 raise ValueError('Method cannot be applied to one-dimensional correlators.')
-1404 if basematrix is None:
-1405 basematrix = self
-1406 if Ntrunc >= basematrix.N:
-1407 raise ValueError(f'Cannot truncate using Ntrunc >= {basematrix.N}')
-1408 if basematrix.N != self.N:
-1409 raise ValueError('basematrix and targetmatrix have to be of the same size.')
-1410
-1411 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
-1412
-1413 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1414 rmat = []
-1415 for t in range(basematrix.T):
-1416 if self.content[t] is None:
-1417 rmat.append(None)
-1418 else:
-1419 for i in range(Ntrunc):
-1420 for j in range(Ntrunc):
-1421 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1422 rmat.append(np.copy(tmpmat))
-1423
-1424 return Corr(rmat)
-1425
-1426
-1427def _sort_vectors(vec_set_in, ts):
-1428 """Helper function used to find a set of Eigenvectors consistent over all timeslices"""
-1429
-1430 if isinstance(vec_set_in[ts][0][0], Obs):
-1431 vec_set = [anp.vectorize(float)(vi) if vi is not None else vi for vi in vec_set_in]
-1432 else:
-1433 vec_set = vec_set_in
-1434 reference_sorting = np.array(vec_set[ts])
-1435 N = reference_sorting.shape[0]
-1436 sorted_vec_set = []
-1437 for t in range(len(vec_set)):
-1438 if vec_set[t] is None:
-1439 sorted_vec_set.append(None)
-1440 elif not t == ts:
-1441 perms = [list(o) for o in permutations([i for i in range(N)], N)]
-1442 best_score = 0
-1443 for perm in perms:
-1444 current_score = 1
-1445 for k in range(N):
-1446 new_sorting = reference_sorting.copy()
-1447 new_sorting[perm[k], :] = vec_set[t][k]
-1448 current_score *= abs(np.linalg.det(new_sorting))
-1449 if current_score > best_score:
-1450 best_score = current_score
-1451 best_perm = perm
-1452 sorted_vec_set.append([vec_set_in[t][k] for k in best_perm])
-1453 else:
-1454 sorted_vec_set.append(vec_set_in[t])
-1455
-1456 return sorted_vec_set
+1365 This method can be used to reduce the size of an (N x N) correlation matrix
+1366 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1367 is still small.
+1368
+1369 Parameters
+1370 ----------
+1371 Ntrunc: int
+1372 Rank of the target matrix.
+1373 tproj: int
+1374 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1375 The default value is 3.
+1376 t0proj: int
+1377 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1378 discouraged for O(a) improved theories, since the correctness of the procedure
+1379 cannot be granted in this case. The default value is 2.
+1380 basematrix : Corr
+1381 Correlation matrix that is used to determine the eigenvectors of the
+1382 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1383 is is not specified.
+1384
+1385 Notes
+1386 -----
+1387 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1388 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
+1389 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1390 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1391 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1392 correlation matrix and to remove some noise that is added by irrelevant operators.
+1393 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1394 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1395 '''
+1396
+1397 if self.N == 1:
+1398 raise ValueError('Method cannot be applied to one-dimensional correlators.')
+1399 if basematrix is None:
+1400 basematrix = self
+1401 if Ntrunc >= basematrix.N:
+1402 raise ValueError(f'Cannot truncate using Ntrunc >= {basematrix.N}')
+1403 if basematrix.N != self.N:
+1404 raise ValueError('basematrix and targetmatrix have to be of the same size.')
+1405
+1406 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
+1407
+1408 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1409 rmat = []
+1410 for t in range(basematrix.T):
+1411 if self.content[t] is None:
+1412 rmat.append(None)
+1413 else:
+1414 for i in range(Ntrunc):
+1415 for j in range(Ntrunc):
+1416 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1417 rmat.append(np.copy(tmpmat))
+1418
+1419 return Corr(rmat)
+1420
+1421
+1422def _sort_vectors(vec_set_in, ts):
+1423 """Helper function used to find a set of Eigenvectors consistent over all timeslices"""
+1424
+1425 if isinstance(vec_set_in[ts][0][0], Obs):
+1426 vec_set = [anp.vectorize(float)(vi) if vi is not None else vi for vi in vec_set_in]
+1427 else:
+1428 vec_set = vec_set_in
+1429 reference_sorting = np.array(vec_set[ts])
+1430 N = reference_sorting.shape[0]
+1431 sorted_vec_set = []
+1432 for t in range(len(vec_set)):
+1433 if vec_set[t] is None:
+1434 sorted_vec_set.append(None)
+1435 elif not t == ts:
+1436 perms = [list(o) for o in permutations([i for i in range(N)], N)]
+1437 best_score = 0
+1438 for perm in perms:
+1439 current_score = 1
+1440 for k in range(N):
+1441 new_sorting = reference_sorting.copy()
+1442 new_sorting[perm[k], :] = vec_set[t][k]
+1443 current_score *= abs(np.linalg.det(new_sorting))
+1444 if current_score > best_score:
+1445 best_score = current_score
+1446 best_perm = perm
+1447 sorted_vec_set.append([vec_set_in[t][k] for k in best_perm])
+1448 else:
+1449 sorted_vec_set.append(vec_set_in[t])
+1450
+1451 return sorted_vec_set
+1452
+1453
+1454def _check_for_none(corr, entry):
+1455 """Checks if entry for correlator corr is None"""
+1456 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2
1457
1458
-1459def _check_for_none(corr, entry):
-1460 """Checks if entry for correlator corr is None"""
-1461 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2
-1462
+1459def _GEVP_solver(Gt, G0, method='eigh', chol_inv=None):
+1460 r"""Helper function for solving the GEVP and sorting the eigenvectors.
+1461
+1462 Solves $G(t)v_i=\lambda_i G(t_0)v_i$ and returns the eigenvectors v_i
1463
-1464def _GEVP_solver(Gt, G0, method='eigh', chol_inv=None):
-1465 r"""Helper function for solving the GEVP and sorting the eigenvectors.
-1466
-1467 Solves $G(t)v_i=\lambda_i G(t_0)v_i$ and returns the eigenvectors v_i
-1468
-1469 The helper function assumes that both provided matrices are symmetric and
-1470 only processes the lower triangular part of both matrices. In case the matrices
-1471 are not symmetric the upper triangular parts are effectively discarded.
-1472
-1473 Parameters
-1474 ----------
-1475 Gt : array
-1476 The correlator at time t for the left hand side of the GEVP
-1477 G0 : array
-1478 The correlator at time t0 for the right hand side of the GEVP
-1479 Method used to solve the GEVP.
-1480 - "eigh": Use scipy.linalg.eigh to solve the GEVP.
-1481 - "cholesky": Use manually implemented solution via the Cholesky decomposition.
-1482 chol_inv : array, optional
-1483 Inverse of the Cholesky decomposition of G0. May be provided to
-1484 speed up the computation in the case of method=='cholesky'
-1485
-1486 """
-1487 if isinstance(G0[0][0], Obs):
-1488 vector_obs = True
-1489 else:
-1490 vector_obs = False
-1491
-1492 if method == 'cholesky':
-1493 if vector_obs:
-1494 cholesky = linalg.cholesky
-1495 inv = linalg.inv
-1496 eigv = linalg.eigv
-1497 matmul = linalg.matmul
-1498 else:
-1499 cholesky = np.linalg.cholesky
-1500 inv = np.linalg.inv
-1501
-1502 def eigv(x, **kwargs):
-1503 return np.linalg.eigh(x)[1]
-1504
-1505 def matmul(*operands):
-1506 return np.linalg.multi_dot(operands)
-1507 N = Gt.shape[0]
-1508 output = [[] for j in range(N)]
-1509 if chol_inv is None:
-1510 chol = cholesky(G0) # This will automatically report if the matrix is not pos-def
-1511 chol_inv = inv(chol)
-1512
-1513 try:
-1514 new_matrix = matmul(chol_inv, Gt, chol_inv.T)
-1515 ev = eigv(new_matrix)
-1516 ev = matmul(chol_inv.T, ev)
-1517 output = np.flip(ev, axis=1).T
-1518 except (np.linalg.LinAlgError, TypeError, ValueError): # The above code can fail because of linalg-errors or because the entry of the corr is None
-1519 for s in range(N):
-1520 output[s] = None
-1521 return output
-1522 elif method == 'eigh':
-1523 return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1]
+1464 The helper function assumes that both provided matrices are symmetric and
+1465 only processes the lower triangular part of both matrices. In case the matrices
+1466 are not symmetric the upper triangular parts are effectively discarded.
+1467
+1468 Parameters
+1469 ----------
+1470 Gt : array
+1471 The correlator at time t for the left hand side of the GEVP
+1472 G0 : array
+1473 The correlator at time t0 for the right hand side of the GEVP
+1474 Method used to solve the GEVP.
+1475 - "eigh": Use scipy.linalg.eigh to solve the GEVP.
+1476 - "cholesky": Use manually implemented solution via the Cholesky decomposition.
+1477 chol_inv : array, optional
+1478 Inverse of the Cholesky decomposition of G0. May be provided to
+1479 speed up the computation in the case of method=='cholesky'
+1480
+1481 """
+1482 if isinstance(G0[0][0], Obs):
+1483 vector_obs = True
+1484 else:
+1485 vector_obs = False
+1486
+1487 if method == 'cholesky':
+1488 if vector_obs:
+1489 cholesky = linalg.cholesky
+1490 inv = linalg.inv
+1491 eigv = linalg.eigv
+1492 matmul = linalg.matmul
+1493 else:
+1494 cholesky = np.linalg.cholesky
+1495 inv = np.linalg.inv
+1496
+1497 def eigv(x, **kwargs):
+1498 return np.linalg.eigh(x)[1]
+1499
+1500 def matmul(*operands):
+1501 return np.linalg.multi_dot(operands)
+1502 N = Gt.shape[0]
+1503 output = [[] for j in range(N)]
+1504 if chol_inv is None:
+1505 chol = cholesky(G0) # This will automatically report if the matrix is not pos-def
+1506 chol_inv = inv(chol)
+1507
+1508 try:
+1509 new_matrix = matmul(chol_inv, Gt, chol_inv.T)
+1510 ev = eigv(new_matrix)
+1511 ev = matmul(chol_inv.T, ev)
+1512 output = np.flip(ev, axis=1).T
+1513 except (np.linalg.LinAlgError, TypeError, ValueError): # The above code can fail because of linalg-errors or because the entry of the corr is None
+1514 for s in range(N):
+1515 output[s] = None
+1516 return output
+1517 elif method == 'eigh':
+1518 return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1]
@@ -2492,701 +2487,696 @@
728 if variant == 'log':
729 newcontent = []
730 for t in range(self.T - 1):
- 731 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
+ 731 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0) or self.content[t][0].value / self.content[t + 1][0].value < 0:
732 newcontent.append(None)
- 733 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
- 734 newcontent.append(None)
- 735 else:
- 736 newcontent.append(self.content[t] / self.content[t + 1])
- 737 if (all([x is None for x in newcontent])):
- 738 raise ValueError('m_eff is undefined at all timeslices')
+ 733 else:
+ 734 newcontent.append(self.content[t] / self.content[t + 1])
+ 735 if (all([x is None for x in newcontent])):
+ 736 raise ValueError('m_eff is undefined at all timeslices')
+ 737
+ 738 return np.log(Corr(newcontent, padding=[0, 1]))
739
- 740 return np.log(Corr(newcontent, padding=[0, 1]))
- 741
- 742 elif variant == 'logsym':
- 743 newcontent = []
- 744 for t in range(1, self.T - 1):
- 745 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
- 746 newcontent.append(None)
- 747 elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
- 748 newcontent.append(None)
- 749 else:
- 750 newcontent.append(self.content[t - 1] / self.content[t + 1])
- 751 if (all([x is None for x in newcontent])):
- 752 raise ValueError('m_eff is undefined at all timeslices')
- 753
- 754 return np.log(Corr(newcontent, padding=[1, 1])) / 2
- 755
- 756 elif variant in ['periodic', 'cosh', 'sinh']:
- 757 if variant in ['periodic', 'cosh']:
- 758 func = anp.cosh
- 759 else:
- 760 func = anp.sinh
- 761
- 762 def root_function(x, d):
- 763 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
- 764
- 765 newcontent = []
- 766 for t in range(self.T - 1):
- 767 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
- 768 newcontent.append(None)
- 769 # Fill the two timeslices in the middle of the lattice with their predecessors
- 770 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
- 771 newcontent.append(newcontent[-1])
- 772 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
- 773 newcontent.append(None)
- 774 else:
- 775 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
- 776 if (all([x is None for x in newcontent])):
- 777 raise ValueError('m_eff is undefined at all timeslices')
- 778
- 779 return Corr(newcontent, padding=[0, 1])
- 780
- 781 elif variant == 'arccosh':
- 782 newcontent = []
- 783 for t in range(1, self.T - 1):
- 784 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0):
- 785 newcontent.append(None)
- 786 else:
- 787 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
- 788 if (all([x is None for x in newcontent])):
- 789 raise ValueError("m_eff is undefined at all timeslices")
- 790 return np.arccosh(Corr(newcontent, padding=[1, 1]))
- 791
- 792 else:
- 793 raise ValueError('Unknown variant.')
- 794
- 795 def fit(self, function, fitrange=None, silent=False, **kwargs):
- 796 r'''Fits function to the data
- 797
- 798 Parameters
- 799 ----------
- 800 function : obj
- 801 function to fit to the data. See fits.least_squares for details.
- 802 fitrange : list
- 803 Two element list containing the timeslices on which the fit is supposed to start and stop.
- 804 Caution: This range is inclusive as opposed to standard python indexing.
- 805 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
- 806 If not specified, self.prange or all timeslices are used.
- 807 silent : bool
- 808 Decides whether output is printed to the standard output.
- 809 '''
- 810 if self.N != 1:
- 811 raise ValueError("Correlator must be projected before fitting")
- 812
- 813 if fitrange is None:
- 814 if self.prange:
- 815 fitrange = self.prange
- 816 else:
- 817 fitrange = [0, self.T - 1]
- 818 else:
- 819 if not isinstance(fitrange, list):
- 820 raise TypeError("fitrange has to be a list with two elements")
- 821 if len(fitrange) != 2:
- 822 raise ValueError("fitrange has to have exactly two elements [fit_start, fit_stop]")
- 823
- 824 xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
- 825 ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
- 826 result = least_squares(xs, ys, function, silent=silent, **kwargs)
- 827 return result
- 828
- 829 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
- 830 """ Extract a plateau value from a Corr object
- 831
- 832 Parameters
- 833 ----------
- 834 plateau_range : list
- 835 list with two entries, indicating the first and the last timeslice
- 836 of the plateau region.
- 837 method : str
- 838 method to extract the plateau.
- 839 'fit' fits a constant to the plateau region
- 840 'avg', 'average' or 'mean' just average over the given timeslices.
- 841 auto_gamma : bool
- 842 apply gamma_method with default parameters to the Corr. Defaults to None
- 843 """
- 844 if not plateau_range:
- 845 if self.prange:
- 846 plateau_range = self.prange
- 847 else:
- 848 raise ValueError("no plateau range provided")
- 849 if self.N != 1:
- 850 raise ValueError("Correlator must be projected before getting a plateau.")
- 851 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
- 852 raise ValueError("plateau is undefined at all timeslices in plateaurange.")
- 853 if auto_gamma:
- 854 self.gamma_method()
- 855 if method == "fit":
- 856 def const_func(a, t):
- 857 return a[0]
- 858 return self.fit(const_func, plateau_range)[0]
- 859 elif method in ["avg", "average", "mean"]:
- 860 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
- 861 return returnvalue
- 862
- 863 else:
- 864 raise ValueError("Unsupported plateau method: " + method)
- 865
- 866 def set_prange(self, prange):
- 867 """Sets the attribute prange of the Corr object."""
- 868 if not len(prange) == 2:
- 869 raise ValueError("prange must be a list or array with two values")
- 870 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
- 871 raise TypeError("Start and end point must be integers")
- 872 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] <= prange[1]):
- 873 raise ValueError("Start and end point must define a range in the interval 0,T")
- 874
- 875 self.prange = prange
- 876 return
- 877
- 878 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
- 879 """Plots the correlator using the tag of the correlator as label if available.
- 880
- 881 Parameters
- 882 ----------
- 883 x_range : list
- 884 list of two values, determining the range of the x-axis e.g. [4, 8].
- 885 comp : Corr or list of Corr
- 886 Correlator or list of correlators which are plotted for comparison.
- 887 The tags of these correlators are used as labels if available.
- 888 logscale : bool
- 889 Sets y-axis to logscale.
- 890 plateau : Obs
- 891 Plateau value to be visualized in the figure.
- 892 fit_res : Fit_result
- 893 Fit_result object to be visualized.
- 894 fit_key : str
- 895 Key for the fit function in Fit_result.fit_function (for combined fits).
- 896 ylabel : str
- 897 Label for the y-axis.
- 898 save : str
- 899 path to file in which the figure should be saved.
- 900 auto_gamma : bool
- 901 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
- 902 hide_sigma : float
- 903 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
- 904 references : list
- 905 List of floating point values that are displayed as horizontal lines for reference.
- 906 title : string
- 907 Optional title of the figure.
- 908 """
- 909 if self.N != 1:
- 910 raise ValueError("Correlator must be projected before plotting")
- 911
- 912 if auto_gamma:
- 913 self.gamma_method()
- 914
- 915 if x_range is None:
- 916 x_range = [0, self.T - 1]
- 917
- 918 fig = plt.figure()
- 919 ax1 = fig.add_subplot(111)
- 920
- 921 x, y, y_err = self.plottable()
- 922 if hide_sigma:
- 923 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+ 740 elif variant == 'logsym':
+ 741 newcontent = []
+ 742 for t in range(1, self.T - 1):
+ 743 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0) or self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
+ 744 newcontent.append(None)
+ 745 else:
+ 746 newcontent.append(self.content[t - 1] / self.content[t + 1])
+ 747 if (all([x is None for x in newcontent])):
+ 748 raise ValueError('m_eff is undefined at all timeslices')
+ 749
+ 750 return np.log(Corr(newcontent, padding=[1, 1])) / 2
+ 751
+ 752 elif variant in ['periodic', 'cosh', 'sinh']:
+ 753 if variant in ['periodic', 'cosh']:
+ 754 func = anp.cosh
+ 755 else:
+ 756 func = anp.sinh
+ 757
+ 758 def root_function(x, d):
+ 759 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
+ 760
+ 761 newcontent = []
+ 762 for t in range(self.T - 1):
+ 763 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
+ 764 newcontent.append(None)
+ 765 # Fill the two timeslices in the middle of the lattice with their predecessors
+ 766 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+ 767 newcontent.append(newcontent[-1])
+ 768 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
+ 769 newcontent.append(None)
+ 770 else:
+ 771 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+ 772 if (all([x is None for x in newcontent])):
+ 773 raise ValueError('m_eff is undefined at all timeslices')
+ 774
+ 775 return Corr(newcontent, padding=[0, 1])
+ 776
+ 777 elif variant == 'arccosh':
+ 778 newcontent = []
+ 779 for t in range(1, self.T - 1):
+ 780 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0):
+ 781 newcontent.append(None)
+ 782 else:
+ 783 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+ 784 if (all([x is None for x in newcontent])):
+ 785 raise ValueError("m_eff is undefined at all timeslices")
+ 786 return np.arccosh(Corr(newcontent, padding=[1, 1]))
+ 787
+ 788 else:
+ 789 raise ValueError('Unknown variant.')
+ 790
+ 791 def fit(self, function, fitrange=None, silent=False, **kwargs):
+ 792 r'''Fits function to the data
+ 793
+ 794 Parameters
+ 795 ----------
+ 796 function : obj
+ 797 function to fit to the data. See fits.least_squares for details.
+ 798 fitrange : list
+ 799 Two element list containing the timeslices on which the fit is supposed to start and stop.
+ 800 Caution: This range is inclusive as opposed to standard python indexing.
+ 801 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+ 802 If not specified, self.prange or all timeslices are used.
+ 803 silent : bool
+ 804 Decides whether output is printed to the standard output.
+ 805 '''
+ 806 if self.N != 1:
+ 807 raise ValueError("Correlator must be projected before fitting")
+ 808
+ 809 if fitrange is None:
+ 810 if self.prange:
+ 811 fitrange = self.prange
+ 812 else:
+ 813 fitrange = [0, self.T - 1]
+ 814 else:
+ 815 if not isinstance(fitrange, list):
+ 816 raise TypeError("fitrange has to be a list with two elements")
+ 817 if len(fitrange) != 2:
+ 818 raise ValueError("fitrange has to have exactly two elements [fit_start, fit_stop]")
+ 819
+ 820 xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
+ 821 ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
+ 822 result = least_squares(xs, ys, function, silent=silent, **kwargs)
+ 823 return result
+ 824
+ 825 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+ 826 """ Extract a plateau value from a Corr object
+ 827
+ 828 Parameters
+ 829 ----------
+ 830 plateau_range : list
+ 831 list with two entries, indicating the first and the last timeslice
+ 832 of the plateau region.
+ 833 method : str
+ 834 method to extract the plateau.
+ 835 'fit' fits a constant to the plateau region
+ 836 'avg', 'average' or 'mean' just average over the given timeslices.
+ 837 auto_gamma : bool
+ 838 apply gamma_method with default parameters to the Corr. Defaults to None
+ 839 """
+ 840 if not plateau_range:
+ 841 if self.prange:
+ 842 plateau_range = self.prange
+ 843 else:
+ 844 raise ValueError("no plateau range provided")
+ 845 if self.N != 1:
+ 846 raise ValueError("Correlator must be projected before getting a plateau.")
+ 847 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+ 848 raise ValueError("plateau is undefined at all timeslices in plateaurange.")
+ 849 if auto_gamma:
+ 850 self.gamma_method()
+ 851 if method == "fit":
+ 852 def const_func(a, t):
+ 853 return a[0]
+ 854 return self.fit(const_func, plateau_range)[0]
+ 855 elif method in ["avg", "average", "mean"]:
+ 856 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+ 857 return returnvalue
+ 858
+ 859 else:
+ 860 raise ValueError("Unsupported plateau method: " + method)
+ 861
+ 862 def set_prange(self, prange):
+ 863 """Sets the attribute prange of the Corr object."""
+ 864 if not len(prange) == 2:
+ 865 raise ValueError("prange must be a list or array with two values")
+ 866 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+ 867 raise TypeError("Start and end point must be integers")
+ 868 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] <= prange[1]):
+ 869 raise ValueError("Start and end point must define a range in the interval 0,T")
+ 870
+ 871 self.prange = prange
+ 872
+ 873 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
+ 874 """Plots the correlator using the tag of the correlator as label if available.
+ 875
+ 876 Parameters
+ 877 ----------
+ 878 x_range : list
+ 879 list of two values, determining the range of the x-axis e.g. [4, 8].
+ 880 comp : Corr or list of Corr
+ 881 Correlator or list of correlators which are plotted for comparison.
+ 882 The tags of these correlators are used as labels if available.
+ 883 logscale : bool
+ 884 Sets y-axis to logscale.
+ 885 plateau : Obs
+ 886 Plateau value to be visualized in the figure.
+ 887 fit_res : Fit_result
+ 888 Fit_result object to be visualized.
+ 889 fit_key : str
+ 890 Key for the fit function in Fit_result.fit_function (for combined fits).
+ 891 ylabel : str
+ 892 Label for the y-axis.
+ 893 save : str
+ 894 path to file in which the figure should be saved.
+ 895 auto_gamma : bool
+ 896 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+ 897 hide_sigma : float
+ 898 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+ 899 references : list
+ 900 List of floating point values that are displayed as horizontal lines for reference.
+ 901 title : string
+ 902 Optional title of the figure.
+ 903 """
+ 904 if self.N != 1:
+ 905 raise ValueError("Correlator must be projected before plotting")
+ 906
+ 907 if auto_gamma:
+ 908 self.gamma_method()
+ 909
+ 910 if x_range is None:
+ 911 x_range = [0, self.T - 1]
+ 912
+ 913 fig = plt.figure()
+ 914 ax1 = fig.add_subplot(111)
+ 915
+ 916 x, y, y_err = self.plottable()
+ 917 if hide_sigma:
+ 918 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+ 919 else:
+ 920 hide_from = None
+ 921 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+ 922 if logscale:
+ 923 ax1.set_yscale('log')
924 else:
- 925 hide_from = None
- 926 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
- 927 if logscale:
- 928 ax1.set_yscale('log')
- 929 else:
- 930 if y_range is None:
- 931 try:
- 932 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
- 933 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
- 934 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
- 935 except Exception:
- 936 pass
- 937 else:
- 938 ax1.set_ylim(y_range)
- 939 if comp:
- 940 if isinstance(comp, (Corr, list)):
- 941 for corr in comp if isinstance(comp, list) else [comp]:
- 942 if auto_gamma:
- 943 corr.gamma_method()
- 944 x, y, y_err = corr.plottable()
- 945 if hide_sigma:
- 946 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
- 947 else:
- 948 hide_from = None
- 949 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
- 950 else:
- 951 raise TypeError("'comp' must be a correlator or a list of correlators.")
- 952
- 953 if plateau:
- 954 if isinstance(plateau, Obs):
- 955 if auto_gamma:
- 956 plateau.gamma_method()
- 957 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
- 958 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
- 959 else:
- 960 raise TypeError("'plateau' must be an Obs")
- 961
- 962 if references:
- 963 if isinstance(references, list):
- 964 for ref in references:
- 965 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
- 966 else:
- 967 raise TypeError("'references' must be a list of floating pint values.")
- 968
- 969 if self.prange:
- 970 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
- 971 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
- 972
- 973 if fit_res:
- 974 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
- 975 if isinstance(fit_res.fit_function, dict):
- 976 if fit_key:
- 977 ax1.plot(x_samples, fit_res.fit_function[fit_key]([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
- 978 else:
- 979 raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
- 980 else:
- 981 ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
+ 925 if y_range is None:
+ 926 try:
+ 927 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+ 928 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+ 929 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+ 930 except Exception:
+ 931 pass
+ 932 else:
+ 933 ax1.set_ylim(y_range)
+ 934 if comp:
+ 935 if isinstance(comp, (Corr, list)):
+ 936 for corr in comp if isinstance(comp, list) else [comp]:
+ 937 if auto_gamma:
+ 938 corr.gamma_method()
+ 939 x, y, y_err = corr.plottable()
+ 940 if hide_sigma:
+ 941 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+ 942 else:
+ 943 hide_from = None
+ 944 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+ 945 else:
+ 946 raise TypeError("'comp' must be a correlator or a list of correlators.")
+ 947
+ 948 if plateau:
+ 949 if isinstance(plateau, Obs):
+ 950 if auto_gamma:
+ 951 plateau.gamma_method()
+ 952 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+ 953 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+ 954 else:
+ 955 raise TypeError("'plateau' must be an Obs")
+ 956
+ 957 if references:
+ 958 if isinstance(references, list):
+ 959 for ref in references:
+ 960 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+ 961 else:
+ 962 raise TypeError("'references' must be a list of floating pint values.")
+ 963
+ 964 if self.prange:
+ 965 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
+ 966 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
+ 967
+ 968 if fit_res:
+ 969 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+ 970 if isinstance(fit_res.fit_function, dict):
+ 971 if fit_key:
+ 972 ax1.plot(x_samples, fit_res.fit_function[fit_key]([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
+ 973 else:
+ 974 raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
+ 975 else:
+ 976 ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
+ 977
+ 978 ax1.set_xlabel(r'$x_0 / a$')
+ 979 if ylabel:
+ 980 ax1.set_ylabel(ylabel)
+ 981 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
982
- 983 ax1.set_xlabel(r'$x_0 / a$')
- 984 if ylabel:
- 985 ax1.set_ylabel(ylabel)
- 986 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
- 987
- 988 _handles, labels = ax1.get_legend_handles_labels()
- 989 if labels:
- 990 ax1.legend()
+ 983 _handles, labels = ax1.get_legend_handles_labels()
+ 984 if labels:
+ 985 ax1.legend()
+ 986
+ 987 if title:
+ 988 plt.title(title)
+ 989
+ 990 plt.draw()
991
- 992 if title:
- 993 plt.title(title)
- 994
- 995 plt.draw()
- 996
- 997 if save:
- 998 if isinstance(save, str):
- 999 fig.savefig(save, bbox_inches='tight')
-1000 else:
-1001 raise TypeError("'save' has to be a string.")
-1002
-1003 def spaghetti_plot(self, logscale=True):
-1004 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-1005
-1006 Parameters
-1007 ----------
-1008 logscale : bool
-1009 Determines whether the scale of the y-axis is logarithmic or standard.
-1010 """
-1011 if self.N != 1:
-1012 raise ValueError("Correlator needs to be projected first.")
-1013
-1014 mc_names = list(set([item for sublist in [list(itertools.chain.from_iterable(map(o[0].e_content.get, o[0].mc_names))) for o in self.content if o is not None] for item in sublist]))
-1015 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content, strict=True) if o is not None]
-1016
-1017 for name in mc_names:
-1018 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+ 992 if save:
+ 993 if isinstance(save, str):
+ 994 fig.savefig(save, bbox_inches='tight')
+ 995 else:
+ 996 raise TypeError("'save' has to be a string.")
+ 997
+ 998 def spaghetti_plot(self, logscale=True):
+ 999 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+1000
+1001 Parameters
+1002 ----------
+1003 logscale : bool
+1004 Determines whether the scale of the y-axis is logarithmic or standard.
+1005 """
+1006 if self.N != 1:
+1007 raise ValueError("Correlator needs to be projected first.")
+1008
+1009 mc_names = list(set([item for sublist in [list(itertools.chain.from_iterable(map(o[0].e_content.get, o[0].mc_names))) for o in self.content if o is not None] for item in sublist]))
+1010 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content, strict=True) if o is not None]
+1011
+1012 for name in mc_names:
+1013 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+1014
+1015 fig = plt.figure()
+1016 ax = fig.add_subplot(111)
+1017 for dat in data:
+1018 ax.plot(x0_vals, dat, ls='-', marker='')
1019
-1020 fig = plt.figure()
-1021 ax = fig.add_subplot(111)
-1022 for dat in data:
-1023 ax.plot(x0_vals, dat, ls='-', marker='')
-1024
-1025 if logscale is True:
-1026 ax.set_yscale('log')
-1027
-1028 ax.set_xlabel(r'$x_0 / a$')
-1029 plt.title(name)
-1030 plt.draw()
-1031
-1032 def dump(self, filename, datatype="json.gz", **kwargs):
-1033 """Dumps the Corr into a file of chosen type
-1034 Parameters
-1035 ----------
-1036 filename : str
-1037 Name of the file to be saved.
-1038 datatype : str
-1039 Format of the exported file. Supported formats include
-1040 "json.gz" and "pickle"
-1041 path : str
-1042 specifies a custom path for the file (default '.')
-1043 """
-1044 if datatype == "json.gz":
-1045 from .input.json import dump_to_json
-1046 if 'path' in kwargs:
-1047 file_name = kwargs.get('path') + '/' + filename
-1048 else:
-1049 file_name = filename
-1050 dump_to_json(self, file_name)
-1051 elif datatype == "pickle":
-1052 dump_object(self, filename, **kwargs)
-1053 else:
-1054 raise ValueError("Unknown datatype " + str(datatype))
-1055
-1056 def print(self, print_range=None):
-1057 print(self.__repr__(print_range))
-1058
-1059 def __repr__(self, print_range=None):
-1060 if print_range is None:
-1061 print_range = [0, None]
-1062
-1063 content_string = ""
-1064 content_string += "Corr T=" + str(self.T) + " N=" + str(self.N) + "\n" # +" filled with"+ str(type(self.content[0][0])) there should be a good solution here
+1020 if logscale is True:
+1021 ax.set_yscale('log')
+1022
+1023 ax.set_xlabel(r'$x_0 / a$')
+1024 plt.title(name)
+1025 plt.draw()
+1026
+1027 def dump(self, filename, datatype="json.gz", **kwargs):
+1028 """Dumps the Corr into a file of chosen type
+1029 Parameters
+1030 ----------
+1031 filename : str
+1032 Name of the file to be saved.
+1033 datatype : str
+1034 Format of the exported file. Supported formats include
+1035 "json.gz" and "pickle"
+1036 path : str
+1037 specifies a custom path for the file (default '.')
+1038 """
+1039 if datatype == "json.gz":
+1040 from .input.json import dump_to_json
+1041 if 'path' in kwargs:
+1042 file_name = kwargs.get('path') + '/' + filename
+1043 else:
+1044 file_name = filename
+1045 dump_to_json(self, file_name)
+1046 elif datatype == "pickle":
+1047 dump_object(self, filename, **kwargs)
+1048 else:
+1049 raise ValueError("Unknown datatype " + str(datatype))
+1050
+1051 def print(self, print_range=None):
+1052 print(self.__repr__(print_range))
+1053
+1054 def __repr__(self, print_range=None):
+1055 if print_range is None:
+1056 print_range = [0, None]
+1057
+1058 content_string = ""
+1059 content_string += "Corr T=" + str(self.T) + " N=" + str(self.N) + "\n" # +" filled with"+ str(type(self.content[0][0])) there should be a good solution here
+1060
+1061 if self.tag is not None:
+1062 content_string += "Description: " + self.tag + "\n"
+1063 if self.N != 1:
+1064 return content_string
1065
-1066 if self.tag is not None:
-1067 content_string += "Description: " + self.tag + "\n"
-1068 if self.N != 1:
-1069 return content_string
-1070
-1071 if print_range[1]:
-1072 print_range[1] += 1
-1073 content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
-1074 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]):
-1075 if sub_corr is None:
-1076 content_string += str(i + print_range[0]) + '\n'
-1077 else:
-1078 content_string += str(i + print_range[0])
-1079 for element in sub_corr:
-1080 content_string += f"\t{element:+2}"
-1081 content_string += '\n'
-1082 return content_string
-1083
-1084 def __str__(self):
-1085 return self.__repr__()
+1066 if print_range[1]:
+1067 print_range[1] += 1
+1068 content_string += 'x0/a\tCorr(x0/a)\n------------------\n'
+1069 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]):
+1070 if sub_corr is None:
+1071 content_string += str(i + print_range[0]) + '\n'
+1072 else:
+1073 content_string += str(i + print_range[0])
+1074 for element in sub_corr:
+1075 content_string += f"\t{element:+2}"
+1076 content_string += '\n'
+1077 return content_string
+1078
+1079 def __str__(self):
+1080 return self.__repr__()
+1081
+1082 # We define the basic operations, that can be performed with correlators.
+1083 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
+1084 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
+1085 # One could try and tell Obs to check if the y in __mul__ is a Corr and
1086
-1087 # We define the basic operations, that can be performed with correlators.
-1088 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr.
-1089 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception.
-1090 # One could try and tell Obs to check if the y in __mul__ is a Corr and
-1091
-1092 __array_priority__ = 10000
-1093
-1094 def __eq__(self, y):
-1095 if isinstance(y, Corr):
-1096 comp = np.asarray(y.content, dtype=object)
-1097 else:
-1098 comp = np.asarray(y)
-1099 return np.asarray(self.content, dtype=object) == comp
-1100
-1101 __hash__ = None
-1102
-1103 def __add__(self, y):
-1104 if isinstance(y, Corr):
-1105 if ((self.N != y.N) or (self.T != y.T)):
-1106 raise ValueError("Addition of Corrs with different shape")
-1107 newcontent = []
-1108 for t in range(self.T):
-1109 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
-1110 newcontent.append(None)
-1111 else:
-1112 newcontent.append(self.content[t] + y.content[t])
-1113 return Corr(newcontent)
-1114
-1115 elif isinstance(y, (Obs, int, float, CObs, complex)):
-1116 newcontent = []
-1117 for t in range(self.T):
-1118 if _check_for_none(self, self.content[t]):
-1119 newcontent.append(None)
-1120 else:
-1121 newcontent.append(self.content[t] + y)
-1122 return Corr(newcontent, prange=self.prange)
-1123 elif isinstance(y, np.ndarray):
-1124 if y.shape == (self.T,):
-1125 return Corr(list((np.array(self.content).T + y).T))
-1126 else:
-1127 raise ValueError("operands could not be broadcast together")
-1128 else:
-1129 raise TypeError("Corr + wrong type")
-1130
-1131 def __mul__(self, y):
-1132 if isinstance(y, Corr):
-1133 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
-1134 raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
-1135 newcontent = []
-1136 for t in range(self.T):
-1137 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
-1138 newcontent.append(None)
-1139 else:
-1140 newcontent.append(self.content[t] * y.content[t])
-1141 return Corr(newcontent)
-1142
-1143 elif isinstance(y, (Obs, int, float, CObs, complex)):
-1144 newcontent = []
-1145 for t in range(self.T):
-1146 if _check_for_none(self, self.content[t]):
-1147 newcontent.append(None)
-1148 else:
-1149 newcontent.append(self.content[t] * y)
-1150 return Corr(newcontent, prange=self.prange)
-1151 elif isinstance(y, np.ndarray):
-1152 if y.shape == (self.T,):
-1153 return Corr(list((np.array(self.content).T * y).T))
-1154 else:
-1155 raise ValueError("operands could not be broadcast together")
-1156 else:
-1157 raise TypeError("Corr * wrong type")
-1158
-1159 def __matmul__(self, y):
-1160 if isinstance(y, np.ndarray):
-1161 if y.ndim != 2 or y.shape[0] != y.shape[1]:
-1162 raise ValueError("Can only multiply correlators by square matrices.")
-1163 if not self.N == y.shape[0]:
-1164 raise ValueError("matmul: mismatch of matrix dimensions")
-1165 newcontent = []
-1166 for t in range(self.T):
-1167 if _check_for_none(self, self.content[t]):
-1168 newcontent.append(None)
-1169 else:
-1170 newcontent.append(self.content[t] @ y)
-1171 return Corr(newcontent)
-1172 elif isinstance(y, Corr):
-1173 if not self.N == y.N:
-1174 raise ValueError("matmul: mismatch of matrix dimensions")
-1175 newcontent = []
-1176 for t in range(self.T):
-1177 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
-1178 newcontent.append(None)
-1179 else:
-1180 newcontent.append(self.content[t] @ y.content[t])
-1181 return Corr(newcontent)
-1182
-1183 else:
-1184 return NotImplemented
-1185
-1186 def __rmatmul__(self, y):
-1187 if isinstance(y, np.ndarray):
-1188 if y.ndim != 2 or y.shape[0] != y.shape[1]:
-1189 raise ValueError("Can only multiply correlators by square matrices.")
-1190 if not self.N == y.shape[0]:
-1191 raise ValueError("matmul: mismatch of matrix dimensions")
-1192 newcontent = []
-1193 for t in range(self.T):
-1194 if _check_for_none(self, self.content[t]):
-1195 newcontent.append(None)
-1196 else:
-1197 newcontent.append(y @ self.content[t])
-1198 return Corr(newcontent)
-1199 else:
-1200 return NotImplemented
-1201
-1202 def __truediv__(self, y):
-1203 if isinstance(y, Corr):
-1204 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
-1205 raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
-1206 newcontent = []
+1087 __array_priority__ = 10000
+1088
+1089 def __eq__(self, y):
+1090 if isinstance(y, Corr):
+1091 comp = np.asarray(y.content, dtype=object)
+1092 else:
+1093 comp = np.asarray(y)
+1094 return np.asarray(self.content, dtype=object) == comp
+1095
+1096 __hash__ = None
+1097
+1098 def __add__(self, y):
+1099 if isinstance(y, Corr):
+1100 if ((self.N != y.N) or (self.T != y.T)):
+1101 raise ValueError("Addition of Corrs with different shape")
+1102 newcontent = []
+1103 for t in range(self.T):
+1104 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
+1105 newcontent.append(None)
+1106 else:
+1107 newcontent.append(self.content[t] + y.content[t])
+1108 return Corr(newcontent)
+1109
+1110 elif isinstance(y, (Obs, int, float, CObs, complex)):
+1111 newcontent = []
+1112 for t in range(self.T):
+1113 if _check_for_none(self, self.content[t]):
+1114 newcontent.append(None)
+1115 else:
+1116 newcontent.append(self.content[t] + y)
+1117 return Corr(newcontent, prange=self.prange)
+1118 elif isinstance(y, np.ndarray):
+1119 if y.shape == (self.T,):
+1120 return Corr(list((np.array(self.content).T + y).T))
+1121 else:
+1122 raise ValueError("operands could not be broadcast together")
+1123 else:
+1124 raise TypeError("Corr + wrong type")
+1125
+1126 def __mul__(self, y):
+1127 if isinstance(y, Corr):
+1128 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
+1129 raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
+1130 newcontent = []
+1131 for t in range(self.T):
+1132 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
+1133 newcontent.append(None)
+1134 else:
+1135 newcontent.append(self.content[t] * y.content[t])
+1136 return Corr(newcontent)
+1137
+1138 elif isinstance(y, (Obs, int, float, CObs, complex)):
+1139 newcontent = []
+1140 for t in range(self.T):
+1141 if _check_for_none(self, self.content[t]):
+1142 newcontent.append(None)
+1143 else:
+1144 newcontent.append(self.content[t] * y)
+1145 return Corr(newcontent, prange=self.prange)
+1146 elif isinstance(y, np.ndarray):
+1147 if y.shape == (self.T,):
+1148 return Corr(list((np.array(self.content).T * y).T))
+1149 else:
+1150 raise ValueError("operands could not be broadcast together")
+1151 else:
+1152 raise TypeError("Corr * wrong type")
+1153
+1154 def __matmul__(self, y):
+1155 if isinstance(y, np.ndarray):
+1156 if y.ndim != 2 or y.shape[0] != y.shape[1]:
+1157 raise ValueError("Can only multiply correlators by square matrices.")
+1158 if not self.N == y.shape[0]:
+1159 raise ValueError("matmul: mismatch of matrix dimensions")
+1160 newcontent = []
+1161 for t in range(self.T):
+1162 if _check_for_none(self, self.content[t]):
+1163 newcontent.append(None)
+1164 else:
+1165 newcontent.append(self.content[t] @ y)
+1166 return Corr(newcontent)
+1167 elif isinstance(y, Corr):
+1168 if not self.N == y.N:
+1169 raise ValueError("matmul: mismatch of matrix dimensions")
+1170 newcontent = []
+1171 for t in range(self.T):
+1172 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
+1173 newcontent.append(None)
+1174 else:
+1175 newcontent.append(self.content[t] @ y.content[t])
+1176 return Corr(newcontent)
+1177
+1178 else:
+1179 return NotImplemented
+1180
+1181 def __rmatmul__(self, y):
+1182 if isinstance(y, np.ndarray):
+1183 if y.ndim != 2 or y.shape[0] != y.shape[1]:
+1184 raise ValueError("Can only multiply correlators by square matrices.")
+1185 if not self.N == y.shape[0]:
+1186 raise ValueError("matmul: mismatch of matrix dimensions")
+1187 newcontent = []
+1188 for t in range(self.T):
+1189 if _check_for_none(self, self.content[t]):
+1190 newcontent.append(None)
+1191 else:
+1192 newcontent.append(y @ self.content[t])
+1193 return Corr(newcontent)
+1194 else:
+1195 return NotImplemented
+1196
+1197 def __truediv__(self, y):
+1198 if isinstance(y, Corr):
+1199 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T):
+1200 raise ValueError("Multiplication of Corr object requires N=N or N=1 and T=T")
+1201 newcontent = []
+1202 for t in range(self.T):
+1203 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
+1204 newcontent.append(None)
+1205 else:
+1206 newcontent.append(self.content[t] / y.content[t])
1207 for t in range(self.T):
-1208 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]):
-1209 newcontent.append(None)
-1210 else:
-1211 newcontent.append(self.content[t] / y.content[t])
-1212 for t in range(self.T):
-1213 if _check_for_none(self, newcontent[t]):
-1214 continue
-1215 if np.isnan(np.sum(newcontent[t]).value):
-1216 newcontent[t] = None
-1217
-1218 if all([item is None for item in newcontent]):
-1219 raise ValueError("Division returns completely undefined correlator")
-1220 return Corr(newcontent)
-1221
-1222 elif isinstance(y, (Obs, CObs)):
-1223 if isinstance(y, Obs):
-1224 if y.value == 0:
-1225 raise ValueError('Division by zero will return undefined correlator')
-1226 if isinstance(y, CObs):
-1227 if y.is_zero():
-1228 raise ValueError('Division by zero will return undefined correlator')
-1229
-1230 newcontent = []
-1231 for t in range(self.T):
-1232 if _check_for_none(self, self.content[t]):
-1233 newcontent.append(None)
-1234 else:
-1235 newcontent.append(self.content[t] / y)
-1236 return Corr(newcontent, prange=self.prange)
-1237
-1238 elif isinstance(y, (int, float)):
-1239 if y == 0:
-1240 raise ValueError('Division by zero will return undefined correlator')
-1241 newcontent = []
-1242 for t in range(self.T):
-1243 if _check_for_none(self, self.content[t]):
-1244 newcontent.append(None)
-1245 else:
-1246 newcontent.append(self.content[t] / y)
-1247 return Corr(newcontent, prange=self.prange)
-1248 elif isinstance(y, np.ndarray):
-1249 if y.shape == (self.T,):
-1250 return Corr(list((np.array(self.content).T / y).T))
-1251 else:
-1252 raise ValueError("operands could not be broadcast together")
-1253 else:
-1254 raise TypeError('Corr / wrong type')
-1255
-1256 def __neg__(self):
-1257 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
-1258 return Corr(newcontent, prange=self.prange)
-1259
-1260 def __sub__(self, y):
-1261 return self + (-y)
-1262
-1263 def __pow__(self, y):
-1264 if isinstance(y, (Obs, int, float, CObs)):
-1265 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
-1266 return Corr(newcontent, prange=self.prange)
-1267 else:
-1268 raise TypeError('Type of exponent not supported')
-1269
-1270 def __abs__(self):
-1271 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
-1272 return Corr(newcontent, prange=self.prange)
-1273
-1274 # The numpy functions:
-1275 def sqrt(self):
-1276 return self ** 0.5
-1277
-1278 def log(self):
-1279 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
-1280 return Corr(newcontent, prange=self.prange)
-1281
-1282 def exp(self):
-1283 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
-1284 return Corr(newcontent, prange=self.prange)
-1285
-1286 def _apply_func_to_corr(self, func):
-1287 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
-1288 for t in range(self.T):
-1289 if _check_for_none(self, newcontent[t]):
-1290 continue
-1291 tmp_sum = np.sum(newcontent[t])
-1292 if hasattr(tmp_sum, "value"):
-1293 if np.isnan(tmp_sum.value):
-1294 newcontent[t] = None
-1295 if all([item is None for item in newcontent]):
-1296 raise ValueError('Operation returns undefined correlator')
-1297 return Corr(newcontent)
-1298
-1299 def sin(self):
-1300 return self._apply_func_to_corr(np.sin)
-1301
-1302 def cos(self):
-1303 return self._apply_func_to_corr(np.cos)
-1304
-1305 def tan(self):
-1306 return self._apply_func_to_corr(np.tan)
-1307
-1308 def sinh(self):
-1309 return self._apply_func_to_corr(np.sinh)
-1310
-1311 def cosh(self):
-1312 return self._apply_func_to_corr(np.cosh)
-1313
-1314 def tanh(self):
-1315 return self._apply_func_to_corr(np.tanh)
-1316
-1317 def arcsin(self):
-1318 return self._apply_func_to_corr(np.arcsin)
-1319
-1320 def arccos(self):
-1321 return self._apply_func_to_corr(np.arccos)
-1322
-1323 def arctan(self):
-1324 return self._apply_func_to_corr(np.arctan)
-1325
-1326 def arcsinh(self):
-1327 return self._apply_func_to_corr(np.arcsinh)
-1328
-1329 def arccosh(self):
-1330 return self._apply_func_to_corr(np.arccosh)
-1331
-1332 def arctanh(self):
-1333 return self._apply_func_to_corr(np.arctanh)
-1334
-1335 # Right hand side operations (require tweak in main module to work)
-1336 def __radd__(self, y):
-1337 return self + y
-1338
-1339 def __rsub__(self, y):
-1340 return -self + y
-1341
-1342 def __rmul__(self, y):
-1343 return self * y
-1344
-1345 def __rtruediv__(self, y):
-1346 return (self / y) ** (-1)
-1347
-1348 @property
-1349 def real(self):
-1350 def return_real(obs_OR_cobs):
-1351 if isinstance(obs_OR_cobs.flatten()[0], CObs):
-1352 return np.vectorize(lambda x: x.real)(obs_OR_cobs)
-1353 else:
-1354 return obs_OR_cobs
-1355
-1356 return self._apply_func_to_corr(return_real)
-1357
-1358 @property
-1359 def imag(self):
-1360 def return_imag(obs_OR_cobs):
-1361 if isinstance(obs_OR_cobs.flatten()[0], CObs):
-1362 return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
-1363 else:
-1364 return obs_OR_cobs * 0 # So it stays the right type
+1208 if _check_for_none(self, newcontent[t]):
+1209 continue
+1210 if np.isnan(np.sum(newcontent[t]).value):
+1211 newcontent[t] = None
+1212
+1213 if all([item is None for item in newcontent]):
+1214 raise ValueError("Division returns completely undefined correlator")
+1215 return Corr(newcontent)
+1216
+1217 elif isinstance(y, (Obs, CObs)):
+1218 if isinstance(y, Obs):
+1219 if y.value == 0:
+1220 raise ValueError('Division by zero will return undefined correlator')
+1221 if isinstance(y, CObs):
+1222 if y.is_zero():
+1223 raise ValueError('Division by zero will return undefined correlator')
+1224
+1225 newcontent = []
+1226 for t in range(self.T):
+1227 if _check_for_none(self, self.content[t]):
+1228 newcontent.append(None)
+1229 else:
+1230 newcontent.append(self.content[t] / y)
+1231 return Corr(newcontent, prange=self.prange)
+1232
+1233 elif isinstance(y, (int, float)):
+1234 if y == 0:
+1235 raise ValueError('Division by zero will return undefined correlator')
+1236 newcontent = []
+1237 for t in range(self.T):
+1238 if _check_for_none(self, self.content[t]):
+1239 newcontent.append(None)
+1240 else:
+1241 newcontent.append(self.content[t] / y)
+1242 return Corr(newcontent, prange=self.prange)
+1243 elif isinstance(y, np.ndarray):
+1244 if y.shape == (self.T,):
+1245 return Corr(list((np.array(self.content).T / y).T))
+1246 else:
+1247 raise ValueError("operands could not be broadcast together")
+1248 else:
+1249 raise TypeError('Corr / wrong type')
+1250
+1251 def __neg__(self):
+1252 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content]
+1253 return Corr(newcontent, prange=self.prange)
+1254
+1255 def __sub__(self, y):
+1256 return self + (-y)
+1257
+1258 def __pow__(self, y):
+1259 if isinstance(y, (Obs, int, float, CObs)):
+1260 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content]
+1261 return Corr(newcontent, prange=self.prange)
+1262 else:
+1263 raise TypeError('Type of exponent not supported')
+1264
+1265 def __abs__(self):
+1266 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content]
+1267 return Corr(newcontent, prange=self.prange)
+1268
+1269 # The numpy functions:
+1270 def sqrt(self):
+1271 return self ** 0.5
+1272
+1273 def log(self):
+1274 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
+1275 return Corr(newcontent, prange=self.prange)
+1276
+1277 def exp(self):
+1278 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
+1279 return Corr(newcontent, prange=self.prange)
+1280
+1281 def _apply_func_to_corr(self, func):
+1282 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content]
+1283 for t in range(self.T):
+1284 if _check_for_none(self, newcontent[t]):
+1285 continue
+1286 tmp_sum = np.sum(newcontent[t])
+1287 if hasattr(tmp_sum, "value"):
+1288 if np.isnan(tmp_sum.value):
+1289 newcontent[t] = None
+1290 if all([item is None for item in newcontent]):
+1291 raise ValueError('Operation returns undefined correlator')
+1292 return Corr(newcontent)
+1293
+1294 def sin(self):
+1295 return self._apply_func_to_corr(np.sin)
+1296
+1297 def cos(self):
+1298 return self._apply_func_to_corr(np.cos)
+1299
+1300 def tan(self):
+1301 return self._apply_func_to_corr(np.tan)
+1302
+1303 def sinh(self):
+1304 return self._apply_func_to_corr(np.sinh)
+1305
+1306 def cosh(self):
+1307 return self._apply_func_to_corr(np.cosh)
+1308
+1309 def tanh(self):
+1310 return self._apply_func_to_corr(np.tanh)
+1311
+1312 def arcsin(self):
+1313 return self._apply_func_to_corr(np.arcsin)
+1314
+1315 def arccos(self):
+1316 return self._apply_func_to_corr(np.arccos)
+1317
+1318 def arctan(self):
+1319 return self._apply_func_to_corr(np.arctan)
+1320
+1321 def arcsinh(self):
+1322 return self._apply_func_to_corr(np.arcsinh)
+1323
+1324 def arccosh(self):
+1325 return self._apply_func_to_corr(np.arccosh)
+1326
+1327 def arctanh(self):
+1328 return self._apply_func_to_corr(np.arctanh)
+1329
+1330 # Right hand side operations (require tweak in main module to work)
+1331 def __radd__(self, y):
+1332 return self + y
+1333
+1334 def __rsub__(self, y):
+1335 return -self + y
+1336
+1337 def __rmul__(self, y):
+1338 return self * y
+1339
+1340 def __rtruediv__(self, y):
+1341 return (self / y) ** (-1)
+1342
+1343 @property
+1344 def real(self):
+1345 def return_real(obs_OR_cobs):
+1346 if isinstance(obs_OR_cobs.flatten()[0], CObs):
+1347 return np.vectorize(lambda x: x.real)(obs_OR_cobs)
+1348 else:
+1349 return obs_OR_cobs
+1350
+1351 return self._apply_func_to_corr(return_real)
+1352
+1353 @property
+1354 def imag(self):
+1355 def return_imag(obs_OR_cobs):
+1356 if isinstance(obs_OR_cobs.flatten()[0], CObs):
+1357 return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
+1358 else:
+1359 return obs_OR_cobs * 0 # So it stays the right type
+1360
+1361 return self._apply_func_to_corr(return_imag)
+1362
+1363 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1364 r''' Project large correlation matrix to lowest states
1365
-1366 return self._apply_func_to_corr(return_imag)
-1367
-1368 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1369 r''' Project large correlation matrix to lowest states
-1370
-1371 This method can be used to reduce the size of an (N x N) correlation matrix
-1372 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1373 is still small.
-1374
-1375 Parameters
-1376 ----------
-1377 Ntrunc: int
-1378 Rank of the target matrix.
-1379 tproj: int
-1380 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1381 The default value is 3.
-1382 t0proj: int
-1383 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1384 discouraged for O(a) improved theories, since the correctness of the procedure
-1385 cannot be granted in this case. The default value is 2.
-1386 basematrix : Corr
-1387 Correlation matrix that is used to determine the eigenvectors of the
-1388 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1389 is is not specified.
-1390
-1391 Notes
-1392 -----
-1393 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1394 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
-1395 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1396 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1397 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1398 correlation matrix and to remove some noise that is added by irrelevant operators.
-1399 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1400 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1401 '''
-1402
-1403 if self.N == 1:
-1404 raise ValueError('Method cannot be applied to one-dimensional correlators.')
-1405 if basematrix is None:
-1406 basematrix = self
-1407 if Ntrunc >= basematrix.N:
-1408 raise ValueError(f'Cannot truncate using Ntrunc >= {basematrix.N}')
-1409 if basematrix.N != self.N:
-1410 raise ValueError('basematrix and targetmatrix have to be of the same size.')
-1411
-1412 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
-1413
-1414 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1415 rmat = []
-1416 for t in range(basematrix.T):
-1417 if self.content[t] is None:
-1418 rmat.append(None)
-1419 else:
-1420 for i in range(Ntrunc):
-1421 for j in range(Ntrunc):
-1422 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1423 rmat.append(np.copy(tmpmat))
-1424
-1425 return Corr(rmat)
+1366 This method can be used to reduce the size of an (N x N) correlation matrix
+1367 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1368 is still small.
+1369
+1370 Parameters
+1371 ----------
+1372 Ntrunc: int
+1373 Rank of the target matrix.
+1374 tproj: int
+1375 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1376 The default value is 3.
+1377 t0proj: int
+1378 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1379 discouraged for O(a) improved theories, since the correctness of the procedure
+1380 cannot be granted in this case. The default value is 2.
+1381 basematrix : Corr
+1382 Correlation matrix that is used to determine the eigenvectors of the
+1383 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1384 is is not specified.
+1385
+1386 Notes
+1387 -----
+1388 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1389 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
+1390 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1391 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1392 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1393 correlation matrix and to remove some noise that is added by irrelevant operators.
+1394 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1395 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1396 '''
+1397
+1398 if self.N == 1:
+1399 raise ValueError('Method cannot be applied to one-dimensional correlators.')
+1400 if basematrix is None:
+1401 basematrix = self
+1402 if Ntrunc >= basematrix.N:
+1403 raise ValueError(f'Cannot truncate using Ntrunc >= {basematrix.N}')
+1404 if basematrix.N != self.N:
+1405 raise ValueError('basematrix and targetmatrix have to be of the same size.')
+1406
+1407 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
+1408
+1409 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1410 rmat = []
+1411 for t in range(basematrix.T):
+1412 if self.content[t] is None:
+1413 rmat.append(None)
+1414 else:
+1415 for i in range(Ntrunc):
+1416 for j in range(Ntrunc):
+1417 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1418 rmat.append(np.copy(tmpmat))
+1419
+1420 return Corr(rmat)
@@ -4535,69 +4525,65 @@ Available choice:
728 if variant == 'log':
729 newcontent = []
730 for t in range(self.T - 1):
-731 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
+731 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0) or self.content[t][0].value / self.content[t + 1][0].value < 0:
732 newcontent.append(None)
-733 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
-734 newcontent.append(None)
-735 else:
-736 newcontent.append(self.content[t] / self.content[t + 1])
-737 if (all([x is None for x in newcontent])):
-738 raise ValueError('m_eff is undefined at all timeslices')
+733 else:
+734 newcontent.append(self.content[t] / self.content[t + 1])
+735 if (all([x is None for x in newcontent])):
+736 raise ValueError('m_eff is undefined at all timeslices')
+737
+738 return np.log(Corr(newcontent, padding=[0, 1]))
739
-740 return np.log(Corr(newcontent, padding=[0, 1]))
-741
-742 elif variant == 'logsym':
-743 newcontent = []
-744 for t in range(1, self.T - 1):
-745 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
-746 newcontent.append(None)
-747 elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
-748 newcontent.append(None)
-749 else:
-750 newcontent.append(self.content[t - 1] / self.content[t + 1])
-751 if (all([x is None for x in newcontent])):
-752 raise ValueError('m_eff is undefined at all timeslices')
-753
-754 return np.log(Corr(newcontent, padding=[1, 1])) / 2
-755
-756 elif variant in ['periodic', 'cosh', 'sinh']:
-757 if variant in ['periodic', 'cosh']:
-758 func = anp.cosh
-759 else:
-760 func = anp.sinh
-761
-762 def root_function(x, d):
-763 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
-764
-765 newcontent = []
-766 for t in range(self.T - 1):
-767 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
-768 newcontent.append(None)
-769 # Fill the two timeslices in the middle of the lattice with their predecessors
-770 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
-771 newcontent.append(newcontent[-1])
-772 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
-773 newcontent.append(None)
-774 else:
-775 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
-776 if (all([x is None for x in newcontent])):
-777 raise ValueError('m_eff is undefined at all timeslices')
-778
-779 return Corr(newcontent, padding=[0, 1])
-780
-781 elif variant == 'arccosh':
-782 newcontent = []
-783 for t in range(1, self.T - 1):
-784 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0):
-785 newcontent.append(None)
-786 else:
-787 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
-788 if (all([x is None for x in newcontent])):
-789 raise ValueError("m_eff is undefined at all timeslices")
-790 return np.arccosh(Corr(newcontent, padding=[1, 1]))
-791
-792 else:
-793 raise ValueError('Unknown variant.')
+740 elif variant == 'logsym':
+741 newcontent = []
+742 for t in range(1, self.T - 1):
+743 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0) or self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
+744 newcontent.append(None)
+745 else:
+746 newcontent.append(self.content[t - 1] / self.content[t + 1])
+747 if (all([x is None for x in newcontent])):
+748 raise ValueError('m_eff is undefined at all timeslices')
+749
+750 return np.log(Corr(newcontent, padding=[1, 1])) / 2
+751
+752 elif variant in ['periodic', 'cosh', 'sinh']:
+753 if variant in ['periodic', 'cosh']:
+754 func = anp.cosh
+755 else:
+756 func = anp.sinh
+757
+758 def root_function(x, d):
+759 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
+760
+761 newcontent = []
+762 for t in range(self.T - 1):
+763 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
+764 newcontent.append(None)
+765 # Fill the two timeslices in the middle of the lattice with their predecessors
+766 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+767 newcontent.append(newcontent[-1])
+768 elif self.content[t][0].value / self.content[t + 1][0].value < 0:
+769 newcontent.append(None)
+770 else:
+771 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+772 if (all([x is None for x in newcontent])):
+773 raise ValueError('m_eff is undefined at all timeslices')
+774
+775 return Corr(newcontent, padding=[0, 1])
+776
+777 elif variant == 'arccosh':
+778 newcontent = []
+779 for t in range(1, self.T - 1):
+780 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None) or (self.content[t][0].value == 0):
+781 newcontent.append(None)
+782 else:
+783 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+784 if (all([x is None for x in newcontent])):
+785 raise ValueError("m_eff is undefined at all timeslices")
+786 return np.arccosh(Corr(newcontent, padding=[1, 1]))
+787
+788 else:
+789 raise ValueError('Unknown variant.')
@@ -4631,39 +4617,39 @@ guess for the root finder, only relevant for the root variant
- 795 def fit(self, function, fitrange=None, silent=False, **kwargs):
-796 r'''Fits function to the data
-797
-798 Parameters
-799 ----------
-800 function : obj
-801 function to fit to the data. See fits.least_squares for details.
-802 fitrange : list
-803 Two element list containing the timeslices on which the fit is supposed to start and stop.
-804 Caution: This range is inclusive as opposed to standard python indexing.
-805 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
-806 If not specified, self.prange or all timeslices are used.
-807 silent : bool
-808 Decides whether output is printed to the standard output.
-809 '''
-810 if self.N != 1:
-811 raise ValueError("Correlator must be projected before fitting")
-812
-813 if fitrange is None:
-814 if self.prange:
-815 fitrange = self.prange
-816 else:
-817 fitrange = [0, self.T - 1]
-818 else:
-819 if not isinstance(fitrange, list):
-820 raise TypeError("fitrange has to be a list with two elements")
-821 if len(fitrange) != 2:
-822 raise ValueError("fitrange has to have exactly two elements [fit_start, fit_stop]")
-823
-824 xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
-825 ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
-826 result = least_squares(xs, ys, function, silent=silent, **kwargs)
-827 return result
+ 791 def fit(self, function, fitrange=None, silent=False, **kwargs):
+792 r'''Fits function to the data
+793
+794 Parameters
+795 ----------
+796 function : obj
+797 function to fit to the data. See fits.least_squares for details.
+798 fitrange : list
+799 Two element list containing the timeslices on which the fit is supposed to start and stop.
+800 Caution: This range is inclusive as opposed to standard python indexing.
+801 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+802 If not specified, self.prange or all timeslices are used.
+803 silent : bool
+804 Decides whether output is printed to the standard output.
+805 '''
+806 if self.N != 1:
+807 raise ValueError("Correlator must be projected before fitting")
+808
+809 if fitrange is None:
+810 if self.prange:
+811 fitrange = self.prange
+812 else:
+813 fitrange = [0, self.T - 1]
+814 else:
+815 if not isinstance(fitrange, list):
+816 raise TypeError("fitrange has to be a list with two elements")
+817 if len(fitrange) != 2:
+818 raise ValueError("fitrange has to have exactly two elements [fit_start, fit_stop]")
+819
+820 xs = np.array([x for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
+821 ys = np.array([self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if self.content[x] is not None])
+822 result = least_squares(xs, ys, function, silent=silent, **kwargs)
+823 return result
829 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
-830 """ Extract a plateau value from a Corr object
-831
-832 Parameters
-833 ----------
-834 plateau_range : list
-835 list with two entries, indicating the first and the last timeslice
-836 of the plateau region.
-837 method : str
-838 method to extract the plateau.
-839 'fit' fits a constant to the plateau region
-840 'avg', 'average' or 'mean' just average over the given timeslices.
-841 auto_gamma : bool
-842 apply gamma_method with default parameters to the Corr. Defaults to None
-843 """
-844 if not plateau_range:
-845 if self.prange:
-846 plateau_range = self.prange
-847 else:
-848 raise ValueError("no plateau range provided")
-849 if self.N != 1:
-850 raise ValueError("Correlator must be projected before getting a plateau.")
-851 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
-852 raise ValueError("plateau is undefined at all timeslices in plateaurange.")
-853 if auto_gamma:
-854 self.gamma_method()
-855 if method == "fit":
-856 def const_func(a, t):
-857 return a[0]
-858 return self.fit(const_func, plateau_range)[0]
-859 elif method in ["avg", "average", "mean"]:
-860 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
-861 return returnvalue
-862
-863 else:
-864 raise ValueError("Unsupported plateau method: " + method)
+ 825 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+826 """ Extract a plateau value from a Corr object
+827
+828 Parameters
+829 ----------
+830 plateau_range : list
+831 list with two entries, indicating the first and the last timeslice
+832 of the plateau region.
+833 method : str
+834 method to extract the plateau.
+835 'fit' fits a constant to the plateau region
+836 'avg', 'average' or 'mean' just average over the given timeslices.
+837 auto_gamma : bool
+838 apply gamma_method with default parameters to the Corr. Defaults to None
+839 """
+840 if not plateau_range:
+841 if self.prange:
+842 plateau_range = self.prange
+843 else:
+844 raise ValueError("no plateau range provided")
+845 if self.N != 1:
+846 raise ValueError("Correlator must be projected before getting a plateau.")
+847 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+848 raise ValueError("plateau is undefined at all timeslices in plateaurange.")
+849 if auto_gamma:
+850 self.gamma_method()
+851 if method == "fit":
+852 def const_func(a, t):
+853 return a[0]
+854 return self.fit(const_func, plateau_range)[0]
+855 elif method in ["avg", "average", "mean"]:
+856 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+857 return returnvalue
+858
+859 else:
+860 raise ValueError("Unsupported plateau method: " + method)
866 def set_prange(self, prange):
-867 """Sets the attribute prange of the Corr object."""
-868 if not len(prange) == 2:
-869 raise ValueError("prange must be a list or array with two values")
-870 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
-871 raise TypeError("Start and end point must be integers")
-872 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] <= prange[1]):
-873 raise ValueError("Start and end point must define a range in the interval 0,T")
-874
-875 self.prange = prange
-876 return
+ 862 def set_prange(self, prange):
+863 """Sets the attribute prange of the Corr object."""
+864 if not len(prange) == 2:
+865 raise ValueError("prange must be a list or array with two values")
+866 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+867 raise TypeError("Start and end point must be integers")
+868 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] <= prange[1]):
+869 raise ValueError("Start and end point must define a range in the interval 0,T")
+870
+871 self.prange = prange
878 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
- 879 """Plots the correlator using the tag of the correlator as label if available.
- 880
- 881 Parameters
- 882 ----------
- 883 x_range : list
- 884 list of two values, determining the range of the x-axis e.g. [4, 8].
- 885 comp : Corr or list of Corr
- 886 Correlator or list of correlators which are plotted for comparison.
- 887 The tags of these correlators are used as labels if available.
- 888 logscale : bool
- 889 Sets y-axis to logscale.
- 890 plateau : Obs
- 891 Plateau value to be visualized in the figure.
- 892 fit_res : Fit_result
- 893 Fit_result object to be visualized.
- 894 fit_key : str
- 895 Key for the fit function in Fit_result.fit_function (for combined fits).
- 896 ylabel : str
- 897 Label for the y-axis.
- 898 save : str
- 899 path to file in which the figure should be saved.
- 900 auto_gamma : bool
- 901 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
- 902 hide_sigma : float
- 903 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
- 904 references : list
- 905 List of floating point values that are displayed as horizontal lines for reference.
- 906 title : string
- 907 Optional title of the figure.
- 908 """
- 909 if self.N != 1:
- 910 raise ValueError("Correlator must be projected before plotting")
- 911
- 912 if auto_gamma:
- 913 self.gamma_method()
- 914
- 915 if x_range is None:
- 916 x_range = [0, self.T - 1]
- 917
- 918 fig = plt.figure()
- 919 ax1 = fig.add_subplot(111)
- 920
- 921 x, y, y_err = self.plottable()
- 922 if hide_sigma:
- 923 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
- 924 else:
- 925 hide_from = None
- 926 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
- 927 if logscale:
- 928 ax1.set_yscale('log')
- 929 else:
- 930 if y_range is None:
- 931 try:
- 932 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
- 933 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
- 934 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
- 935 except Exception:
- 936 pass
- 937 else:
- 938 ax1.set_ylim(y_range)
- 939 if comp:
- 940 if isinstance(comp, (Corr, list)):
- 941 for corr in comp if isinstance(comp, list) else [comp]:
- 942 if auto_gamma:
- 943 corr.gamma_method()
- 944 x, y, y_err = corr.plottable()
- 945 if hide_sigma:
- 946 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
- 947 else:
- 948 hide_from = None
- 949 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
- 950 else:
- 951 raise TypeError("'comp' must be a correlator or a list of correlators.")
- 952
- 953 if plateau:
- 954 if isinstance(plateau, Obs):
- 955 if auto_gamma:
- 956 plateau.gamma_method()
- 957 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
- 958 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
- 959 else:
- 960 raise TypeError("'plateau' must be an Obs")
- 961
- 962 if references:
- 963 if isinstance(references, list):
- 964 for ref in references:
- 965 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
- 966 else:
- 967 raise TypeError("'references' must be a list of floating pint values.")
- 968
- 969 if self.prange:
- 970 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
- 971 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
- 972
- 973 if fit_res:
- 974 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
- 975 if isinstance(fit_res.fit_function, dict):
- 976 if fit_key:
- 977 ax1.plot(x_samples, fit_res.fit_function[fit_key]([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
- 978 else:
- 979 raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
- 980 else:
- 981 ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
- 982
- 983 ax1.set_xlabel(r'$x_0 / a$')
- 984 if ylabel:
- 985 ax1.set_ylabel(ylabel)
- 986 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
- 987
- 988 _handles, labels = ax1.get_legend_handles_labels()
- 989 if labels:
- 990 ax1.legend()
- 991
- 992 if title:
- 993 plt.title(title)
- 994
- 995 plt.draw()
- 996
- 997 if save:
- 998 if isinstance(save, str):
- 999 fig.savefig(save, bbox_inches='tight')
-1000 else:
-1001 raise TypeError("'save' has to be a string.")
+ 873 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, fit_key=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
+874 """Plots the correlator using the tag of the correlator as label if available.
+875
+876 Parameters
+877 ----------
+878 x_range : list
+879 list of two values, determining the range of the x-axis e.g. [4, 8].
+880 comp : Corr or list of Corr
+881 Correlator or list of correlators which are plotted for comparison.
+882 The tags of these correlators are used as labels if available.
+883 logscale : bool
+884 Sets y-axis to logscale.
+885 plateau : Obs
+886 Plateau value to be visualized in the figure.
+887 fit_res : Fit_result
+888 Fit_result object to be visualized.
+889 fit_key : str
+890 Key for the fit function in Fit_result.fit_function (for combined fits).
+891 ylabel : str
+892 Label for the y-axis.
+893 save : str
+894 path to file in which the figure should be saved.
+895 auto_gamma : bool
+896 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+897 hide_sigma : float
+898 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+899 references : list
+900 List of floating point values that are displayed as horizontal lines for reference.
+901 title : string
+902 Optional title of the figure.
+903 """
+904 if self.N != 1:
+905 raise ValueError("Correlator must be projected before plotting")
+906
+907 if auto_gamma:
+908 self.gamma_method()
+909
+910 if x_range is None:
+911 x_range = [0, self.T - 1]
+912
+913 fig = plt.figure()
+914 ax1 = fig.add_subplot(111)
+915
+916 x, y, y_err = self.plottable()
+917 if hide_sigma:
+918 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+919 else:
+920 hide_from = None
+921 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+922 if logscale:
+923 ax1.set_yscale('log')
+924 else:
+925 if y_range is None:
+926 try:
+927 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+928 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+929 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+930 except Exception:
+931 pass
+932 else:
+933 ax1.set_ylim(y_range)
+934 if comp:
+935 if isinstance(comp, (Corr, list)):
+936 for corr in comp if isinstance(comp, list) else [comp]:
+937 if auto_gamma:
+938 corr.gamma_method()
+939 x, y, y_err = corr.plottable()
+940 if hide_sigma:
+941 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+942 else:
+943 hide_from = None
+944 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+945 else:
+946 raise TypeError("'comp' must be a correlator or a list of correlators.")
+947
+948 if plateau:
+949 if isinstance(plateau, Obs):
+950 if auto_gamma:
+951 plateau.gamma_method()
+952 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+953 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+954 else:
+955 raise TypeError("'plateau' must be an Obs")
+956
+957 if references:
+958 if isinstance(references, list):
+959 for ref in references:
+960 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+961 else:
+962 raise TypeError("'references' must be a list of floating pint values.")
+963
+964 if self.prange:
+965 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',', color="black", zorder=0)
+966 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',', color="black", zorder=0)
+967
+968 if fit_res:
+969 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+970 if isinstance(fit_res.fit_function, dict):
+971 if fit_key:
+972 ax1.plot(x_samples, fit_res.fit_function[fit_key]([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
+973 else:
+974 raise ValueError("Please provide a 'fit_key' for visualizing combined fits.")
+975 else:
+976 ax1.plot(x_samples, fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), ls='-', marker=',', lw=2)
+977
+978 ax1.set_xlabel(r'$x_0 / a$')
+979 if ylabel:
+980 ax1.set_ylabel(ylabel)
+981 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
+982
+983 _handles, labels = ax1.get_legend_handles_labels()
+984 if labels:
+985 ax1.legend()
+986
+987 if title:
+988 plt.title(title)
+989
+990 plt.draw()
+991
+992 if save:
+993 if isinstance(save, str):
+994 fig.savefig(save, bbox_inches='tight')
+995 else:
+996 raise TypeError("'save' has to be a string.")
1003 def spaghetti_plot(self, logscale=True):
-1004 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-1005
-1006 Parameters
-1007 ----------
-1008 logscale : bool
-1009 Determines whether the scale of the y-axis is logarithmic or standard.
-1010 """
-1011 if self.N != 1:
-1012 raise ValueError("Correlator needs to be projected first.")
-1013
-1014 mc_names = list(set([item for sublist in [list(itertools.chain.from_iterable(map(o[0].e_content.get, o[0].mc_names))) for o in self.content if o is not None] for item in sublist]))
-1015 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content, strict=True) if o is not None]
-1016
-1017 for name in mc_names:
-1018 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+ 998 def spaghetti_plot(self, logscale=True):
+ 999 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+1000
+1001 Parameters
+1002 ----------
+1003 logscale : bool
+1004 Determines whether the scale of the y-axis is logarithmic or standard.
+1005 """
+1006 if self.N != 1:
+1007 raise ValueError("Correlator needs to be projected first.")
+1008
+1009 mc_names = list(set([item for sublist in [list(itertools.chain.from_iterable(map(o[0].e_content.get, o[0].mc_names))) for o in self.content if o is not None] for item in sublist]))
+1010 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content, strict=True) if o is not None]
+1011
+1012 for name in mc_names:
+1013 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+1014
+1015 fig = plt.figure()
+1016 ax = fig.add_subplot(111)
+1017 for dat in data:
+1018 ax.plot(x0_vals, dat, ls='-', marker='')
1019
-1020 fig = plt.figure()
-1021 ax = fig.add_subplot(111)
-1022 for dat in data:
-1023 ax.plot(x0_vals, dat, ls='-', marker='')
-1024
-1025 if logscale is True:
-1026 ax.set_yscale('log')
-1027
-1028 ax.set_xlabel(r'$x_0 / a$')
-1029 plt.title(name)
-1030 plt.draw()
+1020 if logscale is True:
+1021 ax.set_yscale('log')
+1022
+1023 ax.set_xlabel(r'$x_0 / a$')
+1024 plt.title(name)
+1025 plt.draw()
1032 def dump(self, filename, datatype="json.gz", **kwargs):
-1033 """Dumps the Corr into a file of chosen type
-1034 Parameters
-1035 ----------
-1036 filename : str
-1037 Name of the file to be saved.
-1038 datatype : str
-1039 Format of the exported file. Supported formats include
-1040 "json.gz" and "pickle"
-1041 path : str
-1042 specifies a custom path for the file (default '.')
-1043 """
-1044 if datatype == "json.gz":
-1045 from .input.json import dump_to_json
-1046 if 'path' in kwargs:
-1047 file_name = kwargs.get('path') + '/' + filename
-1048 else:
-1049 file_name = filename
-1050 dump_to_json(self, file_name)
-1051 elif datatype == "pickle":
-1052 dump_object(self, filename, **kwargs)
-1053 else:
-1054 raise ValueError("Unknown datatype " + str(datatype))
+ 1027 def dump(self, filename, datatype="json.gz", **kwargs):
+1028 """Dumps the Corr into a file of chosen type
+1029 Parameters
+1030 ----------
+1031 filename : str
+1032 Name of the file to be saved.
+1033 datatype : str
+1034 Format of the exported file. Supported formats include
+1035 "json.gz" and "pickle"
+1036 path : str
+1037 specifies a custom path for the file (default '.')
+1038 """
+1039 if datatype == "json.gz":
+1040 from .input.json import dump_to_json
+1041 if 'path' in kwargs:
+1042 file_name = kwargs.get('path') + '/' + filename
+1043 else:
+1044 file_name = filename
+1045 dump_to_json(self, file_name)
+1046 elif datatype == "pickle":
+1047 dump_object(self, filename, **kwargs)
+1048 else:
+1049 raise ValueError("Unknown datatype " + str(datatype))
1056 def print(self, print_range=None):
-1057 print(self.__repr__(print_range))
+ 1051 def print(self, print_range=None):
+1052 print(self.__repr__(print_range))
1275 def sqrt(self):
-1276 return self ** 0.5
+ 1270 def sqrt(self):
+1271 return self ** 0.5
1278 def log(self):
-1279 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
-1280 return Corr(newcontent, prange=self.prange)
+ 1273 def log(self):
+1274 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
+1275 return Corr(newcontent, prange=self.prange)
1282 def exp(self):
-1283 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
-1284 return Corr(newcontent, prange=self.prange)
+ 1277 def exp(self):
+1278 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
+1279 return Corr(newcontent, prange=self.prange)
1299 def sin(self):
-1300 return self._apply_func_to_corr(np.sin)
+ 1294 def sin(self):
+1295 return self._apply_func_to_corr(np.sin)
1302 def cos(self):
-1303 return self._apply_func_to_corr(np.cos)
+ 1297 def cos(self):
+1298 return self._apply_func_to_corr(np.cos)
1305 def tan(self):
-1306 return self._apply_func_to_corr(np.tan)
+ 1300 def tan(self):
+1301 return self._apply_func_to_corr(np.tan)
1308 def sinh(self):
-1309 return self._apply_func_to_corr(np.sinh)
+ 1303 def sinh(self):
+1304 return self._apply_func_to_corr(np.sinh)
1311 def cosh(self):
-1312 return self._apply_func_to_corr(np.cosh)
+ 1306 def cosh(self):
+1307 return self._apply_func_to_corr(np.cosh)
1314 def tanh(self):
-1315 return self._apply_func_to_corr(np.tanh)
+ 1309 def tanh(self):
+1310 return self._apply_func_to_corr(np.tanh)
1317 def arcsin(self):
-1318 return self._apply_func_to_corr(np.arcsin)
+ 1312 def arcsin(self):
+1313 return self._apply_func_to_corr(np.arcsin)
1320 def arccos(self):
-1321 return self._apply_func_to_corr(np.arccos)
+ 1315 def arccos(self):
+1316 return self._apply_func_to_corr(np.arccos)
1323 def arctan(self):
-1324 return self._apply_func_to_corr(np.arctan)
+ 1318 def arctan(self):
+1319 return self._apply_func_to_corr(np.arctan)
1326 def arcsinh(self):
-1327 return self._apply_func_to_corr(np.arcsinh)
+ 1321 def arcsinh(self):
+1322 return self._apply_func_to_corr(np.arcsinh)
1329 def arccosh(self):
-1330 return self._apply_func_to_corr(np.arccosh)
+ 1324 def arccosh(self):
+1325 return self._apply_func_to_corr(np.arccosh)
1332 def arctanh(self):
-1333 return self._apply_func_to_corr(np.arctanh)
+ 1327 def arctanh(self):
+1328 return self._apply_func_to_corr(np.arctanh)
1348 @property
-1349 def real(self):
-1350 def return_real(obs_OR_cobs):
-1351 if isinstance(obs_OR_cobs.flatten()[0], CObs):
-1352 return np.vectorize(lambda x: x.real)(obs_OR_cobs)
-1353 else:
-1354 return obs_OR_cobs
-1355
-1356 return self._apply_func_to_corr(return_real)
+ 1343 @property
+1344 def real(self):
+1345 def return_real(obs_OR_cobs):
+1346 if isinstance(obs_OR_cobs.flatten()[0], CObs):
+1347 return np.vectorize(lambda x: x.real)(obs_OR_cobs)
+1348 else:
+1349 return obs_OR_cobs
+1350
+1351 return self._apply_func_to_corr(return_real)
1358 @property
-1359 def imag(self):
-1360 def return_imag(obs_OR_cobs):
-1361 if isinstance(obs_OR_cobs.flatten()[0], CObs):
-1362 return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
-1363 else:
-1364 return obs_OR_cobs * 0 # So it stays the right type
-1365
-1366 return self._apply_func_to_corr(return_imag)
+ 1353 @property
+1354 def imag(self):
+1355 def return_imag(obs_OR_cobs):
+1356 if isinstance(obs_OR_cobs.flatten()[0], CObs):
+1357 return np.vectorize(lambda x: x.imag)(obs_OR_cobs)
+1358 else:
+1359 return obs_OR_cobs * 0 # So it stays the right type
+1360
+1361 return self._apply_func_to_corr(return_imag)
1368 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1369 r''' Project large correlation matrix to lowest states
-1370
-1371 This method can be used to reduce the size of an (N x N) correlation matrix
-1372 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1373 is still small.
-1374
-1375 Parameters
-1376 ----------
-1377 Ntrunc: int
-1378 Rank of the target matrix.
-1379 tproj: int
-1380 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1381 The default value is 3.
-1382 t0proj: int
-1383 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1384 discouraged for O(a) improved theories, since the correctness of the procedure
-1385 cannot be granted in this case. The default value is 2.
-1386 basematrix : Corr
-1387 Correlation matrix that is used to determine the eigenvectors of the
-1388 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1389 is is not specified.
-1390
-1391 Notes
-1392 -----
-1393 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1394 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
-1395 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1396 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1397 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1398 correlation matrix and to remove some noise that is added by irrelevant operators.
-1399 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1400 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1401 '''
-1402
-1403 if self.N == 1:
-1404 raise ValueError('Method cannot be applied to one-dimensional correlators.')
-1405 if basematrix is None:
-1406 basematrix = self
-1407 if Ntrunc >= basematrix.N:
-1408 raise ValueError(f'Cannot truncate using Ntrunc >= {basematrix.N}')
-1409 if basematrix.N != self.N:
-1410 raise ValueError('basematrix and targetmatrix have to be of the same size.')
-1411
-1412 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
-1413
-1414 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1415 rmat = []
-1416 for t in range(basematrix.T):
-1417 if self.content[t] is None:
-1418 rmat.append(None)
-1419 else:
-1420 for i in range(Ntrunc):
-1421 for j in range(Ntrunc):
-1422 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1423 rmat.append(np.copy(tmpmat))
-1424
-1425 return Corr(rmat)
+ 1363 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1364 r''' Project large correlation matrix to lowest states
+1365
+1366 This method can be used to reduce the size of an (N x N) correlation matrix
+1367 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1368 is still small.
+1369
+1370 Parameters
+1371 ----------
+1372 Ntrunc: int
+1373 Rank of the target matrix.
+1374 tproj: int
+1375 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1376 The default value is 3.
+1377 t0proj: int
+1378 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1379 discouraged for O(a) improved theories, since the correctness of the procedure
+1380 cannot be granted in this case. The default value is 2.
+1381 basematrix : Corr
+1382 Correlation matrix that is used to determine the eigenvectors of the
+1383 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1384 is is not specified.
+1385
+1386 Notes
+1387 -----
+1388 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1389 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
+1390 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1391 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1392 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1393 correlation matrix and to remove some noise that is added by irrelevant operators.
+1394 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1395 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1396 '''
+1397
+1398 if self.N == 1:
+1399 raise ValueError('Method cannot be applied to one-dimensional correlators.')
+1400 if basematrix is None:
+1401 basematrix = self
+1402 if Ntrunc >= basematrix.N:
+1403 raise ValueError(f'Cannot truncate using Ntrunc >= {basematrix.N}')
+1404 if basematrix.N != self.N:
+1405 raise ValueError('basematrix and targetmatrix have to be of the same size.')
+1406
+1407 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
+1408
+1409 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1410 rmat = []
+1411 for t in range(basematrix.T):
+1412 if self.content[t] is None:
+1413 rmat.append(None)
+1414 else:
+1415 for i in range(Ntrunc):
+1416 for j in range(Ntrunc):
+1417 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1418 rmat.append(np.copy(tmpmat))
+1419
+1420 return Corr(rmat)
408 def details(self, ens_content=True):
-409 """Output detailed properties of the Obs.
-410
-411 Parameters
-412 ----------
-413 ens_content : bool
-414 print details about the ensembles and replica if true.
-415 """
-416 if self.tag is not None:
-417 print("Description:", self.tag)
-418 if not hasattr(self, 'e_dvalue'):
-419 print(f'Result\t {self.value:3.8e}')
-420 else:
-421 if self.value == 0.0:
-422 percentage = np.nan
-423 else:
-424 percentage = np.abs(self._dvalue / self.value) * 100
-425 print(f'Result\t {self.value:3.8e} +/- {self._dvalue:3.8e} +/- {self.ddvalue:3.8e} ({percentage:3.3f}%)')
-426 if len(self.e_names) > 1:
-427 print(' Ensemble errors:')
-428 e_content = self.e_content
-429 for e_name in self.mc_names:
-430 gap = _determine_gap(self, e_content, e_name)
-431
-432 if len(self.e_names) > 1:
-433 print('', e_name, f'\t {self.e_dvalue[e_name]:3.6e} +/- {self.e_ddvalue[e_name]:3.6e}')
-434 tau_string = " \N{GREEK SMALL LETTER TAU}_int\t " + _format_uncertainty(self.e_tauint[e_name], self.e_dtauint[e_name])
-435 tau_string += f" in units of {gap} config"
-436 if gap > 1:
-437 tau_string += "s"
-438 if self.tau_exp[e_name] > 0:
-439 tau_string = f"{tau_string: <45}" + f'\t(\N{GREEK SMALL LETTER TAU}_exp={self.tau_exp[e_name]:3.2f}, N_\N{GREEK SMALL LETTER SIGMA}={self.N_sigma[e_name]:g})'
-440 else:
-441 tau_string = f"{tau_string: <45}" + f'\t(S={self.S[e_name]:3.2f})'
-442 print(tau_string)
-443 for e_name in self.cov_names:
-444 print('', e_name, f'\t {self.e_dvalue[e_name]:3.8e}')
-445 if ens_content is True:
-446 if len(self.e_names) == 1:
-447 print(self.N, 'samples in', len(self.e_names), 'ensemble:')
-448 else:
-449 print(self.N, 'samples in', len(self.e_names), 'ensembles:')
-450 my_string_list = []
-451 for key, value in sorted(self.e_content.items()):
-452 if key not in self.covobs:
-453 my_string = ' ' + "\u00B7 Ensemble '" + key + "' "
-454 if len(value) == 1:
-455 my_string += f': {self.shape[value[0]]} configurations'
-456 if isinstance(self.idl[value[0]], range):
-457 my_string += f' (from {self.idl[value[0]].start} to {self.idl[value[0]][-1]}' + int(self.idl[value[0]].step != 1) * f' in steps of {self.idl[value[0]].step}' + ')'
-458 else:
-459 my_string += f' (irregular range from {self.idl[value[0]][0]} to {self.idl[value[0]][-1]})'
-460 else:
-461 sublist = []
-462 for v in value:
-463 my_substring = ' ' + "\u00B7 Replicum '" + v[len(key) + 1:] + "' "
-464 my_substring += f': {self.shape[v]} configurations'
-465 if isinstance(self.idl[v], range):
-466 my_substring += f' (from {self.idl[v].start} to {self.idl[v][-1]}' + int(self.idl[v].step != 1) * f' in steps of {self.idl[v].step}' + ')'
-467 else:
-468 my_substring += f' (irregular range from {self.idl[v][0]} to {self.idl[v][-1]})'
-469 sublist.append(my_substring)
-470
-471 my_string += '\n' + '\n'.join(sublist)
-472 else:
-473 my_string = ' ' + "\u00B7 Covobs '" + key + "' "
-474 my_string_list.append(my_string)
-475 print('\n'.join(my_string_list))
+ 407 def details(self, ens_content=True):
+408 """Output detailed properties of the Obs.
+409
+410 Parameters
+411 ----------
+412 ens_content : bool
+413 print details about the ensembles and replica if true.
+414 """
+415 if self.tag is not None:
+416 print("Description:", self.tag)
+417 if not hasattr(self, 'e_dvalue'):
+418 print(f'Result\t {self.value:3.8e}')
+419 else:
+420 if self.value == 0.0:
+421 percentage = np.nan
+422 else:
+423 percentage = np.abs(self._dvalue / self.value) * 100
+424 print(f'Result\t {self.value:3.8e} +/- {self._dvalue:3.8e} +/- {self.ddvalue:3.8e} ({percentage:3.3f}%)')
+425 if len(self.e_names) > 1:
+426 print(' Ensemble errors:')
+427 e_content = self.e_content
+428 for e_name in self.mc_names:
+429 gap = _determine_gap(self, e_content, e_name)
+430
+431 if len(self.e_names) > 1:
+432 print('', e_name, f'\t {self.e_dvalue[e_name]:3.6e} +/- {self.e_ddvalue[e_name]:3.6e}')
+433 tau_string = " \N{GREEK SMALL LETTER TAU}_int\t " + _format_uncertainty(self.e_tauint[e_name], self.e_dtauint[e_name])
+434 tau_string += f" in units of {gap} config"
+435 if gap > 1:
+436 tau_string += "s"
+437 if self.tau_exp[e_name] > 0:
+438 tau_string = f"{tau_string: <45}" + f'\t(\N{GREEK SMALL LETTER TAU}_exp={self.tau_exp[e_name]:3.2f}, N_\N{GREEK SMALL LETTER SIGMA}={self.N_sigma[e_name]:g})'
+439 else:
+440 tau_string = f"{tau_string: <45}" + f'\t(S={self.S[e_name]:3.2f})'
+441 print(tau_string)
+442 for e_name in self.cov_names:
+443 print('', e_name, f'\t {self.e_dvalue[e_name]:3.8e}')
+444 if ens_content is True:
+445 if len(self.e_names) == 1:
+446 print(self.N, 'samples in', len(self.e_names), 'ensemble:')
+447 else:
+448 print(self.N, 'samples in', len(self.e_names), 'ensembles:')
+449 my_string_list = []
+450 for key, value in sorted(self.e_content.items()):
+451 if key not in self.covobs:
+452 my_string = ' ' + "\u00B7 Ensemble '" + key + "' "
+453 if len(value) == 1:
+454 my_string += f': {self.shape[value[0]]} configurations'
+455 if isinstance(self.idl[value[0]], range):
+456 my_string += f' (from {self.idl[value[0]].start} to {self.idl[value[0]][-1]}' + int(self.idl[value[0]].step != 1) * f' in steps of {self.idl[value[0]].step}' + ')'
+457 else:
+458 my_string += f' (irregular range from {self.idl[value[0]][0]} to {self.idl[value[0]][-1]})'
+459 else:
+460 sublist = []
+461 for v in value:
+462 my_substring = ' ' + "\u00B7 Replicum '" + v[len(key) + 1:] + "' "
+463 my_substring += f': {self.shape[v]} configurations'
+464 if isinstance(self.idl[v], range):
+465 my_substring += f' (from {self.idl[v].start} to {self.idl[v][-1]}' + int(self.idl[v].step != 1) * f' in steps of {self.idl[v].step}' + ')'
+466 else:
+467 my_substring += f' (irregular range from {self.idl[v][0]} to {self.idl[v][-1]})'
+468 sublist.append(my_substring)
+469
+470 my_string += '\n' + '\n'.join(sublist)
+471 else:
+472 my_string = ' ' + "\u00B7 Covobs '" + key + "' "
+473 my_string_list.append(my_string)
+474 print('\n'.join(my_string_list))
477 def reweight(self, weight):
-478 """Reweight the obs with given rewighting factors.
-479
-480 Parameters
-481 ----------
-482 weight : Obs
-483 Reweighting factor. An Observable that has to be defined on a superset of the
-484 configurations in obs[i].idl for all i.
-485 all_configs : bool
-486 if True, the reweighted observables are normalized by the average of
-487 the reweighting factor on all configurations in weight.idl and not
-488 on the configurations in obs[i].idl. Default False.
-489 """
-490 return reweight(weight, [self])[0]
+ 476 def reweight(self, weight):
+477 """Reweight the obs with given rewighting factors.
+478
+479 Parameters
+480 ----------
+481 weight : Obs
+482 Reweighting factor. An Observable that has to be defined on a superset of the
+483 configurations in obs[i].idl for all i.
+484 all_configs : bool
+485 if True, the reweighted observables are normalized by the average of
+486 the reweighting factor on all configurations in weight.idl and not
+487 on the configurations in obs[i].idl. Default False.
+488 """
+489 return reweight(weight, [self])[0]
492 def is_zero_within_error(self, sigma=1):
-493 """Checks whether the observable is zero within 'sigma' standard errors.
-494
-495 Parameters
-496 ----------
-497 sigma : int
-498 Number of standard errors used for the check.
-499
-500 Works only properly when the gamma method was run.
-501 """
-502 return self.is_zero() or np.abs(self.value) <= sigma * self._dvalue
+ 491 def is_zero_within_error(self, sigma=1):
+492 """Checks whether the observable is zero within 'sigma' standard errors.
+493
+494 Parameters
+495 ----------
+496 sigma : int
+497 Number of standard errors used for the check.
+498
+499 Works only properly when the gamma method was run.
+500 """
+501 return self.is_zero() or np.abs(self.value) <= sigma * self._dvalue
504 def is_zero(self, atol=1e-10):
-505 """Checks whether the observable is zero within a given tolerance.
-506
-507 Parameters
-508 ----------
-509 atol : float
-510 Absolute tolerance (for details see numpy documentation).
-511 """
-512 return np.isclose(0.0, self.value, 1e-14, atol) and all(np.allclose(0.0, delta, 1e-14, atol) for delta in self.deltas.values()) and all(np.allclose(0.0, delta.errsq(), 1e-14, atol) for delta in self.covobs.values())
+ 503 def is_zero(self, atol=1e-10):
+504 """Checks whether the observable is zero within a given tolerance.
+505
+506 Parameters
+507 ----------
+508 atol : float
+509 Absolute tolerance (for details see numpy documentation).
+510 """
+511 return np.isclose(0.0, self.value, 1e-14, atol) and all(np.allclose(0.0, delta, 1e-14, atol) for delta in self.deltas.values()) and all(np.allclose(0.0, delta.errsq(), 1e-14, atol) for delta in self.covobs.values())
514 def plot_tauint(self, save=None):
-515 """Plot integrated autocorrelation time for each ensemble.
-516
-517 Parameters
-518 ----------
-519 save : str
-520 saves the figure to a file named 'save' if.
-521 """
-522 if not hasattr(self, 'e_dvalue'):
-523 raise Exception('Run the gamma method first.')
-524
-525 for e, e_name in enumerate(self.mc_names):
-526 fig = plt.figure()
-527 plt.xlabel(r'$W$')
-528 plt.ylabel(r'$\tau_\mathrm{int}$')
-529 length = len(self.e_n_tauint[e_name])
-530 if self.tau_exp[e_name] > 0:
-531 base = self.e_n_tauint[e_name][self.e_windowsize[e_name]]
-532 x_help = np.arange(2 * self.tau_exp[e_name])
-533 y_help = (x_help + 1) * np.abs(self.e_rho[e_name][self.e_windowsize[e_name] + 1]) * (1 - x_help / (2 * (2 * self.tau_exp[e_name] - 1))) + base
-534 x_arr = np.arange(self.e_windowsize[e_name] + 1, self.e_windowsize[e_name] + 1 + 2 * self.tau_exp[e_name])
-535 plt.plot(x_arr, y_help, 'C' + str(e), linewidth=1, ls='--', marker=',')
-536 plt.errorbar([self.e_windowsize[e_name] + 2 * self.tau_exp[e_name]], [self.e_tauint[e_name]],
-537 yerr=[self.e_dtauint[e_name]], fmt='C' + str(e), linewidth=1, capsize=2, marker='o', mfc=plt.rcParams['axes.facecolor'])
-538 xmax = self.e_windowsize[e_name] + 2 * self.tau_exp[e_name] + 1.5
-539 label = e_name + r', $\tau_\mathrm{exp}$=' + str(np.around(self.tau_exp[e_name], decimals=2))
-540 else:
-541 label = e_name + ', S=' + str(np.around(self.S[e_name], decimals=2))
-542 xmax = max(10.5, 2 * self.e_windowsize[e_name] - 0.5)
-543
-544 plt.errorbar(np.arange(length)[:int(xmax) + 1], self.e_n_tauint[e_name][:int(xmax) + 1], yerr=self.e_n_dtauint[e_name][:int(xmax) + 1], linewidth=1, capsize=2, label=label)
-545 plt.axvline(x=self.e_windowsize[e_name], color='C' + str(e), alpha=0.5, marker=',', ls='--')
-546 plt.legend()
-547 plt.xlim(-0.5, xmax)
-548 ylim = plt.ylim()
-549 plt.ylim(bottom=0.0, top=max(1.0, ylim[1]))
-550 plt.draw()
-551 if save:
-552 fig.savefig(save + "_" + str(e))
+ 513 def plot_tauint(self, save=None):
+514 """Plot integrated autocorrelation time for each ensemble.
+515
+516 Parameters
+517 ----------
+518 save : str
+519 saves the figure to a file named 'save' if.
+520 """
+521 if not hasattr(self, 'e_dvalue'):
+522 raise Exception('Run the gamma method first.')
+523
+524 for e, e_name in enumerate(self.mc_names):
+525 fig = plt.figure()
+526 plt.xlabel(r'$W$')
+527 plt.ylabel(r'$\tau_\mathrm{int}$')
+528 length = len(self.e_n_tauint[e_name])
+529 if self.tau_exp[e_name] > 0:
+530 base = self.e_n_tauint[e_name][self.e_windowsize[e_name]]
+531 x_help = np.arange(2 * self.tau_exp[e_name])
+532 y_help = (x_help + 1) * np.abs(self.e_rho[e_name][self.e_windowsize[e_name] + 1]) * (1 - x_help / (2 * (2 * self.tau_exp[e_name] - 1))) + base
+533 x_arr = np.arange(self.e_windowsize[e_name] + 1, self.e_windowsize[e_name] + 1 + 2 * self.tau_exp[e_name])
+534 plt.plot(x_arr, y_help, 'C' + str(e), linewidth=1, ls='--', marker=',')
+535 plt.errorbar([self.e_windowsize[e_name] + 2 * self.tau_exp[e_name]], [self.e_tauint[e_name]],
+536 yerr=[self.e_dtauint[e_name]], fmt='C' + str(e), linewidth=1, capsize=2, marker='o', mfc=plt.rcParams['axes.facecolor'])
+537 xmax = self.e_windowsize[e_name] + 2 * self.tau_exp[e_name] + 1.5
+538 label = e_name + r', $\tau_\mathrm{exp}$=' + str(np.around(self.tau_exp[e_name], decimals=2))
+539 else:
+540 label = e_name + ', S=' + str(np.around(self.S[e_name], decimals=2))
+541 xmax = max(10.5, 2 * self.e_windowsize[e_name] - 0.5)
+542
+543 plt.errorbar(np.arange(length)[:int(xmax) + 1], self.e_n_tauint[e_name][:int(xmax) + 1], yerr=self.e_n_dtauint[e_name][:int(xmax) + 1], linewidth=1, capsize=2, label=label)
+544 plt.axvline(x=self.e_windowsize[e_name], color='C' + str(e), alpha=0.5, marker=',', ls='--')
+545 plt.legend()
+546 plt.xlim(-0.5, xmax)
+547 ylim = plt.ylim()
+548 plt.ylim(bottom=0.0, top=max(1.0, ylim[1]))
+549 plt.draw()
+550 if save:
+551 fig.savefig(save + "_" + str(e))
554 def plot_rho(self, save=None):
-555 """Plot normalized autocorrelation function time for each ensemble.
-556
-557 Parameters
-558 ----------
-559 save : str
-560 saves the figure to a file named 'save' if.
-561 """
-562 if not hasattr(self, 'e_dvalue'):
-563 raise Exception('Run the gamma method first.')
-564 for e, e_name in enumerate(self.mc_names):
-565 fig = plt.figure()
-566 plt.xlabel('W')
-567 plt.ylabel('rho')
-568 length = len(self.e_drho[e_name])
-569 plt.errorbar(np.arange(length), self.e_rho[e_name][:length], yerr=self.e_drho[e_name][:], linewidth=1, capsize=2)
-570 plt.axvline(x=self.e_windowsize[e_name], color='r', alpha=0.25, ls='--', marker=',')
-571 if self.tau_exp[e_name] > 0:
-572 plt.plot([self.e_windowsize[e_name] + 1, self.e_windowsize[e_name] + 1 + 2 * self.tau_exp[e_name]],
-573 [self.e_rho[e_name][self.e_windowsize[e_name] + 1], 0], 'k-', lw=1)
-574 xmax = self.e_windowsize[e_name] + 2 * self.tau_exp[e_name] + 1.5
-575 plt.title('Rho ' + e_name + r', tau\_exp=' + str(np.around(self.tau_exp[e_name], decimals=2)))
-576 else:
-577 xmax = max(10.5, 2 * self.e_windowsize[e_name] - 0.5)
-578 plt.title('Rho ' + e_name + ', S=' + str(np.around(self.S[e_name], decimals=2)))
-579 plt.plot([-0.5, xmax], [0, 0], 'k--', lw=1)
-580 plt.xlim(-0.5, xmax)
-581 plt.draw()
-582 if save:
-583 fig.savefig(save + "_" + str(e))
+ 553 def plot_rho(self, save=None):
+554 """Plot normalized autocorrelation function time for each ensemble.
+555
+556 Parameters
+557 ----------
+558 save : str
+559 saves the figure to a file named 'save' if.
+560 """
+561 if not hasattr(self, 'e_dvalue'):
+562 raise Exception('Run the gamma method first.')
+563 for e, e_name in enumerate(self.mc_names):
+564 fig = plt.figure()
+565 plt.xlabel('W')
+566 plt.ylabel('rho')
+567 length = len(self.e_drho[e_name])
+568 plt.errorbar(np.arange(length), self.e_rho[e_name][:length], yerr=self.e_drho[e_name][:], linewidth=1, capsize=2)
+569 plt.axvline(x=self.e_windowsize[e_name], color='r', alpha=0.25, ls='--', marker=',')
+570 if self.tau_exp[e_name] > 0:
+571 plt.plot([self.e_windowsize[e_name] + 1, self.e_windowsize[e_name] + 1 + 2 * self.tau_exp[e_name]],
+572 [self.e_rho[e_name][self.e_windowsize[e_name] + 1], 0], 'k-', lw=1)
+573 xmax = self.e_windowsize[e_name] + 2 * self.tau_exp[e_name] + 1.5
+574 plt.title('Rho ' + e_name + r', tau\_exp=' + str(np.around(self.tau_exp[e_name], decimals=2)))
+575 else:
+576 xmax = max(10.5, 2 * self.e_windowsize[e_name] - 0.5)
+577 plt.title('Rho ' + e_name + ', S=' + str(np.around(self.S[e_name], decimals=2)))
+578 plt.plot([-0.5, xmax], [0, 0], 'k--', lw=1)
+579 plt.xlim(-0.5, xmax)
+580 plt.draw()
+581 if save:
+582 fig.savefig(save + "_" + str(e))
585 def plot_rep_dist(self):
-586 """Plot replica distribution for each ensemble with more than one replicum."""
-587 if not hasattr(self, 'e_dvalue'):
-588 raise Exception('Run the gamma method first.')
-589 for _e, e_name in enumerate(self.mc_names):
-590 if len(self.e_content[e_name]) == 1:
-591 print('No replica distribution for a single replicum (', e_name, ')')
-592 continue
-593 r_length = []
-594 sub_r_mean = 0
-595 for r_name in self.e_content[e_name]:
-596 r_length.append(len(self.deltas[r_name]))
-597 sub_r_mean += self.shape[r_name] * self.r_values[r_name]
-598 e_N = np.sum(r_length)
-599 sub_r_mean /= e_N
-600 arr = np.zeros(len(self.e_content[e_name]))
-601 for r, r_name in enumerate(self.e_content[e_name]):
-602 arr[r] = (self.r_values[r_name] - sub_r_mean) / (self.e_dvalue[e_name] * np.sqrt(e_N / self.shape[r_name] - 1))
-603 plt.hist(arr, rwidth=0.8, bins=len(self.e_content[e_name]))
-604 plt.title('Replica distribution' + e_name + ' (mean=0, var=1)')
-605 plt.draw()
+ 584 def plot_rep_dist(self):
+585 """Plot replica distribution for each ensemble with more than one replicum."""
+586 if not hasattr(self, 'e_dvalue'):
+587 raise Exception('Run the gamma method first.')
+588 for _e, e_name in enumerate(self.mc_names):
+589 if len(self.e_content[e_name]) == 1:
+590 print('No replica distribution for a single replicum (', e_name, ')')
+591 continue
+592 r_length = []
+593 sub_r_mean = 0
+594 for r_name in self.e_content[e_name]:
+595 r_length.append(len(self.deltas[r_name]))
+596 sub_r_mean += self.shape[r_name] * self.r_values[r_name]
+597 e_N = np.sum(r_length)
+598 sub_r_mean /= e_N
+599 arr = np.zeros(len(self.e_content[e_name]))
+600 for r, r_name in enumerate(self.e_content[e_name]):
+601 arr[r] = (self.r_values[r_name] - sub_r_mean) / (self.e_dvalue[e_name] * np.sqrt(e_N / self.shape[r_name] - 1))
+602 plt.hist(arr, rwidth=0.8, bins=len(self.e_content[e_name]))
+603 plt.title('Replica distribution' + e_name + ' (mean=0, var=1)')
+604 plt.draw()
607 def plot_history(self, expand=True):
-608 """Plot derived Monte Carlo history for each ensemble
-609
-610 Parameters
-611 ----------
-612 expand : bool
-613 show expanded history for irregular Monte Carlo chains (default: True).
-614 """
-615 for _e, e_name in enumerate(self.mc_names):
-616 plt.figure()
-617 r_length = []
-618 tmp = []
-619 tmp_expanded = []
-620 for _r, r_name in enumerate(self.e_content[e_name]):
-621 tmp.append(self.deltas[r_name] + self.r_values[r_name])
-622 if expand:
-623 tmp_expanded.append(_expand_deltas(self.deltas[r_name], list(self.idl[r_name]), self.shape[r_name], 1) + self.r_values[r_name])
-624 r_length.append(len(tmp_expanded[-1]))
-625 else:
-626 r_length.append(len(tmp[-1]))
-627 e_N = np.sum(r_length)
-628 x = np.arange(e_N)
-629 y_test = np.concatenate(tmp, axis=0)
-630 if expand:
-631 y = np.concatenate(tmp_expanded, axis=0)
-632 else:
-633 y = y_test
-634 plt.errorbar(x, y, fmt='.', markersize=3)
-635 plt.xlim(-0.5, e_N - 0.5)
-636 plt.title(e_name + f'\nskew: {skew(y_test):.3f} (p={skewtest(y_test).pvalue:.3f}), kurtosis: {kurtosis(y_test):.3f} (p={kurtosistest(y_test).pvalue:.3f})')
-637 plt.draw()
+ 606 def plot_history(self, expand=True):
+607 """Plot derived Monte Carlo history for each ensemble
+608
+609 Parameters
+610 ----------
+611 expand : bool
+612 show expanded history for irregular Monte Carlo chains (default: True).
+613 """
+614 for _e, e_name in enumerate(self.mc_names):
+615 plt.figure()
+616 r_length = []
+617 tmp = []
+618 tmp_expanded = []
+619 for _r, r_name in enumerate(self.e_content[e_name]):
+620 tmp.append(self.deltas[r_name] + self.r_values[r_name])
+621 if expand:
+622 tmp_expanded.append(_expand_deltas(self.deltas[r_name], list(self.idl[r_name]), self.shape[r_name], 1) + self.r_values[r_name])
+623 r_length.append(len(tmp_expanded[-1]))
+624 else:
+625 r_length.append(len(tmp[-1]))
+626 e_N = np.sum(r_length)
+627 x = np.arange(e_N)
+628 y_test = np.concatenate(tmp, axis=0)
+629 if expand:
+630 y = np.concatenate(tmp_expanded, axis=0)
+631 else:
+632 y = y_test
+633 plt.errorbar(x, y, fmt='.', markersize=3)
+634 plt.xlim(-0.5, e_N - 0.5)
+635 plt.title(e_name + f'\nskew: {skew(y_test):.3f} (p={skewtest(y_test).pvalue:.3f}), kurtosis: {kurtosis(y_test):.3f} (p={kurtosistest(y_test).pvalue:.3f})')
+636 plt.draw()
639 def plot_piechart(self, save=None):
-640 """Plot piechart which shows the fractional contribution of each
-641 ensemble to the error and returns a dictionary containing the fractions.
-642
-643 Parameters
-644 ----------
-645 save : str
-646 saves the figure to a file named 'save' if.
-647 """
-648 if not hasattr(self, 'e_dvalue'):
-649 raise Exception('Run the gamma method first.')
-650 if np.isclose(0.0, self._dvalue, atol=1e-15):
-651 raise ValueError('Error is 0.0')
-652 labels = self.e_names
-653 sizes = [self.e_dvalue[name] ** 2 for name in labels] / self._dvalue ** 2
-654 fig1, ax1 = plt.subplots()
-655 ax1.pie(sizes, labels=labels, startangle=90, normalize=True)
-656 ax1.axis('equal')
-657 plt.draw()
-658 if save:
-659 fig1.savefig(save)
-660
-661 return dict(zip(labels, sizes, strict=True))
+ 638 def plot_piechart(self, save=None):
+639 """Plot piechart which shows the fractional contribution of each
+640 ensemble to the error and returns a dictionary containing the fractions.
+641
+642 Parameters
+643 ----------
+644 save : str
+645 saves the figure to a file named 'save' if.
+646 """
+647 if not hasattr(self, 'e_dvalue'):
+648 raise Exception('Run the gamma method first.')
+649 if np.isclose(0.0, self._dvalue, atol=1e-15):
+650 raise ValueError('Error is 0.0')
+651 labels = self.e_names
+652 sizes = [self.e_dvalue[name] ** 2 for name in labels] / self._dvalue ** 2
+653 fig1, ax1 = plt.subplots()
+654 ax1.pie(sizes, labels=labels, startangle=90, normalize=True)
+655 ax1.axis('equal')
+656 plt.draw()
+657 if save:
+658 fig1.savefig(save)
+659
+660 return dict(zip(labels, sizes, strict=True))
663 def dump(self, filename, datatype="json.gz", description="", **kwargs):
-664 """Dump the Obs to a file 'name' of chosen format.
-665
-666 Parameters
-667 ----------
-668 filename : str
-669 name of the file to be saved.
-670 datatype : str
-671 Format of the exported file. Supported formats include
-672 "json.gz" and "pickle"
-673 description : str
-674 Description for output file, only relevant for json.gz format.
-675 path : str
-676 specifies a custom path for the file (default '.')
-677 """
-678 if 'path' in kwargs:
-679 file_name = kwargs.get('path') + '/' + filename
-680 else:
-681 file_name = filename
-682
-683 if datatype == "json.gz":
-684 from .input.json import dump_to_json
-685 dump_to_json([self], file_name, description=description)
-686 elif datatype == "pickle":
-687 with open(file_name + '.p', 'wb') as fb:
-688 pickle.dump(self, fb)
-689 else:
-690 raise TypeError("Unknown datatype " + str(datatype))
+ 662 def dump(self, filename, datatype="json.gz", description="", **kwargs):
+663 """Dump the Obs to a file 'name' of chosen format.
+664
+665 Parameters
+666 ----------
+667 filename : str
+668 name of the file to be saved.
+669 datatype : str
+670 Format of the exported file. Supported formats include
+671 "json.gz" and "pickle"
+672 description : str
+673 Description for output file, only relevant for json.gz format.
+674 path : str
+675 specifies a custom path for the file (default '.')
+676 """
+677 if 'path' in kwargs:
+678 file_name = kwargs.get('path') + '/' + filename
+679 else:
+680 file_name = filename
+681
+682 if datatype == "json.gz":
+683 from .input.json import dump_to_json
+684 dump_to_json([self], file_name, description=description)
+685 elif datatype == "pickle":
+686 with open(file_name + '.p', 'wb') as fb:
+687 pickle.dump(self, fb)
+688 else:
+689 raise TypeError("Unknown datatype " + str(datatype))
692 def export_jackknife(self):
-693 """Export jackknife samples from the Obs
-694
-695 Returns
-696 -------
-697 numpy.ndarray
-698 Returns a numpy array of length N + 1 where N is the number of samples
-699 for the given ensemble and replicum. The zeroth entry of the array contains
-700 the mean value of the Obs, entries 1 to N contain the N jackknife samples
-701 derived from the Obs. The current implementation only works for observables
-702 defined on exactly one ensemble and replicum. The derived jackknife samples
-703 should agree with samples from a full jackknife analysis up to O(1/N).
-704 """
-705
-706 if len(self.names) != 1:
-707 raise ValueError("'export_jackknife' is only implemented for Obs defined on one ensemble and replicum.")
-708
-709 name = self.names[0]
-710 full_data = self.deltas[name] + self.r_values[name]
-711 n = full_data.size
-712 mean = self.value
-713 tmp_jacks = np.zeros(n + 1)
-714 tmp_jacks[0] = mean
-715 tmp_jacks[1:] = (n * mean - full_data) / (n - 1)
-716 return tmp_jacks
+ 691 def export_jackknife(self):
+692 """Export jackknife samples from the Obs
+693
+694 Returns
+695 -------
+696 numpy.ndarray
+697 Returns a numpy array of length N + 1 where N is the number of samples
+698 for the given ensemble and replicum. The zeroth entry of the array contains
+699 the mean value of the Obs, entries 1 to N contain the N jackknife samples
+700 derived from the Obs. The current implementation only works for observables
+701 defined on exactly one ensemble and replicum. The derived jackknife samples
+702 should agree with samples from a full jackknife analysis up to O(1/N).
+703 """
+704
+705 if len(self.names) != 1:
+706 raise ValueError("'export_jackknife' is only implemented for Obs defined on one ensemble and replicum.")
+707
+708 name = self.names[0]
+709 full_data = self.deltas[name] + self.r_values[name]
+710 n = full_data.size
+711 mean = self.value
+712 tmp_jacks = np.zeros(n + 1)
+713 tmp_jacks[0] = mean
+714 tmp_jacks[1:] = (n * mean - full_data) / (n - 1)
+715 return tmp_jacks
718 def export_bootstrap(self, samples=500, random_numbers=None, save_rng=None):
-719 """Export bootstrap samples from the Obs
-720
-721 Parameters
-722 ----------
-723 samples : int
-724 Number of bootstrap samples to generate.
-725 random_numbers : np.ndarray
-726 Array of shape (samples, length) containing the random numbers to generate the bootstrap samples.
-727 If not provided the bootstrap samples are generated bashed on the md5 hash of the enesmble name.
-728 save_rng : str
-729 Save the random numbers to a file if a path is specified.
-730
-731 Returns
-732 -------
-733 numpy.ndarray
-734 Returns a numpy array of length N + 1 where N is the number of samples
-735 for the given ensemble and replicum. The zeroth entry of the array contains
-736 the mean value of the Obs, entries 1 to N contain the N import_bootstrap samples
-737 derived from the Obs. The current implementation only works for observables
-738 defined on exactly one ensemble and replicum. The derived bootstrap samples
-739 should agree with samples from a full bootstrap analysis up to O(1/N).
-740 """
-741 if len(self.names) != 1:
-742 raise ValueError("'export_boostrap' is only implemented for Obs defined on one ensemble and replicum.")
-743
-744 name = self.names[0]
-745 length = self.N
-746
-747 if random_numbers is None:
-748 seed = int(hashlib.md5(name.encode()).hexdigest(), 16) & 0xFFFFFFFF
-749 rng = np.random.default_rng(seed)
-750 random_numbers = rng.integers(0, length, size=(samples, length))
-751
-752 if save_rng is not None:
-753 np.savetxt(save_rng, random_numbers, fmt='%i')
-754
-755 proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
-756 ret = np.zeros(samples + 1)
-757 ret[0] = self.value
-758 ret[1:] = proj @ (self.deltas[name] + self.r_values[name])
-759 return ret
+ 717 def export_bootstrap(self, samples=500, random_numbers=None, save_rng=None):
+718 """Export bootstrap samples from the Obs
+719
+720 Parameters
+721 ----------
+722 samples : int
+723 Number of bootstrap samples to generate.
+724 random_numbers : np.ndarray
+725 Array of shape (samples, length) containing the random numbers to generate the bootstrap samples.
+726 If not provided the bootstrap samples are generated bashed on the md5 hash of the enesmble name.
+727 save_rng : str
+728 Save the random numbers to a file if a path is specified.
+729
+730 Returns
+731 -------
+732 numpy.ndarray
+733 Returns a numpy array of length N + 1 where N is the number of samples
+734 for the given ensemble and replicum. The zeroth entry of the array contains
+735 the mean value of the Obs, entries 1 to N contain the N import_bootstrap samples
+736 derived from the Obs. The current implementation only works for observables
+737 defined on exactly one ensemble and replicum. The derived bootstrap samples
+738 should agree with samples from a full bootstrap analysis up to O(1/N).
+739 """
+740 if len(self.names) != 1:
+741 raise ValueError("'export_boostrap' is only implemented for Obs defined on one ensemble and replicum.")
+742
+743 name = self.names[0]
+744 length = self.N
+745
+746 if random_numbers is None:
+747 seed = int(hashlib.md5(name.encode()).hexdigest(), 16) & 0xFFFFFFFF
+748 rng = np.random.default_rng(seed)
+749 random_numbers = rng.integers(0, length, size=(samples, length))
+750
+751 if save_rng is not None:
+752 np.savetxt(save_rng, random_numbers, fmt='%i')
+753
+754 proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
+755 ret = np.zeros(samples + 1)
+756 ret[0] = self.value
+757 ret[1:] = proj @ (self.deltas[name] + self.r_values[name])
+758 return ret
898 def sqrt(self):
-899 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
+ 897 def sqrt(self):
+898 return derived_observable(lambda x, **kwargs: np.sqrt(x[0]), [self], man_grad=[1 / 2 / np.sqrt(self.value)])
901 def log(self):
-902 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
+ 900 def log(self):
+901 return derived_observable(lambda x, **kwargs: np.log(x[0]), [self], man_grad=[1 / self.value])
904 def exp(self):
-905 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
+ 903 def exp(self):
+904 return derived_observable(lambda x, **kwargs: np.exp(x[0]), [self], man_grad=[np.exp(self.value)])
907 def sin(self):
-908 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
+ 906 def sin(self):
+907 return derived_observable(lambda x, **kwargs: np.sin(x[0]), [self], man_grad=[np.cos(self.value)])
910 def cos(self):
-911 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
+ 909 def cos(self):
+910 return derived_observable(lambda x, **kwargs: np.cos(x[0]), [self], man_grad=[-np.sin(self.value)])
913 def tan(self):
-914 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
+ 912 def tan(self):
+913 return derived_observable(lambda x, **kwargs: np.tan(x[0]), [self], man_grad=[1 / np.cos(self.value) ** 2])
916 def arcsin(self):
-917 return derived_observable(lambda x: anp.arcsin(x[0]), [self])
+ 915 def arcsin(self):
+916 return derived_observable(lambda x: anp.arcsin(x[0]), [self])
919 def arccos(self):
-920 return derived_observable(lambda x: anp.arccos(x[0]), [self])
+ 918 def arccos(self):
+919 return derived_observable(lambda x: anp.arccos(x[0]), [self])
922 def arctan(self):
-923 return derived_observable(lambda x: anp.arctan(x[0]), [self])
+ 921 def arctan(self):
+922 return derived_observable(lambda x: anp.arctan(x[0]), [self])
925 def sinh(self):
-926 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
+ 924 def sinh(self):
+925 return derived_observable(lambda x, **kwargs: np.sinh(x[0]), [self], man_grad=[np.cosh(self.value)])
928 def cosh(self):
-929 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
+ 927 def cosh(self):
+928 return derived_observable(lambda x, **kwargs: np.cosh(x[0]), [self], man_grad=[np.sinh(self.value)])
931 def tanh(self):
-932 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
+ 930 def tanh(self):
+931 return derived_observable(lambda x, **kwargs: np.tanh(x[0]), [self], man_grad=[1 / np.cosh(self.value) ** 2])
934 def arcsinh(self):
-935 return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
+ 933 def arcsinh(self):
+934 return derived_observable(lambda x: anp.arcsinh(x[0]), [self])
937 def arccosh(self):
-938 return derived_observable(lambda x: anp.arccosh(x[0]), [self])
+ 936 def arccosh(self):
+937 return derived_observable(lambda x: anp.arccosh(x[0]), [self])
940 def arctanh(self):
-941 return derived_observable(lambda x: anp.arctanh(x[0]), [self])
+ 939 def arctanh(self):
+940 return derived_observable(lambda x: anp.arctanh(x[0]), [self])
944class CObs:
- 945 """Class for a complex valued observable."""
- 946 __slots__ = ['_imag', '_real', 'tag']
- 947
- 948 def __init__(self, real, imag=0.0):
- 949 self._real = real
- 950 self._imag = imag
- 951 self.tag = None
- 952
- 953 @property
- 954 def real(self):
- 955 return self._real
- 956
- 957 @property
- 958 def imag(self):
- 959 return self._imag
- 960
- 961 def gamma_method(self, **kwargs):
- 962 """Executes the gamma_method for the real and the imaginary part."""
- 963 if isinstance(self.real, Obs):
- 964 self.real.gamma_method(**kwargs)
- 965 if isinstance(self.imag, Obs):
- 966 self.imag.gamma_method(**kwargs)
- 967
- 968 def is_zero(self):
- 969 """Checks whether both real and imaginary part are zero within machine precision."""
- 970 return self.real == 0.0 and self.imag == 0.0
- 971
- 972 def conjugate(self):
- 973 return CObs(self.real, -self.imag)
- 974
- 975 def __add__(self, other):
- 976 if isinstance(other, np.ndarray):
- 977 return other + self
- 978 elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 979 return CObs(self.real + other.real,
- 980 self.imag + other.imag)
- 981 else:
- 982 return CObs(self.real + other, self.imag)
- 983
- 984 def __radd__(self, y):
- 985 return self + y
- 986
- 987 def __sub__(self, other):
- 988 if isinstance(other, np.ndarray):
- 989 return -1 * (other - self)
- 990 elif hasattr(other, 'real') and hasattr(other, 'imag'):
- 991 return CObs(self.real - other.real, self.imag - other.imag)
- 992 else:
- 993 return CObs(self.real - other, self.imag)
- 994
- 995 def __rsub__(self, other):
- 996 return -1 * (self - other)
- 997
- 998 def __mul__(self, other):
- 999 if isinstance(other, np.ndarray):
-1000 return other * self
-1001 elif hasattr(other, 'real') and hasattr(other, 'imag'):
-1002 if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
-1003 return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
-1004 [self.real, other.real, self.imag, other.imag],
-1005 man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
-1006 derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
-1007 [self.real, other.real, self.imag, other.imag],
-1008 man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
-1009 elif getattr(other, 'imag', 0) != 0:
-1010 return CObs(self.real * other.real - self.imag * other.imag,
-1011 self.imag * other.real + self.real * other.imag)
-1012 else:
-1013 return CObs(self.real * other.real, self.imag * other.real)
-1014 else:
-1015 return CObs(self.real * other, self.imag * other)
-1016
-1017 def __rmul__(self, other):
-1018 return self * other
-1019
-1020 def __truediv__(self, other):
-1021 if isinstance(other, np.ndarray):
-1022 return 1 / (other / self)
-1023 elif hasattr(other, 'real') and hasattr(other, 'imag'):
-1024 r = other.real ** 2 + other.imag ** 2
-1025 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
-1026 else:
-1027 return CObs(self.real / other, self.imag / other)
-1028
-1029 def __rtruediv__(self, other):
-1030 r = self.real ** 2 + self.imag ** 2
-1031 if hasattr(other, 'real') and hasattr(other, 'imag'):
-1032 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
-1033 else:
-1034 return CObs(self.real * other / r, -self.imag * other / r)
-1035
-1036 def __abs__(self):
-1037 return np.sqrt(self.real**2 + self.imag**2)
-1038
-1039 def __pos__(self):
-1040 return self
-1041
-1042 def __neg__(self):
-1043 return -1 * self
-1044
-1045 def __eq__(self, other):
-1046 return self.real == other.real and self.imag == other.imag
-1047
-1048 __hash__ = None
-1049
-1050 def __str__(self):
-1051 return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
-1052
-1053 def __repr__(self):
-1054 return 'CObs[' + str(self) + ']'
-1055
-1056 def __format__(self, format_type):
-1057 if format_type == "":
-1058 significance = 2
-1059 format_type = "2"
-1060 else:
-1061 significance = int(float(format_type.replace("+", "").replace("-", "")))
-1062 return f"({self.real:{format_type}}{self.imag:+{significance}}j)"
+ 943class CObs:
+ 944 """Class for a complex valued observable."""
+ 945 __slots__ = ['_imag', '_real', 'tag']
+ 946
+ 947 def __init__(self, real, imag=0.0):
+ 948 self._real = real
+ 949 self._imag = imag
+ 950 self.tag = None
+ 951
+ 952 @property
+ 953 def real(self):
+ 954 return self._real
+ 955
+ 956 @property
+ 957 def imag(self):
+ 958 return self._imag
+ 959
+ 960 def gamma_method(self, **kwargs):
+ 961 """Executes the gamma_method for the real and the imaginary part."""
+ 962 if isinstance(self.real, Obs):
+ 963 self.real.gamma_method(**kwargs)
+ 964 if isinstance(self.imag, Obs):
+ 965 self.imag.gamma_method(**kwargs)
+ 966
+ 967 def is_zero(self):
+ 968 """Checks whether both real and imaginary part are zero within machine precision."""
+ 969 return self.real == 0.0 and self.imag == 0.0
+ 970
+ 971 def conjugate(self):
+ 972 return CObs(self.real, -self.imag)
+ 973
+ 974 def __add__(self, other):
+ 975 if isinstance(other, np.ndarray):
+ 976 return other + self
+ 977 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 978 return CObs(self.real + other.real,
+ 979 self.imag + other.imag)
+ 980 else:
+ 981 return CObs(self.real + other, self.imag)
+ 982
+ 983 def __radd__(self, y):
+ 984 return self + y
+ 985
+ 986 def __sub__(self, other):
+ 987 if isinstance(other, np.ndarray):
+ 988 return -1 * (other - self)
+ 989 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+ 990 return CObs(self.real - other.real, self.imag - other.imag)
+ 991 else:
+ 992 return CObs(self.real - other, self.imag)
+ 993
+ 994 def __rsub__(self, other):
+ 995 return -1 * (self - other)
+ 996
+ 997 def __mul__(self, other):
+ 998 if isinstance(other, np.ndarray):
+ 999 return other * self
+1000 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+1001 if all(isinstance(i, Obs) for i in [self.real, self.imag, other.real, other.imag]):
+1002 return CObs(derived_observable(lambda x, **kwargs: x[0] * x[1] - x[2] * x[3],
+1003 [self.real, other.real, self.imag, other.imag],
+1004 man_grad=[other.real.value, self.real.value, -other.imag.value, -self.imag.value]),
+1005 derived_observable(lambda x, **kwargs: x[2] * x[1] + x[0] * x[3],
+1006 [self.real, other.real, self.imag, other.imag],
+1007 man_grad=[other.imag.value, self.imag.value, other.real.value, self.real.value]))
+1008 elif getattr(other, 'imag', 0) != 0:
+1009 return CObs(self.real * other.real - self.imag * other.imag,
+1010 self.imag * other.real + self.real * other.imag)
+1011 else:
+1012 return CObs(self.real * other.real, self.imag * other.real)
+1013 else:
+1014 return CObs(self.real * other, self.imag * other)
+1015
+1016 def __rmul__(self, other):
+1017 return self * other
+1018
+1019 def __truediv__(self, other):
+1020 if isinstance(other, np.ndarray):
+1021 return 1 / (other / self)
+1022 elif hasattr(other, 'real') and hasattr(other, 'imag'):
+1023 r = other.real ** 2 + other.imag ** 2
+1024 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.imag * other.real - self.real * other.imag) / r)
+1025 else:
+1026 return CObs(self.real / other, self.imag / other)
+1027
+1028 def __rtruediv__(self, other):
+1029 r = self.real ** 2 + self.imag ** 2
+1030 if hasattr(other, 'real') and hasattr(other, 'imag'):
+1031 return CObs((self.real * other.real + self.imag * other.imag) / r, (self.real * other.imag - self.imag * other.real) / r)
+1032 else:
+1033 return CObs(self.real * other / r, -self.imag * other / r)
+1034
+1035 def __abs__(self):
+1036 return np.sqrt(self.real**2 + self.imag**2)
+1037
+1038 def __pos__(self):
+1039 return self
+1040
+1041 def __neg__(self):
+1042 return -1 * self
+1043
+1044 def __eq__(self, other):
+1045 return self.real == other.real and self.imag == other.imag
+1046
+1047 __hash__ = None
+1048
+1049 def __str__(self):
+1050 return '(' + str(self.real) + int(self.imag >= 0.0) * '+' + str(self.imag) + 'j)'
+1051
+1052 def __repr__(self):
+1053 return 'CObs[' + str(self) + ']'
+1054
+1055 def __format__(self, format_type):
+1056 if format_type == "":
+1057 significance = 2
+1058 format_type = "2"
+1059 else:
+1060 significance = int(float(format_type.replace("+", "").replace("-", "")))
+1061 return f"({self.real:{format_type}}{self.imag:+{significance}}j)"
948 def __init__(self, real, imag=0.0):
-949 self._real = real
-950 self._imag = imag
-951 self.tag = None
+ 947 def __init__(self, real, imag=0.0):
+948 self._real = real
+949 self._imag = imag
+950 self.tag = None
953 @property
-954 def real(self):
-955 return self._real
+ 952 @property
+953 def real(self):
+954 return self._real
957 @property
-958 def imag(self):
-959 return self._imag
+ 956 @property
+957 def imag(self):
+958 return self._imag
961 def gamma_method(self, **kwargs):
-962 """Executes the gamma_method for the real and the imaginary part."""
-963 if isinstance(self.real, Obs):
-964 self.real.gamma_method(**kwargs)
-965 if isinstance(self.imag, Obs):
-966 self.imag.gamma_method(**kwargs)
+ 960 def gamma_method(self, **kwargs):
+961 """Executes the gamma_method for the real and the imaginary part."""
+962 if isinstance(self.real, Obs):
+963 self.real.gamma_method(**kwargs)
+964 if isinstance(self.imag, Obs):
+965 self.imag.gamma_method(**kwargs)
968 def is_zero(self):
-969 """Checks whether both real and imaginary part are zero within machine precision."""
-970 return self.real == 0.0 and self.imag == 0.0
+ 967 def is_zero(self):
+968 """Checks whether both real and imaginary part are zero within machine precision."""
+969 return self.real == 0.0 and self.imag == 0.0
972 def conjugate(self):
-973 return CObs(self.real, -self.imag)
+ 971 def conjugate(self):
+972 return CObs(self.real, -self.imag)
1065def gamma_method(x, **kwargs):
-1066 """Vectorized version of the gamma_method applicable to lists or arrays of Obs.
-1067
-1068 See docstring of pe.Obs.gamma_method for details.
-1069 """
-1070 return np.vectorize(lambda o: o.gm(**kwargs))(x)
+ 1064def gamma_method(x, **kwargs):
+1065 """Vectorized version of the gamma_method applicable to lists or arrays of Obs.
+1066
+1067 See docstring of pe.Obs.gamma_method for details.
+1068 """
+1069 return np.vectorize(lambda o: o.gm(**kwargs))(x)
1065def gamma_method(x, **kwargs):
-1066 """Vectorized version of the gamma_method applicable to lists or arrays of Obs.
-1067
-1068 See docstring of pe.Obs.gamma_method for details.
-1069 """
-1070 return np.vectorize(lambda o: o.gm(**kwargs))(x)
+ 1064def gamma_method(x, **kwargs):
+1065 """Vectorized version of the gamma_method applicable to lists or arrays of Obs.
+1066
+1067 See docstring of pe.Obs.gamma_method for details.
+1068 """
+1069 return np.vectorize(lambda o: o.gm(**kwargs))(x)
1200def derived_observable(func, data, array_mode=False, **kwargs):
-1201 """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
-1202
-1203 Parameters
-1204 ----------
-1205 func : object
-1206 arbitrary function of the form func(data, **kwargs). For the
-1207 automatic differentiation to work, all numpy functions have to have
-1208 the autograd wrapper (use 'import autograd.numpy as anp').
-1209 data : list
-1210 list of Obs, e.g. [obs1, obs2, obs3].
-1211 num_grad : bool
-1212 if True, numerical derivatives are used instead of autograd
-1213 (default False). To control the numerical differentiation the
-1214 kwargs of numdifftools.step_generators.MaxStepGenerator
-1215 can be used.
-1216 man_grad : list
-1217 manually supply a list or an array which contains the jacobian
-1218 of func. Use cautiously, supplying the wrong derivative will
-1219 not be intercepted.
-1220
-1221 Notes
-1222 -----
-1223 For simple mathematical operations it can be practical to use anonymous
-1224 functions. For the ratio of two observables one can e.g. use
-1225
-1226 new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
-1227 """
-1228
-1229 data = np.asarray(data)
-1230 raveled_data = data.ravel()
-1231
-1232 # Workaround for matrix operations containing non Obs data
-1233 if not all(isinstance(x, Obs) for x in raveled_data):
-1234 for i in range(len(raveled_data)):
-1235 if isinstance(raveled_data[i], (int, float)):
-1236 raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
-1237
-1238 allcov = {}
-1239 for o in raveled_data:
-1240 for name in o.cov_names:
-1241 if name in allcov:
-1242 if not np.allclose(allcov[name], o.covobs[name].cov):
-1243 raise Exception(f'Inconsistent covariance matrices for {name}!')
-1244 else:
-1245 allcov[name] = o.covobs[name].cov
-1246
-1247 n_obs = len(raveled_data)
-1248 new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
-1249 new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
-1250 new_sample_names = sorted(set(new_names) - set(new_cov_names))
-1251
-1252 reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
-1253
-1254 if data.ndim == 1:
-1255 values = np.array([o.value for o in data])
-1256 else:
-1257 values = np.vectorize(lambda x: x.value)(data)
-1258
-1259 new_values = func(values, **kwargs)
-1260
-1261 multi = int(isinstance(new_values, np.ndarray))
-1262
-1263 new_r_values = {}
-1264 new_idl_d = {}
-1265 for name in new_sample_names:
-1266 idl = []
-1267 tmp_values = np.zeros(n_obs)
-1268 for i, item in enumerate(raveled_data):
-1269 tmp_values[i] = item.r_values.get(name, item.value)
-1270 tmp_idl = item.idl.get(name)
-1271 if tmp_idl is not None:
-1272 idl.append(tmp_idl)
-1273 if multi > 0:
-1274 tmp_values = np.array(tmp_values).reshape(data.shape)
-1275 new_r_values[name] = func(tmp_values, **kwargs)
-1276 new_idl_d[name] = _merge_idx(idl)
-1277
-1278 def _compute_scalefactor_missing_rep(obs):
-1279 """
-1280 Computes the scale factor that is to be multiplied with the deltas
-1281 in the case where Obs with different subsets of replica are merged.
-1282 Returns a dictionary with the scale factor for each Monte Carlo name.
-1283
-1284 Parameters
-1285 ----------
-1286 obs : Obs
-1287 The observable corresponding to the deltas that are to be scaled
-1288 """
-1289 scalef_d = {}
-1290 for mc_name in obs.mc_names:
-1291 mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')]
-1292 new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')]
-1293 if len(mc_idl_d) > 0 and len(mc_idl_d) < len(new_mc_idl_d):
-1294 scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d])
-1295 return scalef_d
-1296
-1297 if 'man_grad' in kwargs:
-1298 deriv = np.asarray(kwargs.get('man_grad'))
-1299 if new_values.shape + data.shape != deriv.shape:
-1300 raise ValueError('Manual derivative does not have correct shape.')
-1301 elif kwargs.get('num_grad') is True:
-1302 if multi > 0:
-1303 raise NotImplementedError('Multi mode currently not supported for numerical derivative')
-1304 options = {
-1305 'base_step': 0.1,
-1306 'step_ratio': 2.5}
-1307 for key in options.keys():
-1308 kwarg = kwargs.get(key)
-1309 if kwarg is not None:
-1310 options[key] = kwarg
-1311 tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
-1312 if tmp_df.size == 1:
-1313 deriv = np.array([tmp_df.real])
-1314 else:
-1315 deriv = tmp_df.real
-1316 else:
-1317 deriv = jacobian(func)(values, **kwargs)
-1318
-1319 final_result = np.zeros(new_values.shape, dtype=object)
-1320
-1321 if array_mode is True:
-1322
-1323 class _Zero_grad:
-1324 def __init__(self, N):
-1325 self.grad = np.zeros((N, 1))
-1326
-1327 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]))
-1328 d_extracted = {}
-1329 g_extracted = {}
-1330 for name in new_sample_names:
-1331 d_extracted[name] = []
-1332 ens_length = len(new_idl_d[name])
-1333 for dat in data:
-1334 d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name], _compute_scalefactor_missing_rep(o).get(name.split('|')[0], 1)) for o in dat.reshape(np.prod(dat.shape))]).reshape((*dat.shape, ens_length)))
-1335 for name in new_cov_names:
-1336 g_extracted[name] = []
-1337 zero_grad = _Zero_grad(new_covobs_lengths[name])
-1338 for dat in data:
-1339 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)))
-1340
-1341 for i_val, new_val in np.ndenumerate(new_values):
-1342 new_deltas = {}
-1343 new_grad = {}
-1344 if array_mode is True:
-1345 for name in new_sample_names:
-1346 ens_length = d_extracted[name][0].shape[-1]
-1347 new_deltas[name] = np.zeros(ens_length)
-1348 for i_dat, dat in enumerate(d_extracted[name]):
-1349 new_deltas[name] += np.tensordot(deriv[(*i_val, i_dat)], dat)
-1350 for name in new_cov_names:
-1351 new_grad[name] = 0
-1352 for i_dat, dat in enumerate(g_extracted[name]):
-1353 new_grad[name] += np.tensordot(deriv[(*i_val, i_dat)], dat)
-1354 else:
-1355 for j_obs, obs in np.ndenumerate(data):
-1356 scalef_d = _compute_scalefactor_missing_rep(obs)
-1357 for name in obs.names:
-1358 if name in obs.cov_names:
-1359 new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
-1360 else:
-1361 new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name], scalef_d.get(name.split('|')[0], 1))
-1362
-1363 new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
-1364
-1365 if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
-1366 raise ValueError('The same name has been used for deltas and covobs!')
-1367 new_samples = []
-1368 new_means = []
-1369 new_idl = []
-1370 new_names_obs = []
-1371 for name in new_names:
-1372 if name not in new_covobs:
-1373 new_samples.append(new_deltas[name])
-1374 new_idl.append(new_idl_d[name])
-1375 new_means.append(new_r_values[name][i_val])
-1376 new_names_obs.append(name)
-1377 final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
-1378 for name in new_covobs:
-1379 final_result[i_val].names.append(name)
-1380 final_result[i_val]._covobs = new_covobs
-1381 final_result[i_val]._value = new_val
-1382 final_result[i_val].reweighted = reweighted
-1383
-1384 if multi == 0:
-1385 final_result = final_result.item()
-1386
-1387 return final_result
+ 1199def derived_observable(func, data, array_mode=False, **kwargs):
+1200 """Construct a derived Obs according to func(data, **kwargs) using automatic differentiation.
+1201
+1202 Parameters
+1203 ----------
+1204 func : object
+1205 arbitrary function of the form func(data, **kwargs). For the
+1206 automatic differentiation to work, all numpy functions have to have
+1207 the autograd wrapper (use 'import autograd.numpy as anp').
+1208 data : list
+1209 list of Obs, e.g. [obs1, obs2, obs3].
+1210 num_grad : bool
+1211 if True, numerical derivatives are used instead of autograd
+1212 (default False). To control the numerical differentiation the
+1213 kwargs of numdifftools.step_generators.MaxStepGenerator
+1214 can be used.
+1215 man_grad : list
+1216 manually supply a list or an array which contains the jacobian
+1217 of func. Use cautiously, supplying the wrong derivative will
+1218 not be intercepted.
+1219
+1220 Notes
+1221 -----
+1222 For simple mathematical operations it can be practical to use anonymous
+1223 functions. For the ratio of two observables one can e.g. use
+1224
+1225 new_obs = derived_observable(lambda x: x[0] / x[1], [obs1, obs2])
+1226 """
+1227
+1228 data = np.asarray(data)
+1229 raveled_data = data.ravel()
+1230
+1231 # Workaround for matrix operations containing non Obs data
+1232 if not all(isinstance(x, Obs) for x in raveled_data):
+1233 for i in range(len(raveled_data)):
+1234 if isinstance(raveled_data[i], (int, float)):
+1235 raveled_data[i] = cov_Obs(raveled_data[i], 0.0, "###dummy_covobs###")
+1236
+1237 allcov = {}
+1238 for o in raveled_data:
+1239 for name in o.cov_names:
+1240 if name in allcov:
+1241 if not np.allclose(allcov[name], o.covobs[name].cov):
+1242 raise Exception(f'Inconsistent covariance matrices for {name}!')
+1243 else:
+1244 allcov[name] = o.covobs[name].cov
+1245
+1246 n_obs = len(raveled_data)
+1247 new_names = sorted(set([y for x in [o.names for o in raveled_data] for y in x]))
+1248 new_cov_names = sorted(set([y for x in [o.cov_names for o in raveled_data] for y in x]))
+1249 new_sample_names = sorted(set(new_names) - set(new_cov_names))
+1250
+1251 reweighted = len(list(filter(lambda o: o.reweighted is True, raveled_data))) > 0
+1252
+1253 if data.ndim == 1:
+1254 values = np.array([o.value for o in data])
+1255 else:
+1256 values = np.vectorize(lambda x: x.value)(data)
+1257
+1258 new_values = func(values, **kwargs)
+1259
+1260 multi = int(isinstance(new_values, np.ndarray))
+1261
+1262 new_r_values = {}
+1263 new_idl_d = {}
+1264 for name in new_sample_names:
+1265 idl = []
+1266 tmp_values = np.zeros(n_obs)
+1267 for i, item in enumerate(raveled_data):
+1268 tmp_values[i] = item.r_values.get(name, item.value)
+1269 tmp_idl = item.idl.get(name)
+1270 if tmp_idl is not None:
+1271 idl.append(tmp_idl)
+1272 if multi > 0:
+1273 tmp_values = np.array(tmp_values).reshape(data.shape)
+1274 new_r_values[name] = func(tmp_values, **kwargs)
+1275 new_idl_d[name] = _merge_idx(idl)
+1276
+1277 def _compute_scalefactor_missing_rep(obs):
+1278 """
+1279 Computes the scale factor that is to be multiplied with the deltas
+1280 in the case where Obs with different subsets of replica are merged.
+1281 Returns a dictionary with the scale factor for each Monte Carlo name.
+1282
+1283 Parameters
+1284 ----------
+1285 obs : Obs
+1286 The observable corresponding to the deltas that are to be scaled
+1287 """
+1288 scalef_d = {}
+1289 for mc_name in obs.mc_names:
+1290 mc_idl_d = [name for name in obs.idl if name.startswith(mc_name + '|')]
+1291 new_mc_idl_d = [name for name in new_idl_d if name.startswith(mc_name + '|')]
+1292 if len(mc_idl_d) > 0 and len(mc_idl_d) < len(new_mc_idl_d):
+1293 scalef_d[mc_name] = sum([len(new_idl_d[name]) for name in new_mc_idl_d]) / sum([len(new_idl_d[name]) for name in mc_idl_d])
+1294 return scalef_d
+1295
+1296 if 'man_grad' in kwargs:
+1297 deriv = np.asarray(kwargs.get('man_grad'))
+1298 if new_values.shape + data.shape != deriv.shape:
+1299 raise ValueError('Manual derivative does not have correct shape.')
+1300 elif kwargs.get('num_grad') is True:
+1301 if multi > 0:
+1302 raise NotImplementedError('Multi mode currently not supported for numerical derivative')
+1303 options = {
+1304 'base_step': 0.1,
+1305 'step_ratio': 2.5}
+1306 for key in options:
+1307 kwarg = kwargs.get(key)
+1308 if kwarg is not None:
+1309 options[key] = kwarg
+1310 tmp_df = nd.Gradient(func, order=4, **{k: v for k, v in options.items() if v is not None})(values, **kwargs)
+1311 if tmp_df.size == 1:
+1312 deriv = np.array([tmp_df.real])
+1313 else:
+1314 deriv = tmp_df.real
+1315 else:
+1316 deriv = jacobian(func)(values, **kwargs)
+1317
+1318 final_result = np.zeros(new_values.shape, dtype=object)
+1319
+1320 if array_mode is True:
+1321
+1322 class _Zero_grad:
+1323 def __init__(self, N):
+1324 self.grad = np.zeros((N, 1))
+1325
+1326 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]))
+1327 d_extracted = {}
+1328 g_extracted = {}
+1329 for name in new_sample_names:
+1330 d_extracted[name] = []
+1331 ens_length = len(new_idl_d[name])
+1332 for dat in data:
+1333 d_extracted[name].append(np.array([_expand_deltas_for_merge(o.deltas.get(name, np.zeros(ens_length)), o.idl.get(name, new_idl_d[name]), o.shape.get(name, ens_length), new_idl_d[name], _compute_scalefactor_missing_rep(o).get(name.split('|')[0], 1)) for o in dat.reshape(np.prod(dat.shape))]).reshape((*dat.shape, ens_length)))
+1334 for name in new_cov_names:
+1335 g_extracted[name] = []
+1336 zero_grad = _Zero_grad(new_covobs_lengths[name])
+1337 for dat in data:
+1338 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)))
+1339
+1340 for i_val, new_val in np.ndenumerate(new_values):
+1341 new_deltas = {}
+1342 new_grad = {}
+1343 if array_mode is True:
+1344 for name in new_sample_names:
+1345 ens_length = d_extracted[name][0].shape[-1]
+1346 new_deltas[name] = np.zeros(ens_length)
+1347 for i_dat, dat in enumerate(d_extracted[name]):
+1348 new_deltas[name] += np.tensordot(deriv[(*i_val, i_dat)], dat)
+1349 for name in new_cov_names:
+1350 new_grad[name] = 0
+1351 for i_dat, dat in enumerate(g_extracted[name]):
+1352 new_grad[name] += np.tensordot(deriv[(*i_val, i_dat)], dat)
+1353 else:
+1354 for j_obs, obs in np.ndenumerate(data):
+1355 scalef_d = _compute_scalefactor_missing_rep(obs)
+1356 for name in obs.names:
+1357 if name in obs.cov_names:
+1358 new_grad[name] = new_grad.get(name, 0) + deriv[i_val + j_obs] * obs.covobs[name].grad
+1359 else:
+1360 new_deltas[name] = new_deltas.get(name, 0) + deriv[i_val + j_obs] * _expand_deltas_for_merge(obs.deltas[name], obs.idl[name], obs.shape[name], new_idl_d[name], scalef_d.get(name.split('|')[0], 1))
+1361
+1362 new_covobs = {name: Covobs(0, allcov[name], name, grad=new_grad[name]) for name in new_grad}
+1363
+1364 if not set(new_covobs.keys()).isdisjoint(new_deltas.keys()):
+1365 raise ValueError('The same name has been used for deltas and covobs!')
+1366 new_samples = []
+1367 new_means = []
+1368 new_idl = []
+1369 new_names_obs = []
+1370 for name in new_names:
+1371 if name not in new_covobs:
+1372 new_samples.append(new_deltas[name])
+1373 new_idl.append(new_idl_d[name])
+1374 new_means.append(new_r_values[name][i_val])
+1375 new_names_obs.append(name)
+1376 final_result[i_val] = Obs(new_samples, new_names_obs, means=new_means, idl=new_idl)
+1377 for name in new_covobs:
+1378 final_result[i_val].names.append(name)
+1379 final_result[i_val]._covobs = new_covobs
+1380 final_result[i_val]._value = new_val
+1381 final_result[i_val].reweighted = reweighted
+1382
+1383 if multi == 0:
+1384 final_result = final_result.item()
+1385
+1386 return final_result
1419def reweight(weight, obs, **kwargs):
-1420 """Reweight a list of observables.
-1421
-1422 Parameters
-1423 ----------
-1424 weight : Obs
-1425 Reweighting factor. An Observable that has to be defined on a superset of the
-1426 configurations in obs[i].idl for all i.
-1427 obs : list
-1428 list of Obs, e.g. [obs1, obs2, obs3].
-1429 all_configs : bool
-1430 if True, the reweighted observables are normalized by the average of
-1431 the reweighting factor on all configurations in weight.idl and not
-1432 on the configurations in obs[i].idl. Default False.
-1433 """
-1434 result = []
-1435 for i in range(len(obs)):
-1436 if len(obs[i].cov_names):
-1437 raise ValueError('Error: Not possible to reweight an Obs that contains covobs!')
-1438 if not set(obs[i].names).issubset(weight.names):
-1439 raise ValueError('Error: Ensembles do not fit')
-1440 if len(obs[i].mc_names) > 1 or len(weight.mc_names) > 1:
-1441 raise ValueError('Error: Cannot reweight an Obs that contains multiple ensembles.')
-1442 for name in obs[i].names:
-1443 if not set(obs[i].idl[name]).issubset(weight.idl[name]):
-1444 raise ValueError(f'obs[{i}] has to be defined on a subset of the configs in weight.idl[{name}]!')
-1445 new_samples = []
-1446 w_deltas = {}
-1447 for name in sorted(obs[i].names):
-1448 w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
-1449 new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
-1450 tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
-1451
-1452 if kwargs.get('all_configs'):
-1453 new_weight = weight
-1454 else:
-1455 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)])
-1456
-1457 result.append(tmp_obs / new_weight)
-1458 result[-1].reweighted = True
-1459
-1460 return result
+ 1418def reweight(weight, obs, **kwargs):
+1419 """Reweight a list of observables.
+1420
+1421 Parameters
+1422 ----------
+1423 weight : Obs
+1424 Reweighting factor. An Observable that has to be defined on a superset of the
+1425 configurations in obs[i].idl for all i.
+1426 obs : list
+1427 list of Obs, e.g. [obs1, obs2, obs3].
+1428 all_configs : bool
+1429 if True, the reweighted observables are normalized by the average of
+1430 the reweighting factor on all configurations in weight.idl and not
+1431 on the configurations in obs[i].idl. Default False.
+1432 """
+1433 result = []
+1434 for i in range(len(obs)):
+1435 if len(obs[i].cov_names):
+1436 raise ValueError('Error: Not possible to reweight an Obs that contains covobs!')
+1437 if not set(obs[i].names).issubset(weight.names):
+1438 raise ValueError('Error: Ensembles do not fit')
+1439 if len(obs[i].mc_names) > 1 or len(weight.mc_names) > 1:
+1440 raise ValueError('Error: Cannot reweight an Obs that contains multiple ensembles.')
+1441 for name in obs[i].names:
+1442 if not set(obs[i].idl[name]).issubset(weight.idl[name]):
+1443 raise ValueError(f'obs[{i}] has to be defined on a subset of the configs in weight.idl[{name}]!')
+1444 new_samples = []
+1445 w_deltas = {}
+1446 for name in sorted(obs[i].names):
+1447 w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name])
+1448 new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name]))
+1449 tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)])
+1450
+1451 if kwargs.get('all_configs'):
+1452 new_weight = weight
+1453 else:
+1454 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)])
+1455
+1456 result.append(tmp_obs / new_weight)
+1457 result[-1].reweighted = True
+1458
+1459 return result
1463def correlate(obs_a, obs_b):
-1464 """Correlate two observables.
-1465
-1466 Parameters
-1467 ----------
-1468 obs_a : Obs
-1469 First observable
-1470 obs_b : Obs
-1471 Second observable
-1472
-1473 Notes
-1474 -----
-1475 Keep in mind to only correlate primary observables which have not been reweighted
-1476 yet. The reweighting has to be applied after correlating the observables.
-1477 Only works if a single ensemble is present in the Obs.
-1478 Currently only works if ensemble content is identical (this is not strictly necessary).
-1479 """
-1480
-1481 if len(obs_a.mc_names) > 1 or len(obs_b.mc_names) > 1:
-1482 raise ValueError('Error: Cannot correlate Obs that contain multiple ensembles.')
-1483 if sorted(obs_a.names) != sorted(obs_b.names):
-1484 raise ValueError(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
-1485 if len(obs_a.cov_names) or len(obs_b.cov_names):
-1486 raise ValueError('Error: Not possible to correlate Obs that contain covobs!')
-1487 for name in obs_a.names:
-1488 if obs_a.shape[name] != obs_b.shape[name]:
-1489 raise ValueError('Shapes of ensemble', name, 'do not fit')
-1490 if obs_a.idl[name] != obs_b.idl[name]:
-1491 raise ValueError('idl of ensemble', name, 'do not fit')
-1492
-1493 if obs_a.reweighted is True:
-1494 warnings.warn("The first observable is already reweighted.", RuntimeWarning, stacklevel=2)
-1495 if obs_b.reweighted is True:
-1496 warnings.warn("The second observable is already reweighted.", RuntimeWarning, stacklevel=2)
-1497
-1498 new_samples = []
-1499 new_idl = []
-1500 for name in sorted(obs_a.names):
-1501 new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
-1502 new_idl.append(obs_a.idl[name])
-1503
-1504 o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
-1505 o.reweighted = obs_a.reweighted or obs_b.reweighted
-1506 return o
+ 1462def correlate(obs_a, obs_b):
+1463 """Correlate two observables.
+1464
+1465 Parameters
+1466 ----------
+1467 obs_a : Obs
+1468 First observable
+1469 obs_b : Obs
+1470 Second observable
+1471
+1472 Notes
+1473 -----
+1474 Keep in mind to only correlate primary observables which have not been reweighted
+1475 yet. The reweighting has to be applied after correlating the observables.
+1476 Only works if a single ensemble is present in the Obs.
+1477 Currently only works if ensemble content is identical (this is not strictly necessary).
+1478 """
+1479
+1480 if len(obs_a.mc_names) > 1 or len(obs_b.mc_names) > 1:
+1481 raise ValueError('Error: Cannot correlate Obs that contain multiple ensembles.')
+1482 if sorted(obs_a.names) != sorted(obs_b.names):
+1483 raise ValueError(f"Ensembles do not fit {set(sorted(obs_a.names)) ^ set(sorted(obs_b.names))}")
+1484 if len(obs_a.cov_names) or len(obs_b.cov_names):
+1485 raise ValueError('Error: Not possible to correlate Obs that contain covobs!')
+1486 for name in obs_a.names:
+1487 if obs_a.shape[name] != obs_b.shape[name]:
+1488 raise ValueError('Shapes of ensemble', name, 'do not fit')
+1489 if obs_a.idl[name] != obs_b.idl[name]:
+1490 raise ValueError('idl of ensemble', name, 'do not fit')
+1491
+1492 if obs_a.reweighted is True:
+1493 warnings.warn("The first observable is already reweighted.", RuntimeWarning, stacklevel=2)
+1494 if obs_b.reweighted is True:
+1495 warnings.warn("The second observable is already reweighted.", RuntimeWarning, stacklevel=2)
+1496
+1497 new_samples = []
+1498 new_idl = []
+1499 for name in sorted(obs_a.names):
+1500 new_samples.append((obs_a.deltas[name] + obs_a.r_values[name]) * (obs_b.deltas[name] + obs_b.r_values[name]))
+1501 new_idl.append(obs_a.idl[name])
+1502
+1503 o = Obs(new_samples, sorted(obs_a.names), idl=new_idl)
+1504 o.reweighted = obs_a.reweighted or obs_b.reweighted
+1505 return o
1509def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
-1510 r'''Calculates the error covariance matrix of a set of observables.
-1511
-1512 WARNING: This function should be used with care, especially for observables with support on multiple
-1513 ensembles with differing autocorrelations. See the notes below for details.
-1514
-1515 The gamma method has to be applied first to all observables.
-1516
-1517 Parameters
-1518 ----------
-1519 obs : list or numpy.ndarray
-1520 List or one dimensional array of Obs
-1521 visualize : bool
-1522 If True plots the corresponding normalized correlation matrix (default False).
-1523 correlation : bool
-1524 If True the correlation matrix instead of the error covariance matrix is returned (default False).
-1525 smooth : None or int
-1526 If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
-1527 smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
-1528 largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
-1529 small ones.
-1530
-1531 Notes
-1532 -----
-1533 The error covariance is defined such that it agrees with the squared standard error for two identical observables
-1534 $$\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$$
-1535 in the absence of autocorrelation.
-1536 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
-1537 $$\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.
-1538 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.
-1539 $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
-1540 This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
-1541 '''
-1542
-1543 length = len(obs)
-1544
-1545 max_samples = np.max([o.N for o in obs])
-1546 if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
-1547 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, stacklevel=2)
-1548
-1549 cov = np.zeros((length, length))
-1550 for i in range(length):
-1551 for j in range(i, length):
-1552 cov[i, j] = _covariance_element(obs[i], obs[j])
-1553 cov = cov + cov.T - np.diag(np.diag(cov))
-1554
-1555 corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
-1556
-1557 if isinstance(smooth, int):
-1558 corr = _smooth_eigenvalues(corr, smooth)
-1559
-1560 if visualize:
-1561 plt.matshow(corr, vmin=-1, vmax=1)
-1562 plt.set_cmap('RdBu')
-1563 plt.colorbar()
-1564 plt.draw()
-1565
-1566 if correlation is True:
-1567 return corr
-1568
-1569 errors = [o.dvalue for o in obs]
-1570 cov = np.diag(errors) @ corr @ np.diag(errors)
-1571
-1572 eigenvalues = np.linalg.eigh(cov)[0]
-1573 if not np.all(eigenvalues >= 0):
-1574 warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning, stacklevel=2)
-1575
-1576 return cov
+ 1508def covariance(obs, visualize=False, correlation=False, smooth=None, **kwargs):
+1509 r'''Calculates the error covariance matrix of a set of observables.
+1510
+1511 WARNING: This function should be used with care, especially for observables with support on multiple
+1512 ensembles with differing autocorrelations. See the notes below for details.
+1513
+1514 The gamma method has to be applied first to all observables.
+1515
+1516 Parameters
+1517 ----------
+1518 obs : list or numpy.ndarray
+1519 List or one dimensional array of Obs
+1520 visualize : bool
+1521 If True plots the corresponding normalized correlation matrix (default False).
+1522 correlation : bool
+1523 If True the correlation matrix instead of the error covariance matrix is returned (default False).
+1524 smooth : None or int
+1525 If smooth is an integer 'E' between 2 and the dimension of the matrix minus 1 the eigenvalue
+1526 smoothing procedure of hep-lat/9412087 is applied to the correlation matrix which leaves the
+1527 largest E eigenvalues essentially unchanged and smoothes the smaller eigenvalues to avoid extremely
+1528 small ones.
+1529
+1530 Notes
+1531 -----
+1532 The error covariance is defined such that it agrees with the squared standard error for two identical observables
+1533 $$\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$$
+1534 in the absence of autocorrelation.
+1535 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
+1536 $$\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.
+1537 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.
+1538 $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
+1539 This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
+1540 '''
+1541
+1542 length = len(obs)
+1543
+1544 max_samples = np.max([o.N for o in obs])
+1545 if max_samples <= length and not [item for sublist in [o.cov_names for o in obs] for item in sublist]:
+1546 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, stacklevel=2)
+1547
+1548 cov = np.zeros((length, length))
+1549 for i in range(length):
+1550 for j in range(i, length):
+1551 cov[i, j] = _covariance_element(obs[i], obs[j])
+1552 cov = cov + cov.T - np.diag(np.diag(cov))
+1553
+1554 corr = np.diag(1 / np.sqrt(np.diag(cov))) @ cov @ np.diag(1 / np.sqrt(np.diag(cov)))
+1555
+1556 if isinstance(smooth, int):
+1557 corr = _smooth_eigenvalues(corr, smooth)
+1558
+1559 if visualize:
+1560 plt.matshow(corr, vmin=-1, vmax=1)
+1561 plt.set_cmap('RdBu')
+1562 plt.colorbar()
+1563 plt.draw()
+1564
+1565 if correlation is True:
+1566 return corr
+1567
+1568 errors = [o.dvalue for o in obs]
+1569 cov = np.diag(errors) @ corr @ np.diag(errors)
+1570
+1571 eigenvalues = np.linalg.eigh(cov)[0]
+1572 if not np.all(eigenvalues >= 0):
+1573 warnings.warn("Covariance matrix is not positive semi-definite (Eigenvalues: " + str(eigenvalues) + ")", RuntimeWarning, stacklevel=2)
+1574
+1575 return cov
1579def invert_corr_cov_cholesky(corr, inverrdiag):
-1580 """Constructs a lower triangular matrix `chol` via the Cholesky decomposition of the correlation matrix `corr`
-1581 and then returns the inverse covariance matrix `chol_inv` as a lower triangular matrix by solving `chol * x = inverrdiag`.
-1582
-1583 Parameters
-1584 ----------
-1585 corr : np.ndarray
-1586 correlation matrix
-1587 inverrdiag : np.ndarray
-1588 diagonal matrix, the entries are the inverse errors of the data points considered
-1589 """
-1590
-1591 condn = np.linalg.cond(corr)
-1592 if condn > 0.1 / np.finfo(float).eps:
-1593 raise ValueError(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
-1594 if condn > 1e13:
-1595 warnings.warn(f"Correlation matrix may be ill-conditioned, condition number: {{{condn:1.2e}}}", RuntimeWarning, stacklevel=2)
-1596 chol = np.linalg.cholesky(corr)
-1597 chol_inv = scipy.linalg.solve_triangular(chol, inverrdiag, lower=True)
-1598
-1599 return chol_inv
+ 1578def invert_corr_cov_cholesky(corr, inverrdiag):
+1579 """Constructs a lower triangular matrix `chol` via the Cholesky decomposition of the correlation matrix `corr`
+1580 and then returns the inverse covariance matrix `chol_inv` as a lower triangular matrix by solving `chol * x = inverrdiag`.
+1581
+1582 Parameters
+1583 ----------
+1584 corr : np.ndarray
+1585 correlation matrix
+1586 inverrdiag : np.ndarray
+1587 diagonal matrix, the entries are the inverse errors of the data points considered
+1588 """
+1589
+1590 condn = np.linalg.cond(corr)
+1591 if condn > 0.1 / np.finfo(float).eps:
+1592 raise ValueError(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
+1593 if condn > 1e13:
+1594 warnings.warn(f"Correlation matrix may be ill-conditioned, condition number: {{{condn:1.2e}}}", RuntimeWarning, stacklevel=2)
+1595 chol = np.linalg.cholesky(corr)
+1596 chol_inv = scipy.linalg.solve_triangular(chol, inverrdiag, lower=True)
+1597
+1598 return chol_inv
1602def sort_corr(corr, kl, yd):
-1603 """ Reorders a correlation matrix to match the alphabetical order of its underlying y data.
-1604
-1605 The ordering of the input correlation matrix `corr` is given by the list of keys `kl`.
-1606 The input dictionary `yd` (with the same keys `kl`) must contain the corresponding y data
-1607 that the correlation matrix is based on.
-1608 This function sorts the list of keys `kl` alphabetically and sorts the matrix `corr`
-1609 according to this alphabetical order such that the sorted matrix `corr_sorted` corresponds
-1610 to the y data `yd` when arranged in an alphabetical order by its keys.
-1611
-1612 Parameters
-1613 ----------
-1614 corr : np.ndarray
-1615 A square correlation matrix constructed using the order of the y data specified by `kl`.
-1616 The dimensions of `corr` should match the total number of y data points in `yd` combined.
-1617 kl : list of str
-1618 A list of keys that denotes the order in which the y data from `yd` was used to build the
-1619 input correlation matrix `corr`.
-1620 yd : dict of list
-1621 A dictionary where each key corresponds to a unique identifier, and its value is a list of
-1622 y data points. The total number of y data points across all keys must match the dimensions
-1623 of `corr`. The lists in the dictionary can be lists of Obs.
-1624
-1625 Returns
-1626 -------
-1627 np.ndarray
-1628 A new, sorted correlation matrix that corresponds to the y data from `yd` when arranged alphabetically by its keys.
-1629
-1630 Example
-1631 -------
-1632 >>> import numpy as np
-1633 >>> import pyerrors as pe
-1634 >>> corr = np.array([[1, 0.2, 0.3], [0.2, 1, 0.4], [0.3, 0.4, 1]])
-1635 >>> kl = ['b', 'a']
-1636 >>> yd = {'a': [1, 2], 'b': [3]}
-1637 >>> sorted_corr = pe.obs.sort_corr(corr, kl, yd)
-1638 >>> print(sorted_corr)
-1639 array([[1. , 0.3, 0.4],
-1640 [0.3, 1. , 0.2],
-1641 [0.4, 0.2, 1. ]])
-1642
-1643 """
-1644 kl_sorted = sorted(kl)
-1645
-1646 posd = {}
-1647 ofs = 0
-1648 for _ki, k in enumerate(kl):
-1649 posd[k] = [i + ofs for i in range(len(yd[k]))]
-1650 ofs += len(posd[k])
-1651
-1652 mapping = []
-1653 for k in kl_sorted:
-1654 for i in range(len(yd[k])):
-1655 mapping.append(posd[k][i])
-1656
-1657 corr_sorted = np.zeros_like(corr)
-1658 for i in range(corr.shape[0]):
-1659 for j in range(corr.shape[0]):
-1660 corr_sorted[i][j] = corr[mapping[i]][mapping[j]]
-1661
-1662 return corr_sorted
+ 1601def sort_corr(corr, kl, yd):
+1602 """ Reorders a correlation matrix to match the alphabetical order of its underlying y data.
+1603
+1604 The ordering of the input correlation matrix `corr` is given by the list of keys `kl`.
+1605 The input dictionary `yd` (with the same keys `kl`) must contain the corresponding y data
+1606 that the correlation matrix is based on.
+1607 This function sorts the list of keys `kl` alphabetically and sorts the matrix `corr`
+1608 according to this alphabetical order such that the sorted matrix `corr_sorted` corresponds
+1609 to the y data `yd` when arranged in an alphabetical order by its keys.
+1610
+1611 Parameters
+1612 ----------
+1613 corr : np.ndarray
+1614 A square correlation matrix constructed using the order of the y data specified by `kl`.
+1615 The dimensions of `corr` should match the total number of y data points in `yd` combined.
+1616 kl : list of str
+1617 A list of keys that denotes the order in which the y data from `yd` was used to build the
+1618 input correlation matrix `corr`.
+1619 yd : dict of list
+1620 A dictionary where each key corresponds to a unique identifier, and its value is a list of
+1621 y data points. The total number of y data points across all keys must match the dimensions
+1622 of `corr`. The lists in the dictionary can be lists of Obs.
+1623
+1624 Returns
+1625 -------
+1626 np.ndarray
+1627 A new, sorted correlation matrix that corresponds to the y data from `yd` when arranged alphabetically by its keys.
+1628
+1629 Example
+1630 -------
+1631 >>> import numpy as np
+1632 >>> import pyerrors as pe
+1633 >>> corr = np.array([[1, 0.2, 0.3], [0.2, 1, 0.4], [0.3, 0.4, 1]])
+1634 >>> kl = ['b', 'a']
+1635 >>> yd = {'a': [1, 2], 'b': [3]}
+1636 >>> sorted_corr = pe.obs.sort_corr(corr, kl, yd)
+1637 >>> print(sorted_corr)
+1638 array([[1. , 0.3, 0.4],
+1639 [0.3, 1. , 0.2],
+1640 [0.4, 0.2, 1. ]])
+1641
+1642 """
+1643 kl_sorted = sorted(kl)
+1644
+1645 posd = {}
+1646 ofs = 0
+1647 for _ki, k in enumerate(kl):
+1648 posd[k] = [i + ofs for i in range(len(yd[k]))]
+1649 ofs += len(posd[k])
+1650
+1651 mapping = []
+1652 for k in kl_sorted:
+1653 for i in range(len(yd[k])):
+1654 mapping.append(posd[k][i])
+1655
+1656 corr_sorted = np.zeros_like(corr)
+1657 for i in range(corr.shape[0]):
+1658 for j in range(corr.shape[0]):
+1659 corr_sorted[i][j] = corr[mapping[i]][mapping[j]]
+1660
+1661 return corr_sorted
corr. The lists in the dictionary can be lists of Obs.
1742def import_jackknife(jacks, name, idl=None):
-1743 """Imports jackknife samples and returns an Obs
-1744
-1745 Parameters
-1746 ----------
-1747 jacks : numpy.ndarray
-1748 numpy array containing the mean value as zeroth entry and
-1749 the N jackknife samples as first to Nth entry.
-1750 name : str
-1751 name of the ensemble the samples are defined on.
-1752 """
-1753 length = len(jacks) - 1
-1754 prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
-1755 samples = jacks[1:] @ prj
-1756 mean = np.mean(samples)
-1757 new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
-1758 new_obs._value = jacks[0]
-1759 return new_obs
+ 1741def import_jackknife(jacks, name, idl=None):
+1742 """Imports jackknife samples and returns an Obs
+1743
+1744 Parameters
+1745 ----------
+1746 jacks : numpy.ndarray
+1747 numpy array containing the mean value as zeroth entry and
+1748 the N jackknife samples as first to Nth entry.
+1749 name : str
+1750 name of the ensemble the samples are defined on.
+1751 """
+1752 length = len(jacks) - 1
+1753 prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
+1754 samples = jacks[1:] @ prj
+1755 mean = np.mean(samples)
+1756 new_obs = Obs([samples - mean], [name], idl=idl, means=[mean])
+1757 new_obs._value = jacks[0]
+1758 return new_obs
1762def import_bootstrap(boots, name, random_numbers):
-1763 """Imports bootstrap samples and returns an Obs
-1764
-1765 Parameters
-1766 ----------
-1767 boots : numpy.ndarray
-1768 numpy array containing the mean value as zeroth entry and
-1769 the N bootstrap samples as first to Nth entry.
-1770 name : str
-1771 name of the ensemble the samples are defined on.
-1772 random_numbers : np.ndarray
-1773 Array of shape (samples, length) containing the random numbers to generate the bootstrap samples,
-1774 where samples is the number of bootstrap samples and length is the length of the original Monte Carlo
-1775 chain to be reconstructed.
-1776 """
-1777 samples, length = random_numbers.shape
-1778 if samples != len(boots) - 1:
-1779 raise ValueError("Random numbers do not have the correct shape.")
-1780
-1781 if samples < length:
-1782 raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.")
-1783
-1784 proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
-1785
-1786 samples = scipy.linalg.lstsq(proj, boots[1:])[0]
-1787 ret = Obs([samples], [name])
-1788 ret._value = boots[0]
-1789 return ret
+ 1761def import_bootstrap(boots, name, random_numbers):
+1762 """Imports bootstrap samples and returns an Obs
+1763
+1764 Parameters
+1765 ----------
+1766 boots : numpy.ndarray
+1767 numpy array containing the mean value as zeroth entry and
+1768 the N bootstrap samples as first to Nth entry.
+1769 name : str
+1770 name of the ensemble the samples are defined on.
+1771 random_numbers : np.ndarray
+1772 Array of shape (samples, length) containing the random numbers to generate the bootstrap samples,
+1773 where samples is the number of bootstrap samples and length is the length of the original Monte Carlo
+1774 chain to be reconstructed.
+1775 """
+1776 samples, length = random_numbers.shape
+1777 if samples != len(boots) - 1:
+1778 raise ValueError("Random numbers do not have the correct shape.")
+1779
+1780 if samples < length:
+1781 raise ValueError("Obs can't be reconstructed if there are fewer bootstrap samples than Monte Carlo data points.")
+1782
+1783 proj = np.vstack([np.bincount(o, minlength=length) for o in random_numbers]) / length
+1784
+1785 samples = scipy.linalg.lstsq(proj, boots[1:])[0]
+1786 ret = Obs([samples], [name])
+1787 ret._value = boots[0]
+1788 return ret
1792def merge_obs(list_of_obs):
-1793 """Combine all observables in list_of_obs into one new observable.
-1794 This allows to merge Obs that have been computed on multiple replica
-1795 of the same ensemble.
-1796 If you like to merge Obs that are based on several ensembles, please
-1797 average them yourself.
-1798
-1799 Parameters
-1800 ----------
-1801 list_of_obs : list
-1802 list of the Obs object to be combined
-1803
-1804 Notes
-1805 -----
-1806 It is not possible to combine obs which are based on the same replicum
-1807 """
-1808 replist = [item for obs in list_of_obs for item in obs.names]
-1809 if (len(replist) == len(set(replist))) is False:
-1810 raise ValueError(f'list_of_obs contains duplicate replica: {replist!s}')
-1811 if any([len(o.cov_names) for o in list_of_obs]):
-1812 raise ValueError('Not possible to merge data that contains covobs!')
-1813 new_dict = {}
-1814 idl_dict = {}
-1815 for o in list_of_obs:
-1816 new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
-1817 for key in set(o.deltas) | set(o.r_values)})
-1818 idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
-1819
-1820 names = sorted(new_dict.keys())
-1821 o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
-1822 o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
-1823 return o
+ 1791def merge_obs(list_of_obs):
+1792 """Combine all observables in list_of_obs into one new observable.
+1793 This allows to merge Obs that have been computed on multiple replica
+1794 of the same ensemble.
+1795 If you like to merge Obs that are based on several ensembles, please
+1796 average them yourself.
+1797
+1798 Parameters
+1799 ----------
+1800 list_of_obs : list
+1801 list of the Obs object to be combined
+1802
+1803 Notes
+1804 -----
+1805 It is not possible to combine obs which are based on the same replicum
+1806 """
+1807 replist = [item for obs in list_of_obs for item in obs.names]
+1808 if (len(replist) == len(set(replist))) is False:
+1809 raise ValueError(f'list_of_obs contains duplicate replica: {replist!s}')
+1810 if any([len(o.cov_names) for o in list_of_obs]):
+1811 raise ValueError('Not possible to merge data that contains covobs!')
+1812 new_dict = {}
+1813 idl_dict = {}
+1814 for o in list_of_obs:
+1815 new_dict.update({key: o.deltas.get(key, 0) + o.r_values.get(key, 0)
+1816 for key in set(o.deltas) | set(o.r_values)})
+1817 idl_dict.update({key: o.idl.get(key, 0) for key in set(o.deltas)})
+1818
+1819 names = sorted(new_dict.keys())
+1820 o = Obs([new_dict[name] for name in names], names, idl=[idl_dict[name] for name in names])
+1821 o.reweighted = np.max([oi.reweighted for oi in list_of_obs])
+1822 return o
1826def cov_Obs(means, cov, name, grad=None):
-1827 """Create an Obs based on mean(s) and a covariance matrix
-1828
-1829 Parameters
-1830 ----------
-1831 mean : list of floats or float
-1832 N mean value(s) of the new Obs
-1833 cov : list or array
-1834 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
-1835 name : str
-1836 identifier for the covariance matrix
-1837 grad : list or array
-1838 Gradient of the Covobs wrt. the means belonging to cov.
-1839 """
-1840
-1841 def covobs_to_obs(co):
-1842 """Make an Obs out of a Covobs
-1843
-1844 Parameters
-1845 ----------
-1846 co : Covobs
-1847 Covobs to be embedded into the Obs
-1848 """
-1849 o = Obs([], [], means=[])
-1850 o._value = co.value
-1851 o.names.append(co.name)
-1852 o._covobs[co.name] = co
-1853 o._dvalue = np.sqrt(co.errsq())
-1854 return o
-1855
-1856 ol = []
-1857 if isinstance(means, (float, int)):
-1858 means = [means]
-1859
-1860 for i in range(len(means)):
-1861 ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
-1862 if ol[0].covobs[name].N != len(means):
-1863 raise ValueError(f'You have to provide {ol[0].N} mean values!')
-1864 if len(ol) == 1:
-1865 return ol[0]
-1866 return ol
+ 1825def cov_Obs(means, cov, name, grad=None):
+1826 """Create an Obs based on mean(s) and a covariance matrix
+1827
+1828 Parameters
+1829 ----------
+1830 mean : list of floats or float
+1831 N mean value(s) of the new Obs
+1832 cov : list or array
+1833 2d (NxN) Covariance matrix, 1d diagonal entries or 0d covariance
+1834 name : str
+1835 identifier for the covariance matrix
+1836 grad : list or array
+1837 Gradient of the Covobs wrt. the means belonging to cov.
+1838 """
+1839
+1840 def covobs_to_obs(co):
+1841 """Make an Obs out of a Covobs
+1842
+1843 Parameters
+1844 ----------
+1845 co : Covobs
+1846 Covobs to be embedded into the Obs
+1847 """
+1848 o = Obs([], [], means=[])
+1849 o._value = co.value
+1850 o.names.append(co.name)
+1851 o._covobs[co.name] = co
+1852 o._dvalue = np.sqrt(co.errsq())
+1853 return o
+1854
+1855 ol = []
+1856 if isinstance(means, (float, int)):
+1857 means = [means]
+1858
+1859 for i in range(len(means)):
+1860 ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
+1861 if ol[0].covobs[name].N != len(means):
+1862 raise ValueError(f'You have to provide {ol[0].N} mean values!')
+1863 if len(ol) == 1:
+1864 return ol[0]
+1865 return ol