From 11766b5762f0ee5e681308add8c934fe75e89ce0 Mon Sep 17 00:00:00 2001 From: fjosw Date: Tue, 28 Jun 2022 13:42:42 +0000 Subject: [PATCH] Documentation updated --- docs/pyerrors/correlators.html | 3232 ++++++++++++++++---------------- 1 file changed, 1622 insertions(+), 1610 deletions(-) diff --git a/docs/pyerrors/correlators.html b/docs/pyerrors/correlators.html index 2cfd37dc..3a73ce43 100644 --- a/docs/pyerrors/correlators.html +++ b/docs/pyerrors/correlators.html @@ -818,658 +818,662 @@ 604 for t in range(self.T - 1): 605 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): 606 newcontent.append(None) - 607 else: - 608 newcontent.append(self.content[t] / self.content[t + 1]) - 609 if(all([x is None for x in newcontent])): - 610 raise Exception('m_eff is undefined at all timeslices') - 611 - 612 return np.log(Corr(newcontent, padding=[0, 1])) + 607 elif self.content[t][0].value / self.content[t + 1][0].value < 0: + 608 newcontent.append(None) + 609 else: + 610 newcontent.append(self.content[t] / self.content[t + 1]) + 611 if(all([x is None for x in newcontent])): + 612 raise Exception('m_eff is undefined at all timeslices') 613 - 614 elif variant in ['periodic', 'cosh', 'sinh']: - 615 if variant in ['periodic', 'cosh']: - 616 func = anp.cosh - 617 else: - 618 func = anp.sinh - 619 - 620 def root_function(x, d): - 621 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d - 622 - 623 newcontent = [] - 624 for t in range(self.T - 1): - 625 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): - 626 newcontent.append(None) - 627 # Fill the two timeslices in the middle of the lattice with their predecessors - 628 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: - 629 newcontent.append(newcontent[-1]) - 630 else: - 631 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) - 632 if(all([x is None for x in newcontent])): - 633 raise Exception('m_eff is undefined at all timeslices') - 634 - 635 return Corr(newcontent, padding=[0, 1]) - 636 - 637 elif variant == 'arccosh': - 638 newcontent = [] - 639 for t in range(1, self.T - 1): - 640 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): - 641 newcontent.append(None) - 642 else: - 643 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) - 644 if(all([x is None for x in newcontent])): - 645 raise Exception("m_eff is undefined at all timeslices") - 646 return np.arccosh(Corr(newcontent, padding=[1, 1])) - 647 - 648 else: - 649 raise Exception('Unknown variant.') - 650 - 651 def fit(self, function, fitrange=None, silent=False, **kwargs): - 652 r'''Fits function to the data - 653 - 654 Parameters - 655 ---------- - 656 function : obj - 657 function to fit to the data. See fits.least_squares for details. - 658 fitrange : list - 659 Two element list containing the timeslices on which the fit is supposed to start and stop. - 660 Caution: This range is inclusive as opposed to standard python indexing. - 661 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. - 662 If not specified, self.prange or all timeslices are used. - 663 silent : bool - 664 Decides whether output is printed to the standard output. - 665 ''' - 666 if self.N != 1: - 667 raise Exception("Correlator must be projected before fitting") - 668 - 669 if fitrange is None: - 670 if self.prange: - 671 fitrange = self.prange - 672 else: - 673 fitrange = [0, self.T - 1] - 674 else: - 675 if not isinstance(fitrange, list): - 676 raise Exception("fitrange has to be a list with two elements") - 677 if len(fitrange) != 2: - 678 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") - 679 - 680 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] - 681 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] - 682 result = least_squares(xs, ys, function, silent=silent, **kwargs) - 683 return result - 684 - 685 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): - 686 """ Extract a plateau value from a Corr object - 687 - 688 Parameters - 689 ---------- - 690 plateau_range : list - 691 list with two entries, indicating the first and the last timeslice - 692 of the plateau region. - 693 method : str - 694 method to extract the plateau. - 695 'fit' fits a constant to the plateau region - 696 'avg', 'average' or 'mean' just average over the given timeslices. - 697 auto_gamma : bool - 698 apply gamma_method with default parameters to the Corr. Defaults to None - 699 """ - 700 if not plateau_range: - 701 if self.prange: - 702 plateau_range = self.prange - 703 else: - 704 raise Exception("no plateau range provided") - 705 if self.N != 1: - 706 raise Exception("Correlator must be projected before getting a plateau.") - 707 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): - 708 raise Exception("plateau is undefined at all timeslices in plateaurange.") - 709 if auto_gamma: - 710 self.gamma_method() - 711 if method == "fit": - 712 def const_func(a, t): - 713 return a[0] - 714 return self.fit(const_func, plateau_range)[0] - 715 elif method in ["avg", "average", "mean"]: - 716 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) - 717 return returnvalue - 718 - 719 else: - 720 raise Exception("Unsupported plateau method: " + method) - 721 - 722 def set_prange(self, prange): - 723 """Sets the attribute prange of the Corr object.""" - 724 if not len(prange) == 2: - 725 raise Exception("prange must be a list or array with two values") - 726 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): - 727 raise Exception("Start and end point must be integers") - 728 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): - 729 raise Exception("Start and end point must define a range in the interval 0,T") - 730 - 731 self.prange = prange - 732 return - 733 - 734 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None): - 735 """Plots the correlator using the tag of the correlator as label if available. - 736 - 737 Parameters - 738 ---------- - 739 x_range : list - 740 list of two values, determining the range of the x-axis e.g. [4, 8]. - 741 comp : Corr or list of Corr - 742 Correlator or list of correlators which are plotted for comparison. - 743 The tags of these correlators are used as labels if available. - 744 logscale : bool - 745 Sets y-axis to logscale. - 746 plateau : Obs - 747 Plateau value to be visualized in the figure. - 748 fit_res : Fit_result - 749 Fit_result object to be visualized. - 750 ylabel : str - 751 Label for the y-axis. - 752 save : str - 753 path to file in which the figure should be saved. - 754 auto_gamma : bool - 755 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. - 756 hide_sigma : float - 757 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. - 758 references : list - 759 List of floating point values that are displayed as horizontal lines for reference. - 760 title : string - 761 Optional title of the figure. - 762 """ - 763 if self.N != 1: - 764 raise Exception("Correlator must be projected before plotting") - 765 - 766 if auto_gamma: - 767 self.gamma_method() - 768 - 769 if x_range is None: - 770 x_range = [0, self.T - 1] - 771 - 772 fig = plt.figure() - 773 ax1 = fig.add_subplot(111) - 774 - 775 x, y, y_err = self.plottable() - 776 if hide_sigma: - 777 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 - 778 else: - 779 hide_from = None - 780 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) - 781 if logscale: - 782 ax1.set_yscale('log') - 783 else: - 784 if y_range is None: - 785 try: - 786 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)]) - 787 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)]) - 788 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) - 789 except Exception: - 790 pass - 791 else: - 792 ax1.set_ylim(y_range) - 793 if comp: - 794 if isinstance(comp, (Corr, list)): - 795 for corr in comp if isinstance(comp, list) else [comp]: - 796 if auto_gamma: - 797 corr.gamma_method() - 798 x, y, y_err = corr.plottable() - 799 if hide_sigma: - 800 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 - 801 else: - 802 hide_from = None - 803 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) - 804 else: - 805 raise Exception("'comp' must be a correlator or a list of correlators.") - 806 - 807 if plateau: - 808 if isinstance(plateau, Obs): - 809 if auto_gamma: - 810 plateau.gamma_method() - 811 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) - 812 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') - 813 else: - 814 raise Exception("'plateau' must be an Obs") - 815 - 816 if references: - 817 if isinstance(references, list): - 818 for ref in references: - 819 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') - 820 else: - 821 raise Exception("'references' must be a list of floating pint values.") - 822 - 823 if self.prange: - 824 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') - 825 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') + 614 return np.log(Corr(newcontent, padding=[0, 1])) + 615 + 616 elif variant in ['periodic', 'cosh', 'sinh']: + 617 if variant in ['periodic', 'cosh']: + 618 func = anp.cosh + 619 else: + 620 func = anp.sinh + 621 + 622 def root_function(x, d): + 623 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d + 624 + 625 newcontent = [] + 626 for t in range(self.T - 1): + 627 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): + 628 newcontent.append(None) + 629 # Fill the two timeslices in the middle of the lattice with their predecessors + 630 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: + 631 newcontent.append(newcontent[-1]) + 632 elif self.content[t][0].value / self.content[t + 1][0].value < 0: + 633 newcontent.append(None) + 634 else: + 635 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) + 636 if(all([x is None for x in newcontent])): + 637 raise Exception('m_eff is undefined at all timeslices') + 638 + 639 return Corr(newcontent, padding=[0, 1]) + 640 + 641 elif variant == 'arccosh': + 642 newcontent = [] + 643 for t in range(1, self.T - 1): + 644 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): + 645 newcontent.append(None) + 646 else: + 647 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) + 648 if(all([x is None for x in newcontent])): + 649 raise Exception("m_eff is undefined at all timeslices") + 650 return np.arccosh(Corr(newcontent, padding=[1, 1])) + 651 + 652 else: + 653 raise Exception('Unknown variant.') + 654 + 655 def fit(self, function, fitrange=None, silent=False, **kwargs): + 656 r'''Fits function to the data + 657 + 658 Parameters + 659 ---------- + 660 function : obj + 661 function to fit to the data. See fits.least_squares for details. + 662 fitrange : list + 663 Two element list containing the timeslices on which the fit is supposed to start and stop. + 664 Caution: This range is inclusive as opposed to standard python indexing. + 665 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. + 666 If not specified, self.prange or all timeslices are used. + 667 silent : bool + 668 Decides whether output is printed to the standard output. + 669 ''' + 670 if self.N != 1: + 671 raise Exception("Correlator must be projected before fitting") + 672 + 673 if fitrange is None: + 674 if self.prange: + 675 fitrange = self.prange + 676 else: + 677 fitrange = [0, self.T - 1] + 678 else: + 679 if not isinstance(fitrange, list): + 680 raise Exception("fitrange has to be a list with two elements") + 681 if len(fitrange) != 2: + 682 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") + 683 + 684 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] + 685 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] + 686 result = least_squares(xs, ys, function, silent=silent, **kwargs) + 687 return result + 688 + 689 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): + 690 """ Extract a plateau value from a Corr object + 691 + 692 Parameters + 693 ---------- + 694 plateau_range : list + 695 list with two entries, indicating the first and the last timeslice + 696 of the plateau region. + 697 method : str + 698 method to extract the plateau. + 699 'fit' fits a constant to the plateau region + 700 'avg', 'average' or 'mean' just average over the given timeslices. + 701 auto_gamma : bool + 702 apply gamma_method with default parameters to the Corr. Defaults to None + 703 """ + 704 if not plateau_range: + 705 if self.prange: + 706 plateau_range = self.prange + 707 else: + 708 raise Exception("no plateau range provided") + 709 if self.N != 1: + 710 raise Exception("Correlator must be projected before getting a plateau.") + 711 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): + 712 raise Exception("plateau is undefined at all timeslices in plateaurange.") + 713 if auto_gamma: + 714 self.gamma_method() + 715 if method == "fit": + 716 def const_func(a, t): + 717 return a[0] + 718 return self.fit(const_func, plateau_range)[0] + 719 elif method in ["avg", "average", "mean"]: + 720 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) + 721 return returnvalue + 722 + 723 else: + 724 raise Exception("Unsupported plateau method: " + method) + 725 + 726 def set_prange(self, prange): + 727 """Sets the attribute prange of the Corr object.""" + 728 if not len(prange) == 2: + 729 raise Exception("prange must be a list or array with two values") + 730 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): + 731 raise Exception("Start and end point must be integers") + 732 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): + 733 raise Exception("Start and end point must define a range in the interval 0,T") + 734 + 735 self.prange = prange + 736 return + 737 + 738 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None): + 739 """Plots the correlator using the tag of the correlator as label if available. + 740 + 741 Parameters + 742 ---------- + 743 x_range : list + 744 list of two values, determining the range of the x-axis e.g. [4, 8]. + 745 comp : Corr or list of Corr + 746 Correlator or list of correlators which are plotted for comparison. + 747 The tags of these correlators are used as labels if available. + 748 logscale : bool + 749 Sets y-axis to logscale. + 750 plateau : Obs + 751 Plateau value to be visualized in the figure. + 752 fit_res : Fit_result + 753 Fit_result object to be visualized. + 754 ylabel : str + 755 Label for the y-axis. + 756 save : str + 757 path to file in which the figure should be saved. + 758 auto_gamma : bool + 759 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. + 760 hide_sigma : float + 761 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. + 762 references : list + 763 List of floating point values that are displayed as horizontal lines for reference. + 764 title : string + 765 Optional title of the figure. + 766 """ + 767 if self.N != 1: + 768 raise Exception("Correlator must be projected before plotting") + 769 + 770 if auto_gamma: + 771 self.gamma_method() + 772 + 773 if x_range is None: + 774 x_range = [0, self.T - 1] + 775 + 776 fig = plt.figure() + 777 ax1 = fig.add_subplot(111) + 778 + 779 x, y, y_err = self.plottable() + 780 if hide_sigma: + 781 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 + 782 else: + 783 hide_from = None + 784 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) + 785 if logscale: + 786 ax1.set_yscale('log') + 787 else: + 788 if y_range is None: + 789 try: + 790 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)]) + 791 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)]) + 792 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) + 793 except Exception: + 794 pass + 795 else: + 796 ax1.set_ylim(y_range) + 797 if comp: + 798 if isinstance(comp, (Corr, list)): + 799 for corr in comp if isinstance(comp, list) else [comp]: + 800 if auto_gamma: + 801 corr.gamma_method() + 802 x, y, y_err = corr.plottable() + 803 if hide_sigma: + 804 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 + 805 else: + 806 hide_from = None + 807 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) + 808 else: + 809 raise Exception("'comp' must be a correlator or a list of correlators.") + 810 + 811 if plateau: + 812 if isinstance(plateau, Obs): + 813 if auto_gamma: + 814 plateau.gamma_method() + 815 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) + 816 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') + 817 else: + 818 raise Exception("'plateau' must be an Obs") + 819 + 820 if references: + 821 if isinstance(references, list): + 822 for ref in references: + 823 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') + 824 else: + 825 raise Exception("'references' must be a list of floating pint values.") 826 - 827 if fit_res: - 828 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) - 829 ax1.plot(x_samples, - 830 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), - 831 ls='-', marker=',', lw=2) - 832 - 833 ax1.set_xlabel(r'$x_0 / a$') - 834 if ylabel: - 835 ax1.set_ylabel(ylabel) - 836 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) - 837 - 838 handles, labels = ax1.get_legend_handles_labels() - 839 if labels: - 840 ax1.legend() + 827 if self.prange: + 828 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') + 829 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') + 830 + 831 if fit_res: + 832 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) + 833 ax1.plot(x_samples, + 834 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), + 835 ls='-', marker=',', lw=2) + 836 + 837 ax1.set_xlabel(r'$x_0 / a$') + 838 if ylabel: + 839 ax1.set_ylabel(ylabel) + 840 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 841 - 842 if title: - 843 plt.title(title) - 844 - 845 plt.draw() - 846 - 847 if save: - 848 if isinstance(save, str): - 849 fig.savefig(save, bbox_inches='tight') - 850 else: - 851 raise Exception("'save' has to be a string.") - 852 - 853 def spaghetti_plot(self, logscale=True): - 854 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. - 855 - 856 Parameters - 857 ---------- - 858 logscale : bool - 859 Determines whether the scale of the y-axis is logarithmic or standard. - 860 """ - 861 if self.N != 1: - 862 raise Exception("Correlator needs to be projected first.") - 863 - 864 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist])) - 865 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] - 866 - 867 for name in mc_names: - 868 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T - 869 - 870 fig = plt.figure() - 871 ax = fig.add_subplot(111) - 872 for dat in data: - 873 ax.plot(x0_vals, dat, ls='-', marker='') - 874 - 875 if logscale is True: - 876 ax.set_yscale('log') - 877 - 878 ax.set_xlabel(r'$x_0 / a$') - 879 plt.title(name) - 880 plt.draw() + 842 handles, labels = ax1.get_legend_handles_labels() + 843 if labels: + 844 ax1.legend() + 845 + 846 if title: + 847 plt.title(title) + 848 + 849 plt.draw() + 850 + 851 if save: + 852 if isinstance(save, str): + 853 fig.savefig(save, bbox_inches='tight') + 854 else: + 855 raise Exception("'save' has to be a string.") + 856 + 857 def spaghetti_plot(self, logscale=True): + 858 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. + 859 + 860 Parameters + 861 ---------- + 862 logscale : bool + 863 Determines whether the scale of the y-axis is logarithmic or standard. + 864 """ + 865 if self.N != 1: + 866 raise Exception("Correlator needs to be projected first.") + 867 + 868 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist])) + 869 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] + 870 + 871 for name in mc_names: + 872 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T + 873 + 874 fig = plt.figure() + 875 ax = fig.add_subplot(111) + 876 for dat in data: + 877 ax.plot(x0_vals, dat, ls='-', marker='') + 878 + 879 if logscale is True: + 880 ax.set_yscale('log') 881 - 882 def dump(self, filename, datatype="json.gz", **kwargs): - 883 """Dumps the Corr into a file of chosen type - 884 Parameters - 885 ---------- - 886 filename : str - 887 Name of the file to be saved. - 888 datatype : str - 889 Format of the exported file. Supported formats include - 890 "json.gz" and "pickle" - 891 path : str - 892 specifies a custom path for the file (default '.') - 893 """ - 894 if datatype == "json.gz": - 895 from .input.json import dump_to_json - 896 if 'path' in kwargs: - 897 file_name = kwargs.get('path') + '/' + filename - 898 else: - 899 file_name = filename - 900 dump_to_json(self, file_name) - 901 elif datatype == "pickle": - 902 dump_object(self, filename, **kwargs) - 903 else: - 904 raise Exception("Unknown datatype " + str(datatype)) - 905 - 906 def print(self, print_range=None): - 907 print(self.__repr__(print_range)) - 908 - 909 def __repr__(self, print_range=None): - 910 if print_range is None: - 911 print_range = [0, None] + 882 ax.set_xlabel(r'$x_0 / a$') + 883 plt.title(name) + 884 plt.draw() + 885 + 886 def dump(self, filename, datatype="json.gz", **kwargs): + 887 """Dumps the Corr into a file of chosen type + 888 Parameters + 889 ---------- + 890 filename : str + 891 Name of the file to be saved. + 892 datatype : str + 893 Format of the exported file. Supported formats include + 894 "json.gz" and "pickle" + 895 path : str + 896 specifies a custom path for the file (default '.') + 897 """ + 898 if datatype == "json.gz": + 899 from .input.json import dump_to_json + 900 if 'path' in kwargs: + 901 file_name = kwargs.get('path') + '/' + filename + 902 else: + 903 file_name = filename + 904 dump_to_json(self, file_name) + 905 elif datatype == "pickle": + 906 dump_object(self, filename, **kwargs) + 907 else: + 908 raise Exception("Unknown datatype " + str(datatype)) + 909 + 910 def print(self, print_range=None): + 911 print(self.__repr__(print_range)) 912 - 913 content_string = "" - 914 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 - 915 - 916 if self.tag is not None: - 917 content_string += "Description: " + self.tag + "\n" - 918 if self.N != 1: - 919 return content_string - 920 - 921 if print_range[1]: - 922 print_range[1] += 1 - 923 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' - 924 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): - 925 if sub_corr is None: - 926 content_string += str(i + print_range[0]) + '\n' - 927 else: - 928 content_string += str(i + print_range[0]) - 929 for element in sub_corr: - 930 content_string += '\t' + ' ' * int(element >= 0) + str(element) - 931 content_string += '\n' - 932 return content_string - 933 - 934 def __str__(self): - 935 return self.__repr__() - 936 - 937 # We define the basic operations, that can be performed with correlators. - 938 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. - 939 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. - 940 # One could try and tell Obs to check if the y in __mul__ is a Corr and - 941 - 942 def __add__(self, y): - 943 if isinstance(y, Corr): - 944 if ((self.N != y.N) or (self.T != y.T)): - 945 raise Exception("Addition of Corrs with different shape") - 946 newcontent = [] - 947 for t in range(self.T): - 948 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): - 949 newcontent.append(None) - 950 else: - 951 newcontent.append(self.content[t] + y.content[t]) - 952 return Corr(newcontent) - 953 - 954 elif isinstance(y, (Obs, int, float, CObs)): - 955 newcontent = [] - 956 for t in range(self.T): - 957 if _check_for_none(self, self.content[t]): - 958 newcontent.append(None) - 959 else: - 960 newcontent.append(self.content[t] + y) - 961 return Corr(newcontent, prange=self.prange) - 962 elif isinstance(y, np.ndarray): - 963 if y.shape == (self.T,): - 964 return Corr(list((np.array(self.content).T + y).T)) - 965 else: - 966 raise ValueError("operands could not be broadcast together") - 967 else: - 968 raise TypeError("Corr + wrong type") - 969 - 970 def __mul__(self, y): - 971 if isinstance(y, Corr): - 972 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): - 973 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") - 974 newcontent = [] - 975 for t in range(self.T): - 976 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): - 977 newcontent.append(None) - 978 else: - 979 newcontent.append(self.content[t] * y.content[t]) - 980 return Corr(newcontent) - 981 - 982 elif isinstance(y, (Obs, int, float, CObs)): - 983 newcontent = [] - 984 for t in range(self.T): - 985 if _check_for_none(self, self.content[t]): - 986 newcontent.append(None) - 987 else: - 988 newcontent.append(self.content[t] * y) - 989 return Corr(newcontent, prange=self.prange) - 990 elif isinstance(y, np.ndarray): - 991 if y.shape == (self.T,): - 992 return Corr(list((np.array(self.content).T * y).T)) - 993 else: - 994 raise ValueError("operands could not be broadcast together") - 995 else: - 996 raise TypeError("Corr * wrong type") - 997 - 998 def __truediv__(self, y): - 999 if isinstance(y, Corr): -1000 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1001 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1002 newcontent = [] -1003 for t in range(self.T): -1004 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1005 newcontent.append(None) -1006 else: -1007 newcontent.append(self.content[t] / y.content[t]) -1008 for t in range(self.T): -1009 if _check_for_none(self, newcontent[t]): -1010 continue -1011 if np.isnan(np.sum(newcontent[t]).value): -1012 newcontent[t] = None -1013 -1014 if all([item is None for item in newcontent]): -1015 raise Exception("Division returns completely undefined correlator") -1016 return Corr(newcontent) + 913 def __repr__(self, print_range=None): + 914 if print_range is None: + 915 print_range = [0, None] + 916 + 917 content_string = "" + 918 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 + 919 + 920 if self.tag is not None: + 921 content_string += "Description: " + self.tag + "\n" + 922 if self.N != 1: + 923 return content_string + 924 + 925 if print_range[1]: + 926 print_range[1] += 1 + 927 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' + 928 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): + 929 if sub_corr is None: + 930 content_string += str(i + print_range[0]) + '\n' + 931 else: + 932 content_string += str(i + print_range[0]) + 933 for element in sub_corr: + 934 content_string += '\t' + ' ' * int(element >= 0) + str(element) + 935 content_string += '\n' + 936 return content_string + 937 + 938 def __str__(self): + 939 return self.__repr__() + 940 + 941 # We define the basic operations, that can be performed with correlators. + 942 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. + 943 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. + 944 # One could try and tell Obs to check if the y in __mul__ is a Corr and + 945 + 946 def __add__(self, y): + 947 if isinstance(y, Corr): + 948 if ((self.N != y.N) or (self.T != y.T)): + 949 raise Exception("Addition of Corrs with different shape") + 950 newcontent = [] + 951 for t in range(self.T): + 952 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): + 953 newcontent.append(None) + 954 else: + 955 newcontent.append(self.content[t] + y.content[t]) + 956 return Corr(newcontent) + 957 + 958 elif isinstance(y, (Obs, int, float, CObs)): + 959 newcontent = [] + 960 for t in range(self.T): + 961 if _check_for_none(self, self.content[t]): + 962 newcontent.append(None) + 963 else: + 964 newcontent.append(self.content[t] + y) + 965 return Corr(newcontent, prange=self.prange) + 966 elif isinstance(y, np.ndarray): + 967 if y.shape == (self.T,): + 968 return Corr(list((np.array(self.content).T + y).T)) + 969 else: + 970 raise ValueError("operands could not be broadcast together") + 971 else: + 972 raise TypeError("Corr + wrong type") + 973 + 974 def __mul__(self, y): + 975 if isinstance(y, Corr): + 976 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): + 977 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") + 978 newcontent = [] + 979 for t in range(self.T): + 980 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): + 981 newcontent.append(None) + 982 else: + 983 newcontent.append(self.content[t] * y.content[t]) + 984 return Corr(newcontent) + 985 + 986 elif isinstance(y, (Obs, int, float, CObs)): + 987 newcontent = [] + 988 for t in range(self.T): + 989 if _check_for_none(self, self.content[t]): + 990 newcontent.append(None) + 991 else: + 992 newcontent.append(self.content[t] * y) + 993 return Corr(newcontent, prange=self.prange) + 994 elif isinstance(y, np.ndarray): + 995 if y.shape == (self.T,): + 996 return Corr(list((np.array(self.content).T * y).T)) + 997 else: + 998 raise ValueError("operands could not be broadcast together") + 999 else: +1000 raise TypeError("Corr * wrong type") +1001 +1002 def __truediv__(self, y): +1003 if isinstance(y, Corr): +1004 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1005 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1006 newcontent = [] +1007 for t in range(self.T): +1008 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1009 newcontent.append(None) +1010 else: +1011 newcontent.append(self.content[t] / y.content[t]) +1012 for t in range(self.T): +1013 if _check_for_none(self, newcontent[t]): +1014 continue +1015 if np.isnan(np.sum(newcontent[t]).value): +1016 newcontent[t] = None 1017 -1018 elif isinstance(y, (Obs, CObs)): -1019 if isinstance(y, Obs): -1020 if y.value == 0: -1021 raise Exception('Division by zero will return undefined correlator') -1022 if isinstance(y, CObs): -1023 if y.is_zero(): -1024 raise Exception('Division by zero will return undefined correlator') -1025 -1026 newcontent = [] -1027 for t in range(self.T): -1028 if _check_for_none(self, self.content[t]): -1029 newcontent.append(None) -1030 else: -1031 newcontent.append(self.content[t] / y) -1032 return Corr(newcontent, prange=self.prange) -1033 -1034 elif isinstance(y, (int, float)): -1035 if y == 0: -1036 raise Exception('Division by zero will return undefined correlator') -1037 newcontent = [] -1038 for t in range(self.T): -1039 if _check_for_none(self, self.content[t]): -1040 newcontent.append(None) -1041 else: -1042 newcontent.append(self.content[t] / y) -1043 return Corr(newcontent, prange=self.prange) -1044 elif isinstance(y, np.ndarray): -1045 if y.shape == (self.T,): -1046 return Corr(list((np.array(self.content).T / y).T)) -1047 else: -1048 raise ValueError("operands could not be broadcast together") -1049 else: -1050 raise TypeError('Corr / wrong type') -1051 -1052 def __neg__(self): -1053 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] -1054 return Corr(newcontent, prange=self.prange) +1018 if all([item is None for item in newcontent]): +1019 raise Exception("Division returns completely undefined correlator") +1020 return Corr(newcontent) +1021 +1022 elif isinstance(y, (Obs, CObs)): +1023 if isinstance(y, Obs): +1024 if y.value == 0: +1025 raise Exception('Division by zero will return undefined correlator') +1026 if isinstance(y, CObs): +1027 if y.is_zero(): +1028 raise Exception('Division by zero will return undefined correlator') +1029 +1030 newcontent = [] +1031 for t in range(self.T): +1032 if _check_for_none(self, self.content[t]): +1033 newcontent.append(None) +1034 else: +1035 newcontent.append(self.content[t] / y) +1036 return Corr(newcontent, prange=self.prange) +1037 +1038 elif isinstance(y, (int, float)): +1039 if y == 0: +1040 raise Exception('Division by zero will return undefined correlator') +1041 newcontent = [] +1042 for t in range(self.T): +1043 if _check_for_none(self, self.content[t]): +1044 newcontent.append(None) +1045 else: +1046 newcontent.append(self.content[t] / y) +1047 return Corr(newcontent, prange=self.prange) +1048 elif isinstance(y, np.ndarray): +1049 if y.shape == (self.T,): +1050 return Corr(list((np.array(self.content).T / y).T)) +1051 else: +1052 raise ValueError("operands could not be broadcast together") +1053 else: +1054 raise TypeError('Corr / wrong type') 1055 -1056 def __sub__(self, y): -1057 return self + (-y) -1058 -1059 def __pow__(self, y): -1060 if isinstance(y, (Obs, int, float, CObs)): -1061 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] -1062 return Corr(newcontent, prange=self.prange) -1063 else: -1064 raise TypeError('Type of exponent not supported') -1065 -1066 def __abs__(self): -1067 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] -1068 return Corr(newcontent, prange=self.prange) +1056 def __neg__(self): +1057 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] +1058 return Corr(newcontent, prange=self.prange) +1059 +1060 def __sub__(self, y): +1061 return self + (-y) +1062 +1063 def __pow__(self, y): +1064 if isinstance(y, (Obs, int, float, CObs)): +1065 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] +1066 return Corr(newcontent, prange=self.prange) +1067 else: +1068 raise TypeError('Type of exponent not supported') 1069 -1070 # The numpy functions: -1071 def sqrt(self): -1072 return self ** 0.5 +1070 def __abs__(self): +1071 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] +1072 return Corr(newcontent, prange=self.prange) 1073 -1074 def log(self): -1075 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] -1076 return Corr(newcontent, prange=self.prange) +1074 # The numpy functions: +1075 def sqrt(self): +1076 return self ** 0.5 1077 -1078 def exp(self): -1079 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] +1078 def log(self): +1079 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] 1080 return Corr(newcontent, prange=self.prange) 1081 -1082 def _apply_func_to_corr(self, func): -1083 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] -1084 for t in range(self.T): -1085 if _check_for_none(self, newcontent[t]): -1086 continue -1087 if np.isnan(np.sum(newcontent[t]).value): -1088 newcontent[t] = None -1089 if all([item is None for item in newcontent]): -1090 raise Exception('Operation returns undefined correlator') -1091 return Corr(newcontent) -1092 -1093 def sin(self): -1094 return self._apply_func_to_corr(np.sin) -1095 -1096 def cos(self): -1097 return self._apply_func_to_corr(np.cos) -1098 -1099 def tan(self): -1100 return self._apply_func_to_corr(np.tan) -1101 -1102 def sinh(self): -1103 return self._apply_func_to_corr(np.sinh) -1104 -1105 def cosh(self): -1106 return self._apply_func_to_corr(np.cosh) -1107 -1108 def tanh(self): -1109 return self._apply_func_to_corr(np.tanh) -1110 -1111 def arcsin(self): -1112 return self._apply_func_to_corr(np.arcsin) -1113 -1114 def arccos(self): -1115 return self._apply_func_to_corr(np.arccos) -1116 -1117 def arctan(self): -1118 return self._apply_func_to_corr(np.arctan) -1119 -1120 def arcsinh(self): -1121 return self._apply_func_to_corr(np.arcsinh) -1122 -1123 def arccosh(self): -1124 return self._apply_func_to_corr(np.arccosh) -1125 -1126 def arctanh(self): -1127 return self._apply_func_to_corr(np.arctanh) -1128 -1129 # Right hand side operations (require tweak in main module to work) -1130 def __radd__(self, y): -1131 return self + y +1082 def exp(self): +1083 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] +1084 return Corr(newcontent, prange=self.prange) +1085 +1086 def _apply_func_to_corr(self, func): +1087 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] +1088 for t in range(self.T): +1089 if _check_for_none(self, newcontent[t]): +1090 continue +1091 if np.isnan(np.sum(newcontent[t]).value): +1092 newcontent[t] = None +1093 if all([item is None for item in newcontent]): +1094 raise Exception('Operation returns undefined correlator') +1095 return Corr(newcontent) +1096 +1097 def sin(self): +1098 return self._apply_func_to_corr(np.sin) +1099 +1100 def cos(self): +1101 return self._apply_func_to_corr(np.cos) +1102 +1103 def tan(self): +1104 return self._apply_func_to_corr(np.tan) +1105 +1106 def sinh(self): +1107 return self._apply_func_to_corr(np.sinh) +1108 +1109 def cosh(self): +1110 return self._apply_func_to_corr(np.cosh) +1111 +1112 def tanh(self): +1113 return self._apply_func_to_corr(np.tanh) +1114 +1115 def arcsin(self): +1116 return self._apply_func_to_corr(np.arcsin) +1117 +1118 def arccos(self): +1119 return self._apply_func_to_corr(np.arccos) +1120 +1121 def arctan(self): +1122 return self._apply_func_to_corr(np.arctan) +1123 +1124 def arcsinh(self): +1125 return self._apply_func_to_corr(np.arcsinh) +1126 +1127 def arccosh(self): +1128 return self._apply_func_to_corr(np.arccosh) +1129 +1130 def arctanh(self): +1131 return self._apply_func_to_corr(np.arctanh) 1132 -1133 def __rsub__(self, y): -1134 return -self + y -1135 -1136 def __rmul__(self, y): -1137 return self * y -1138 -1139 def __rtruediv__(self, y): -1140 return (self / y) ** (-1) -1141 -1142 @property -1143 def real(self): -1144 def return_real(obs_OR_cobs): -1145 if isinstance(obs_OR_cobs, CObs): -1146 return obs_OR_cobs.real -1147 else: -1148 return obs_OR_cobs -1149 -1150 return self._apply_func_to_corr(return_real) -1151 -1152 @property -1153 def imag(self): -1154 def return_imag(obs_OR_cobs): -1155 if isinstance(obs_OR_cobs, CObs): -1156 return obs_OR_cobs.imag -1157 else: -1158 return obs_OR_cobs * 0 # So it stays the right type -1159 -1160 return self._apply_func_to_corr(return_imag) -1161 -1162 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): -1163 r''' Project large correlation matrix to lowest states -1164 -1165 This method can be used to reduce the size of an (N x N) correlation matrix -1166 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise -1167 is still small. +1133 # Right hand side operations (require tweak in main module to work) +1134 def __radd__(self, y): +1135 return self + y +1136 +1137 def __rsub__(self, y): +1138 return -self + y +1139 +1140 def __rmul__(self, y): +1141 return self * y +1142 +1143 def __rtruediv__(self, y): +1144 return (self / y) ** (-1) +1145 +1146 @property +1147 def real(self): +1148 def return_real(obs_OR_cobs): +1149 if isinstance(obs_OR_cobs, CObs): +1150 return obs_OR_cobs.real +1151 else: +1152 return obs_OR_cobs +1153 +1154 return self._apply_func_to_corr(return_real) +1155 +1156 @property +1157 def imag(self): +1158 def return_imag(obs_OR_cobs): +1159 if isinstance(obs_OR_cobs, CObs): +1160 return obs_OR_cobs.imag +1161 else: +1162 return obs_OR_cobs * 0 # So it stays the right type +1163 +1164 return self._apply_func_to_corr(return_imag) +1165 +1166 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): +1167 r''' Project large correlation matrix to lowest states 1168 -1169 Parameters -1170 ---------- -1171 Ntrunc: int -1172 Rank of the target matrix. -1173 tproj: int -1174 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. -1175 The default value is 3. -1176 t0proj: int -1177 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly -1178 discouraged for O(a) improved theories, since the correctness of the procedure -1179 cannot be granted in this case. The default value is 2. -1180 basematrix : Corr -1181 Correlation matrix that is used to determine the eigenvectors of the -1182 lowest states based on a GEVP. basematrix is taken to be the Corr itself if -1183 is is not specified. -1184 -1185 Notes -1186 ----- -1187 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving -1188 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}$ -1189 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the -1190 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via -1191 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large -1192 correlation matrix and to remove some noise that is added by irrelevant operators. -1193 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated -1194 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. -1195 ''' -1196 -1197 if self.N == 1: -1198 raise Exception('Method cannot be applied to one-dimensional correlators.') -1199 if basematrix is None: -1200 basematrix = self -1201 if Ntrunc >= basematrix.N: -1202 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) -1203 if basematrix.N != self.N: -1204 raise Exception('basematrix and targetmatrix have to be of the same size.') -1205 -1206 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] -1207 -1208 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) -1209 rmat = [] -1210 for t in range(basematrix.T): -1211 for i in range(Ntrunc): -1212 for j in range(Ntrunc): -1213 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] -1214 rmat.append(np.copy(tmpmat)) -1215 -1216 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] -1217 return Corr(newcontent) -1218 +1169 This method can be used to reduce the size of an (N x N) correlation matrix +1170 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise +1171 is still small. +1172 +1173 Parameters +1174 ---------- +1175 Ntrunc: int +1176 Rank of the target matrix. +1177 tproj: int +1178 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. +1179 The default value is 3. +1180 t0proj: int +1181 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly +1182 discouraged for O(a) improved theories, since the correctness of the procedure +1183 cannot be granted in this case. The default value is 2. +1184 basematrix : Corr +1185 Correlation matrix that is used to determine the eigenvectors of the +1186 lowest states based on a GEVP. basematrix is taken to be the Corr itself if +1187 is is not specified. +1188 +1189 Notes +1190 ----- +1191 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving +1192 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}$ +1193 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the +1194 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via +1195 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large +1196 correlation matrix and to remove some noise that is added by irrelevant operators. +1197 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated +1198 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. +1199 ''' +1200 +1201 if self.N == 1: +1202 raise Exception('Method cannot be applied to one-dimensional correlators.') +1203 if basematrix is None: +1204 basematrix = self +1205 if Ntrunc >= basematrix.N: +1206 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) +1207 if basematrix.N != self.N: +1208 raise Exception('basematrix and targetmatrix have to be of the same size.') +1209 +1210 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] +1211 +1212 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) +1213 rmat = [] +1214 for t in range(basematrix.T): +1215 for i in range(Ntrunc): +1216 for j in range(Ntrunc): +1217 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] +1218 rmat.append(np.copy(tmpmat)) 1219 -1220def _sort_vectors(vec_set, ts): -1221 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" -1222 reference_sorting = np.array(vec_set[ts]) -1223 N = reference_sorting.shape[0] -1224 sorted_vec_set = [] -1225 for t in range(len(vec_set)): -1226 if vec_set[t] is None: -1227 sorted_vec_set.append(None) -1228 elif not t == ts: -1229 perms = [list(o) for o in permutations([i for i in range(N)], N)] -1230 best_score = 0 -1231 for perm in perms: -1232 current_score = 1 -1233 for k in range(N): -1234 new_sorting = reference_sorting.copy() -1235 new_sorting[perm[k], :] = vec_set[t][k] -1236 current_score *= abs(np.linalg.det(new_sorting)) -1237 if current_score > best_score: -1238 best_score = current_score -1239 best_perm = perm -1240 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) -1241 else: -1242 sorted_vec_set.append(vec_set[t]) -1243 -1244 return sorted_vec_set -1245 -1246 -1247def _check_for_none(corr, entry): -1248 """Checks if entry for correlator corr is None""" -1249 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2 +1220 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] +1221 return Corr(newcontent) +1222 +1223 +1224def _sort_vectors(vec_set, ts): +1225 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" +1226 reference_sorting = np.array(vec_set[ts]) +1227 N = reference_sorting.shape[0] +1228 sorted_vec_set = [] +1229 for t in range(len(vec_set)): +1230 if vec_set[t] is None: +1231 sorted_vec_set.append(None) +1232 elif not t == ts: +1233 perms = [list(o) for o in permutations([i for i in range(N)], N)] +1234 best_score = 0 +1235 for perm in perms: +1236 current_score = 1 +1237 for k in range(N): +1238 new_sorting = reference_sorting.copy() +1239 new_sorting[perm[k], :] = vec_set[t][k] +1240 current_score *= abs(np.linalg.det(new_sorting)) +1241 if current_score > best_score: +1242 best_score = current_score +1243 best_perm = perm +1244 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) +1245 else: +1246 sorted_vec_set.append(vec_set[t]) +1247 +1248 return sorted_vec_set +1249 1250 -1251 -1252def _GEVP_solver(Gt, G0): -1253 """Helper function for solving the GEVP and sorting the eigenvectors. +1251def _check_for_none(corr, entry): +1252 """Checks if entry for correlator corr is None""" +1253 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2 1254 -1255 The helper function assumes that both provided matrices are symmetric and -1256 only processes the lower triangular part of both matrices. In case the matrices -1257 are not symmetric the upper triangular parts are effectively discarded.""" -1258 return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1] +1255 +1256def _GEVP_solver(Gt, G0): +1257 """Helper function for solving the GEVP and sorting the eigenvectors. +1258 +1259 The helper function assumes that both provided matrices are symmetric and +1260 only processes the lower triangular part of both matrices. In case the matrices +1261 are not symmetric the upper triangular parts are effectively discarded.""" +1262 return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1] @@ -2079,617 +2083,621 @@ 605 for t in range(self.T - 1): 606 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): 607 newcontent.append(None) - 608 else: - 609 newcontent.append(self.content[t] / self.content[t + 1]) - 610 if(all([x is None for x in newcontent])): - 611 raise Exception('m_eff is undefined at all timeslices') - 612 - 613 return np.log(Corr(newcontent, padding=[0, 1])) + 608 elif self.content[t][0].value / self.content[t + 1][0].value < 0: + 609 newcontent.append(None) + 610 else: + 611 newcontent.append(self.content[t] / self.content[t + 1]) + 612 if(all([x is None for x in newcontent])): + 613 raise Exception('m_eff is undefined at all timeslices') 614 - 615 elif variant in ['periodic', 'cosh', 'sinh']: - 616 if variant in ['periodic', 'cosh']: - 617 func = anp.cosh - 618 else: - 619 func = anp.sinh - 620 - 621 def root_function(x, d): - 622 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d - 623 - 624 newcontent = [] - 625 for t in range(self.T - 1): - 626 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): - 627 newcontent.append(None) - 628 # Fill the two timeslices in the middle of the lattice with their predecessors - 629 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: - 630 newcontent.append(newcontent[-1]) - 631 else: - 632 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) - 633 if(all([x is None for x in newcontent])): - 634 raise Exception('m_eff is undefined at all timeslices') - 635 - 636 return Corr(newcontent, padding=[0, 1]) - 637 - 638 elif variant == 'arccosh': - 639 newcontent = [] - 640 for t in range(1, self.T - 1): - 641 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): - 642 newcontent.append(None) - 643 else: - 644 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) - 645 if(all([x is None for x in newcontent])): - 646 raise Exception("m_eff is undefined at all timeslices") - 647 return np.arccosh(Corr(newcontent, padding=[1, 1])) - 648 - 649 else: - 650 raise Exception('Unknown variant.') - 651 - 652 def fit(self, function, fitrange=None, silent=False, **kwargs): - 653 r'''Fits function to the data - 654 - 655 Parameters - 656 ---------- - 657 function : obj - 658 function to fit to the data. See fits.least_squares for details. - 659 fitrange : list - 660 Two element list containing the timeslices on which the fit is supposed to start and stop. - 661 Caution: This range is inclusive as opposed to standard python indexing. - 662 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. - 663 If not specified, self.prange or all timeslices are used. - 664 silent : bool - 665 Decides whether output is printed to the standard output. - 666 ''' - 667 if self.N != 1: - 668 raise Exception("Correlator must be projected before fitting") - 669 - 670 if fitrange is None: - 671 if self.prange: - 672 fitrange = self.prange - 673 else: - 674 fitrange = [0, self.T - 1] - 675 else: - 676 if not isinstance(fitrange, list): - 677 raise Exception("fitrange has to be a list with two elements") - 678 if len(fitrange) != 2: - 679 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") - 680 - 681 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] - 682 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] - 683 result = least_squares(xs, ys, function, silent=silent, **kwargs) - 684 return result - 685 - 686 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): - 687 """ Extract a plateau value from a Corr object - 688 - 689 Parameters - 690 ---------- - 691 plateau_range : list - 692 list with two entries, indicating the first and the last timeslice - 693 of the plateau region. - 694 method : str - 695 method to extract the plateau. - 696 'fit' fits a constant to the plateau region - 697 'avg', 'average' or 'mean' just average over the given timeslices. - 698 auto_gamma : bool - 699 apply gamma_method with default parameters to the Corr. Defaults to None - 700 """ - 701 if not plateau_range: - 702 if self.prange: - 703 plateau_range = self.prange - 704 else: - 705 raise Exception("no plateau range provided") - 706 if self.N != 1: - 707 raise Exception("Correlator must be projected before getting a plateau.") - 708 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): - 709 raise Exception("plateau is undefined at all timeslices in plateaurange.") - 710 if auto_gamma: - 711 self.gamma_method() - 712 if method == "fit": - 713 def const_func(a, t): - 714 return a[0] - 715 return self.fit(const_func, plateau_range)[0] - 716 elif method in ["avg", "average", "mean"]: - 717 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) - 718 return returnvalue - 719 - 720 else: - 721 raise Exception("Unsupported plateau method: " + method) - 722 - 723 def set_prange(self, prange): - 724 """Sets the attribute prange of the Corr object.""" - 725 if not len(prange) == 2: - 726 raise Exception("prange must be a list or array with two values") - 727 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): - 728 raise Exception("Start and end point must be integers") - 729 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): - 730 raise Exception("Start and end point must define a range in the interval 0,T") - 731 - 732 self.prange = prange - 733 return - 734 - 735 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None): - 736 """Plots the correlator using the tag of the correlator as label if available. - 737 - 738 Parameters - 739 ---------- - 740 x_range : list - 741 list of two values, determining the range of the x-axis e.g. [4, 8]. - 742 comp : Corr or list of Corr - 743 Correlator or list of correlators which are plotted for comparison. - 744 The tags of these correlators are used as labels if available. - 745 logscale : bool - 746 Sets y-axis to logscale. - 747 plateau : Obs - 748 Plateau value to be visualized in the figure. - 749 fit_res : Fit_result - 750 Fit_result object to be visualized. - 751 ylabel : str - 752 Label for the y-axis. - 753 save : str - 754 path to file in which the figure should be saved. - 755 auto_gamma : bool - 756 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. - 757 hide_sigma : float - 758 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. - 759 references : list - 760 List of floating point values that are displayed as horizontal lines for reference. - 761 title : string - 762 Optional title of the figure. - 763 """ - 764 if self.N != 1: - 765 raise Exception("Correlator must be projected before plotting") - 766 - 767 if auto_gamma: - 768 self.gamma_method() - 769 - 770 if x_range is None: - 771 x_range = [0, self.T - 1] - 772 - 773 fig = plt.figure() - 774 ax1 = fig.add_subplot(111) - 775 - 776 x, y, y_err = self.plottable() - 777 if hide_sigma: - 778 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 - 779 else: - 780 hide_from = None - 781 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) - 782 if logscale: - 783 ax1.set_yscale('log') - 784 else: - 785 if y_range is None: - 786 try: - 787 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)]) - 788 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)]) - 789 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) - 790 except Exception: - 791 pass - 792 else: - 793 ax1.set_ylim(y_range) - 794 if comp: - 795 if isinstance(comp, (Corr, list)): - 796 for corr in comp if isinstance(comp, list) else [comp]: - 797 if auto_gamma: - 798 corr.gamma_method() - 799 x, y, y_err = corr.plottable() - 800 if hide_sigma: - 801 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 - 802 else: - 803 hide_from = None - 804 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) - 805 else: - 806 raise Exception("'comp' must be a correlator or a list of correlators.") - 807 - 808 if plateau: - 809 if isinstance(plateau, Obs): - 810 if auto_gamma: - 811 plateau.gamma_method() - 812 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) - 813 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') - 814 else: - 815 raise Exception("'plateau' must be an Obs") - 816 - 817 if references: - 818 if isinstance(references, list): - 819 for ref in references: - 820 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') - 821 else: - 822 raise Exception("'references' must be a list of floating pint values.") - 823 - 824 if self.prange: - 825 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') - 826 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') + 615 return np.log(Corr(newcontent, padding=[0, 1])) + 616 + 617 elif variant in ['periodic', 'cosh', 'sinh']: + 618 if variant in ['periodic', 'cosh']: + 619 func = anp.cosh + 620 else: + 621 func = anp.sinh + 622 + 623 def root_function(x, d): + 624 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d + 625 + 626 newcontent = [] + 627 for t in range(self.T - 1): + 628 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): + 629 newcontent.append(None) + 630 # Fill the two timeslices in the middle of the lattice with their predecessors + 631 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: + 632 newcontent.append(newcontent[-1]) + 633 elif self.content[t][0].value / self.content[t + 1][0].value < 0: + 634 newcontent.append(None) + 635 else: + 636 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) + 637 if(all([x is None for x in newcontent])): + 638 raise Exception('m_eff is undefined at all timeslices') + 639 + 640 return Corr(newcontent, padding=[0, 1]) + 641 + 642 elif variant == 'arccosh': + 643 newcontent = [] + 644 for t in range(1, self.T - 1): + 645 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): + 646 newcontent.append(None) + 647 else: + 648 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) + 649 if(all([x is None for x in newcontent])): + 650 raise Exception("m_eff is undefined at all timeslices") + 651 return np.arccosh(Corr(newcontent, padding=[1, 1])) + 652 + 653 else: + 654 raise Exception('Unknown variant.') + 655 + 656 def fit(self, function, fitrange=None, silent=False, **kwargs): + 657 r'''Fits function to the data + 658 + 659 Parameters + 660 ---------- + 661 function : obj + 662 function to fit to the data. See fits.least_squares for details. + 663 fitrange : list + 664 Two element list containing the timeslices on which the fit is supposed to start and stop. + 665 Caution: This range is inclusive as opposed to standard python indexing. + 666 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. + 667 If not specified, self.prange or all timeslices are used. + 668 silent : bool + 669 Decides whether output is printed to the standard output. + 670 ''' + 671 if self.N != 1: + 672 raise Exception("Correlator must be projected before fitting") + 673 + 674 if fitrange is None: + 675 if self.prange: + 676 fitrange = self.prange + 677 else: + 678 fitrange = [0, self.T - 1] + 679 else: + 680 if not isinstance(fitrange, list): + 681 raise Exception("fitrange has to be a list with two elements") + 682 if len(fitrange) != 2: + 683 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") + 684 + 685 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] + 686 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] + 687 result = least_squares(xs, ys, function, silent=silent, **kwargs) + 688 return result + 689 + 690 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): + 691 """ Extract a plateau value from a Corr object + 692 + 693 Parameters + 694 ---------- + 695 plateau_range : list + 696 list with two entries, indicating the first and the last timeslice + 697 of the plateau region. + 698 method : str + 699 method to extract the plateau. + 700 'fit' fits a constant to the plateau region + 701 'avg', 'average' or 'mean' just average over the given timeslices. + 702 auto_gamma : bool + 703 apply gamma_method with default parameters to the Corr. Defaults to None + 704 """ + 705 if not plateau_range: + 706 if self.prange: + 707 plateau_range = self.prange + 708 else: + 709 raise Exception("no plateau range provided") + 710 if self.N != 1: + 711 raise Exception("Correlator must be projected before getting a plateau.") + 712 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): + 713 raise Exception("plateau is undefined at all timeslices in plateaurange.") + 714 if auto_gamma: + 715 self.gamma_method() + 716 if method == "fit": + 717 def const_func(a, t): + 718 return a[0] + 719 return self.fit(const_func, plateau_range)[0] + 720 elif method in ["avg", "average", "mean"]: + 721 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) + 722 return returnvalue + 723 + 724 else: + 725 raise Exception("Unsupported plateau method: " + method) + 726 + 727 def set_prange(self, prange): + 728 """Sets the attribute prange of the Corr object.""" + 729 if not len(prange) == 2: + 730 raise Exception("prange must be a list or array with two values") + 731 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): + 732 raise Exception("Start and end point must be integers") + 733 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): + 734 raise Exception("Start and end point must define a range in the interval 0,T") + 735 + 736 self.prange = prange + 737 return + 738 + 739 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None): + 740 """Plots the correlator using the tag of the correlator as label if available. + 741 + 742 Parameters + 743 ---------- + 744 x_range : list + 745 list of two values, determining the range of the x-axis e.g. [4, 8]. + 746 comp : Corr or list of Corr + 747 Correlator or list of correlators which are plotted for comparison. + 748 The tags of these correlators are used as labels if available. + 749 logscale : bool + 750 Sets y-axis to logscale. + 751 plateau : Obs + 752 Plateau value to be visualized in the figure. + 753 fit_res : Fit_result + 754 Fit_result object to be visualized. + 755 ylabel : str + 756 Label for the y-axis. + 757 save : str + 758 path to file in which the figure should be saved. + 759 auto_gamma : bool + 760 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. + 761 hide_sigma : float + 762 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. + 763 references : list + 764 List of floating point values that are displayed as horizontal lines for reference. + 765 title : string + 766 Optional title of the figure. + 767 """ + 768 if self.N != 1: + 769 raise Exception("Correlator must be projected before plotting") + 770 + 771 if auto_gamma: + 772 self.gamma_method() + 773 + 774 if x_range is None: + 775 x_range = [0, self.T - 1] + 776 + 777 fig = plt.figure() + 778 ax1 = fig.add_subplot(111) + 779 + 780 x, y, y_err = self.plottable() + 781 if hide_sigma: + 782 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 + 783 else: + 784 hide_from = None + 785 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) + 786 if logscale: + 787 ax1.set_yscale('log') + 788 else: + 789 if y_range is None: + 790 try: + 791 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)]) + 792 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)]) + 793 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) + 794 except Exception: + 795 pass + 796 else: + 797 ax1.set_ylim(y_range) + 798 if comp: + 799 if isinstance(comp, (Corr, list)): + 800 for corr in comp if isinstance(comp, list) else [comp]: + 801 if auto_gamma: + 802 corr.gamma_method() + 803 x, y, y_err = corr.plottable() + 804 if hide_sigma: + 805 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 + 806 else: + 807 hide_from = None + 808 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) + 809 else: + 810 raise Exception("'comp' must be a correlator or a list of correlators.") + 811 + 812 if plateau: + 813 if isinstance(plateau, Obs): + 814 if auto_gamma: + 815 plateau.gamma_method() + 816 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) + 817 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') + 818 else: + 819 raise Exception("'plateau' must be an Obs") + 820 + 821 if references: + 822 if isinstance(references, list): + 823 for ref in references: + 824 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') + 825 else: + 826 raise Exception("'references' must be a list of floating pint values.") 827 - 828 if fit_res: - 829 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) - 830 ax1.plot(x_samples, - 831 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), - 832 ls='-', marker=',', lw=2) - 833 - 834 ax1.set_xlabel(r'$x_0 / a$') - 835 if ylabel: - 836 ax1.set_ylabel(ylabel) - 837 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) - 838 - 839 handles, labels = ax1.get_legend_handles_labels() - 840 if labels: - 841 ax1.legend() + 828 if self.prange: + 829 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') + 830 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') + 831 + 832 if fit_res: + 833 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) + 834 ax1.plot(x_samples, + 835 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), + 836 ls='-', marker=',', lw=2) + 837 + 838 ax1.set_xlabel(r'$x_0 / a$') + 839 if ylabel: + 840 ax1.set_ylabel(ylabel) + 841 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 842 - 843 if title: - 844 plt.title(title) - 845 - 846 plt.draw() - 847 - 848 if save: - 849 if isinstance(save, str): - 850 fig.savefig(save, bbox_inches='tight') - 851 else: - 852 raise Exception("'save' has to be a string.") - 853 - 854 def spaghetti_plot(self, logscale=True): - 855 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. - 856 - 857 Parameters - 858 ---------- - 859 logscale : bool - 860 Determines whether the scale of the y-axis is logarithmic or standard. - 861 """ - 862 if self.N != 1: - 863 raise Exception("Correlator needs to be projected first.") - 864 - 865 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist])) - 866 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] - 867 - 868 for name in mc_names: - 869 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T - 870 - 871 fig = plt.figure() - 872 ax = fig.add_subplot(111) - 873 for dat in data: - 874 ax.plot(x0_vals, dat, ls='-', marker='') - 875 - 876 if logscale is True: - 877 ax.set_yscale('log') - 878 - 879 ax.set_xlabel(r'$x_0 / a$') - 880 plt.title(name) - 881 plt.draw() + 843 handles, labels = ax1.get_legend_handles_labels() + 844 if labels: + 845 ax1.legend() + 846 + 847 if title: + 848 plt.title(title) + 849 + 850 plt.draw() + 851 + 852 if save: + 853 if isinstance(save, str): + 854 fig.savefig(save, bbox_inches='tight') + 855 else: + 856 raise Exception("'save' has to be a string.") + 857 + 858 def spaghetti_plot(self, logscale=True): + 859 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. + 860 + 861 Parameters + 862 ---------- + 863 logscale : bool + 864 Determines whether the scale of the y-axis is logarithmic or standard. + 865 """ + 866 if self.N != 1: + 867 raise Exception("Correlator needs to be projected first.") + 868 + 869 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist])) + 870 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] + 871 + 872 for name in mc_names: + 873 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T + 874 + 875 fig = plt.figure() + 876 ax = fig.add_subplot(111) + 877 for dat in data: + 878 ax.plot(x0_vals, dat, ls='-', marker='') + 879 + 880 if logscale is True: + 881 ax.set_yscale('log') 882 - 883 def dump(self, filename, datatype="json.gz", **kwargs): - 884 """Dumps the Corr into a file of chosen type - 885 Parameters - 886 ---------- - 887 filename : str - 888 Name of the file to be saved. - 889 datatype : str - 890 Format of the exported file. Supported formats include - 891 "json.gz" and "pickle" - 892 path : str - 893 specifies a custom path for the file (default '.') - 894 """ - 895 if datatype == "json.gz": - 896 from .input.json import dump_to_json - 897 if 'path' in kwargs: - 898 file_name = kwargs.get('path') + '/' + filename - 899 else: - 900 file_name = filename - 901 dump_to_json(self, file_name) - 902 elif datatype == "pickle": - 903 dump_object(self, filename, **kwargs) - 904 else: - 905 raise Exception("Unknown datatype " + str(datatype)) - 906 - 907 def print(self, print_range=None): - 908 print(self.__repr__(print_range)) - 909 - 910 def __repr__(self, print_range=None): - 911 if print_range is None: - 912 print_range = [0, None] + 883 ax.set_xlabel(r'$x_0 / a$') + 884 plt.title(name) + 885 plt.draw() + 886 + 887 def dump(self, filename, datatype="json.gz", **kwargs): + 888 """Dumps the Corr into a file of chosen type + 889 Parameters + 890 ---------- + 891 filename : str + 892 Name of the file to be saved. + 893 datatype : str + 894 Format of the exported file. Supported formats include + 895 "json.gz" and "pickle" + 896 path : str + 897 specifies a custom path for the file (default '.') + 898 """ + 899 if datatype == "json.gz": + 900 from .input.json import dump_to_json + 901 if 'path' in kwargs: + 902 file_name = kwargs.get('path') + '/' + filename + 903 else: + 904 file_name = filename + 905 dump_to_json(self, file_name) + 906 elif datatype == "pickle": + 907 dump_object(self, filename, **kwargs) + 908 else: + 909 raise Exception("Unknown datatype " + str(datatype)) + 910 + 911 def print(self, print_range=None): + 912 print(self.__repr__(print_range)) 913 - 914 content_string = "" - 915 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 - 916 - 917 if self.tag is not None: - 918 content_string += "Description: " + self.tag + "\n" - 919 if self.N != 1: - 920 return content_string - 921 - 922 if print_range[1]: - 923 print_range[1] += 1 - 924 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' - 925 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): - 926 if sub_corr is None: - 927 content_string += str(i + print_range[0]) + '\n' - 928 else: - 929 content_string += str(i + print_range[0]) - 930 for element in sub_corr: - 931 content_string += '\t' + ' ' * int(element >= 0) + str(element) - 932 content_string += '\n' - 933 return content_string - 934 - 935 def __str__(self): - 936 return self.__repr__() - 937 - 938 # We define the basic operations, that can be performed with correlators. - 939 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. - 940 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. - 941 # One could try and tell Obs to check if the y in __mul__ is a Corr and - 942 - 943 def __add__(self, y): - 944 if isinstance(y, Corr): - 945 if ((self.N != y.N) or (self.T != y.T)): - 946 raise Exception("Addition of Corrs with different shape") - 947 newcontent = [] - 948 for t in range(self.T): - 949 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): - 950 newcontent.append(None) - 951 else: - 952 newcontent.append(self.content[t] + y.content[t]) - 953 return Corr(newcontent) - 954 - 955 elif isinstance(y, (Obs, int, float, CObs)): - 956 newcontent = [] - 957 for t in range(self.T): - 958 if _check_for_none(self, self.content[t]): - 959 newcontent.append(None) - 960 else: - 961 newcontent.append(self.content[t] + y) - 962 return Corr(newcontent, prange=self.prange) - 963 elif isinstance(y, np.ndarray): - 964 if y.shape == (self.T,): - 965 return Corr(list((np.array(self.content).T + y).T)) - 966 else: - 967 raise ValueError("operands could not be broadcast together") - 968 else: - 969 raise TypeError("Corr + wrong type") - 970 - 971 def __mul__(self, y): - 972 if isinstance(y, Corr): - 973 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): - 974 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") - 975 newcontent = [] - 976 for t in range(self.T): - 977 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): - 978 newcontent.append(None) - 979 else: - 980 newcontent.append(self.content[t] * y.content[t]) - 981 return Corr(newcontent) - 982 - 983 elif isinstance(y, (Obs, int, float, CObs)): - 984 newcontent = [] - 985 for t in range(self.T): - 986 if _check_for_none(self, self.content[t]): - 987 newcontent.append(None) - 988 else: - 989 newcontent.append(self.content[t] * y) - 990 return Corr(newcontent, prange=self.prange) - 991 elif isinstance(y, np.ndarray): - 992 if y.shape == (self.T,): - 993 return Corr(list((np.array(self.content).T * y).T)) - 994 else: - 995 raise ValueError("operands could not be broadcast together") - 996 else: - 997 raise TypeError("Corr * wrong type") - 998 - 999 def __truediv__(self, y): -1000 if isinstance(y, Corr): -1001 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1002 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1003 newcontent = [] -1004 for t in range(self.T): -1005 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1006 newcontent.append(None) -1007 else: -1008 newcontent.append(self.content[t] / y.content[t]) -1009 for t in range(self.T): -1010 if _check_for_none(self, newcontent[t]): -1011 continue -1012 if np.isnan(np.sum(newcontent[t]).value): -1013 newcontent[t] = None -1014 -1015 if all([item is None for item in newcontent]): -1016 raise Exception("Division returns completely undefined correlator") -1017 return Corr(newcontent) + 914 def __repr__(self, print_range=None): + 915 if print_range is None: + 916 print_range = [0, None] + 917 + 918 content_string = "" + 919 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 + 920 + 921 if self.tag is not None: + 922 content_string += "Description: " + self.tag + "\n" + 923 if self.N != 1: + 924 return content_string + 925 + 926 if print_range[1]: + 927 print_range[1] += 1 + 928 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' + 929 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): + 930 if sub_corr is None: + 931 content_string += str(i + print_range[0]) + '\n' + 932 else: + 933 content_string += str(i + print_range[0]) + 934 for element in sub_corr: + 935 content_string += '\t' + ' ' * int(element >= 0) + str(element) + 936 content_string += '\n' + 937 return content_string + 938 + 939 def __str__(self): + 940 return self.__repr__() + 941 + 942 # We define the basic operations, that can be performed with correlators. + 943 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. + 944 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. + 945 # One could try and tell Obs to check if the y in __mul__ is a Corr and + 946 + 947 def __add__(self, y): + 948 if isinstance(y, Corr): + 949 if ((self.N != y.N) or (self.T != y.T)): + 950 raise Exception("Addition of Corrs with different shape") + 951 newcontent = [] + 952 for t in range(self.T): + 953 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): + 954 newcontent.append(None) + 955 else: + 956 newcontent.append(self.content[t] + y.content[t]) + 957 return Corr(newcontent) + 958 + 959 elif isinstance(y, (Obs, int, float, CObs)): + 960 newcontent = [] + 961 for t in range(self.T): + 962 if _check_for_none(self, self.content[t]): + 963 newcontent.append(None) + 964 else: + 965 newcontent.append(self.content[t] + y) + 966 return Corr(newcontent, prange=self.prange) + 967 elif isinstance(y, np.ndarray): + 968 if y.shape == (self.T,): + 969 return Corr(list((np.array(self.content).T + y).T)) + 970 else: + 971 raise ValueError("operands could not be broadcast together") + 972 else: + 973 raise TypeError("Corr + wrong type") + 974 + 975 def __mul__(self, y): + 976 if isinstance(y, Corr): + 977 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): + 978 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") + 979 newcontent = [] + 980 for t in range(self.T): + 981 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): + 982 newcontent.append(None) + 983 else: + 984 newcontent.append(self.content[t] * y.content[t]) + 985 return Corr(newcontent) + 986 + 987 elif isinstance(y, (Obs, int, float, CObs)): + 988 newcontent = [] + 989 for t in range(self.T): + 990 if _check_for_none(self, self.content[t]): + 991 newcontent.append(None) + 992 else: + 993 newcontent.append(self.content[t] * y) + 994 return Corr(newcontent, prange=self.prange) + 995 elif isinstance(y, np.ndarray): + 996 if y.shape == (self.T,): + 997 return Corr(list((np.array(self.content).T * y).T)) + 998 else: + 999 raise ValueError("operands could not be broadcast together") +1000 else: +1001 raise TypeError("Corr * wrong type") +1002 +1003 def __truediv__(self, y): +1004 if isinstance(y, Corr): +1005 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1006 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1007 newcontent = [] +1008 for t in range(self.T): +1009 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1010 newcontent.append(None) +1011 else: +1012 newcontent.append(self.content[t] / y.content[t]) +1013 for t in range(self.T): +1014 if _check_for_none(self, newcontent[t]): +1015 continue +1016 if np.isnan(np.sum(newcontent[t]).value): +1017 newcontent[t] = None 1018 -1019 elif isinstance(y, (Obs, CObs)): -1020 if isinstance(y, Obs): -1021 if y.value == 0: -1022 raise Exception('Division by zero will return undefined correlator') -1023 if isinstance(y, CObs): -1024 if y.is_zero(): -1025 raise Exception('Division by zero will return undefined correlator') -1026 -1027 newcontent = [] -1028 for t in range(self.T): -1029 if _check_for_none(self, self.content[t]): -1030 newcontent.append(None) -1031 else: -1032 newcontent.append(self.content[t] / y) -1033 return Corr(newcontent, prange=self.prange) -1034 -1035 elif isinstance(y, (int, float)): -1036 if y == 0: -1037 raise Exception('Division by zero will return undefined correlator') -1038 newcontent = [] -1039 for t in range(self.T): -1040 if _check_for_none(self, self.content[t]): -1041 newcontent.append(None) -1042 else: -1043 newcontent.append(self.content[t] / y) -1044 return Corr(newcontent, prange=self.prange) -1045 elif isinstance(y, np.ndarray): -1046 if y.shape == (self.T,): -1047 return Corr(list((np.array(self.content).T / y).T)) -1048 else: -1049 raise ValueError("operands could not be broadcast together") -1050 else: -1051 raise TypeError('Corr / wrong type') -1052 -1053 def __neg__(self): -1054 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] -1055 return Corr(newcontent, prange=self.prange) +1019 if all([item is None for item in newcontent]): +1020 raise Exception("Division returns completely undefined correlator") +1021 return Corr(newcontent) +1022 +1023 elif isinstance(y, (Obs, CObs)): +1024 if isinstance(y, Obs): +1025 if y.value == 0: +1026 raise Exception('Division by zero will return undefined correlator') +1027 if isinstance(y, CObs): +1028 if y.is_zero(): +1029 raise Exception('Division by zero will return undefined correlator') +1030 +1031 newcontent = [] +1032 for t in range(self.T): +1033 if _check_for_none(self, self.content[t]): +1034 newcontent.append(None) +1035 else: +1036 newcontent.append(self.content[t] / y) +1037 return Corr(newcontent, prange=self.prange) +1038 +1039 elif isinstance(y, (int, float)): +1040 if y == 0: +1041 raise Exception('Division by zero will return undefined correlator') +1042 newcontent = [] +1043 for t in range(self.T): +1044 if _check_for_none(self, self.content[t]): +1045 newcontent.append(None) +1046 else: +1047 newcontent.append(self.content[t] / y) +1048 return Corr(newcontent, prange=self.prange) +1049 elif isinstance(y, np.ndarray): +1050 if y.shape == (self.T,): +1051 return Corr(list((np.array(self.content).T / y).T)) +1052 else: +1053 raise ValueError("operands could not be broadcast together") +1054 else: +1055 raise TypeError('Corr / wrong type') 1056 -1057 def __sub__(self, y): -1058 return self + (-y) -1059 -1060 def __pow__(self, y): -1061 if isinstance(y, (Obs, int, float, CObs)): -1062 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] -1063 return Corr(newcontent, prange=self.prange) -1064 else: -1065 raise TypeError('Type of exponent not supported') -1066 -1067 def __abs__(self): -1068 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] -1069 return Corr(newcontent, prange=self.prange) +1057 def __neg__(self): +1058 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] +1059 return Corr(newcontent, prange=self.prange) +1060 +1061 def __sub__(self, y): +1062 return self + (-y) +1063 +1064 def __pow__(self, y): +1065 if isinstance(y, (Obs, int, float, CObs)): +1066 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] +1067 return Corr(newcontent, prange=self.prange) +1068 else: +1069 raise TypeError('Type of exponent not supported') 1070 -1071 # The numpy functions: -1072 def sqrt(self): -1073 return self ** 0.5 +1071 def __abs__(self): +1072 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] +1073 return Corr(newcontent, prange=self.prange) 1074 -1075 def log(self): -1076 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] -1077 return Corr(newcontent, prange=self.prange) +1075 # The numpy functions: +1076 def sqrt(self): +1077 return self ** 0.5 1078 -1079 def exp(self): -1080 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] +1079 def log(self): +1080 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] 1081 return Corr(newcontent, prange=self.prange) 1082 -1083 def _apply_func_to_corr(self, func): -1084 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] -1085 for t in range(self.T): -1086 if _check_for_none(self, newcontent[t]): -1087 continue -1088 if np.isnan(np.sum(newcontent[t]).value): -1089 newcontent[t] = None -1090 if all([item is None for item in newcontent]): -1091 raise Exception('Operation returns undefined correlator') -1092 return Corr(newcontent) -1093 -1094 def sin(self): -1095 return self._apply_func_to_corr(np.sin) -1096 -1097 def cos(self): -1098 return self._apply_func_to_corr(np.cos) -1099 -1100 def tan(self): -1101 return self._apply_func_to_corr(np.tan) -1102 -1103 def sinh(self): -1104 return self._apply_func_to_corr(np.sinh) -1105 -1106 def cosh(self): -1107 return self._apply_func_to_corr(np.cosh) -1108 -1109 def tanh(self): -1110 return self._apply_func_to_corr(np.tanh) -1111 -1112 def arcsin(self): -1113 return self._apply_func_to_corr(np.arcsin) -1114 -1115 def arccos(self): -1116 return self._apply_func_to_corr(np.arccos) -1117 -1118 def arctan(self): -1119 return self._apply_func_to_corr(np.arctan) -1120 -1121 def arcsinh(self): -1122 return self._apply_func_to_corr(np.arcsinh) -1123 -1124 def arccosh(self): -1125 return self._apply_func_to_corr(np.arccosh) -1126 -1127 def arctanh(self): -1128 return self._apply_func_to_corr(np.arctanh) -1129 -1130 # Right hand side operations (require tweak in main module to work) -1131 def __radd__(self, y): -1132 return self + y +1083 def exp(self): +1084 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] +1085 return Corr(newcontent, prange=self.prange) +1086 +1087 def _apply_func_to_corr(self, func): +1088 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] +1089 for t in range(self.T): +1090 if _check_for_none(self, newcontent[t]): +1091 continue +1092 if np.isnan(np.sum(newcontent[t]).value): +1093 newcontent[t] = None +1094 if all([item is None for item in newcontent]): +1095 raise Exception('Operation returns undefined correlator') +1096 return Corr(newcontent) +1097 +1098 def sin(self): +1099 return self._apply_func_to_corr(np.sin) +1100 +1101 def cos(self): +1102 return self._apply_func_to_corr(np.cos) +1103 +1104 def tan(self): +1105 return self._apply_func_to_corr(np.tan) +1106 +1107 def sinh(self): +1108 return self._apply_func_to_corr(np.sinh) +1109 +1110 def cosh(self): +1111 return self._apply_func_to_corr(np.cosh) +1112 +1113 def tanh(self): +1114 return self._apply_func_to_corr(np.tanh) +1115 +1116 def arcsin(self): +1117 return self._apply_func_to_corr(np.arcsin) +1118 +1119 def arccos(self): +1120 return self._apply_func_to_corr(np.arccos) +1121 +1122 def arctan(self): +1123 return self._apply_func_to_corr(np.arctan) +1124 +1125 def arcsinh(self): +1126 return self._apply_func_to_corr(np.arcsinh) +1127 +1128 def arccosh(self): +1129 return self._apply_func_to_corr(np.arccosh) +1130 +1131 def arctanh(self): +1132 return self._apply_func_to_corr(np.arctanh) 1133 -1134 def __rsub__(self, y): -1135 return -self + y -1136 -1137 def __rmul__(self, y): -1138 return self * y -1139 -1140 def __rtruediv__(self, y): -1141 return (self / y) ** (-1) -1142 -1143 @property -1144 def real(self): -1145 def return_real(obs_OR_cobs): -1146 if isinstance(obs_OR_cobs, CObs): -1147 return obs_OR_cobs.real -1148 else: -1149 return obs_OR_cobs -1150 -1151 return self._apply_func_to_corr(return_real) -1152 -1153 @property -1154 def imag(self): -1155 def return_imag(obs_OR_cobs): -1156 if isinstance(obs_OR_cobs, CObs): -1157 return obs_OR_cobs.imag -1158 else: -1159 return obs_OR_cobs * 0 # So it stays the right type -1160 -1161 return self._apply_func_to_corr(return_imag) -1162 -1163 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): -1164 r''' Project large correlation matrix to lowest states -1165 -1166 This method can be used to reduce the size of an (N x N) correlation matrix -1167 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise -1168 is still small. +1134 # Right hand side operations (require tweak in main module to work) +1135 def __radd__(self, y): +1136 return self + y +1137 +1138 def __rsub__(self, y): +1139 return -self + y +1140 +1141 def __rmul__(self, y): +1142 return self * y +1143 +1144 def __rtruediv__(self, y): +1145 return (self / y) ** (-1) +1146 +1147 @property +1148 def real(self): +1149 def return_real(obs_OR_cobs): +1150 if isinstance(obs_OR_cobs, CObs): +1151 return obs_OR_cobs.real +1152 else: +1153 return obs_OR_cobs +1154 +1155 return self._apply_func_to_corr(return_real) +1156 +1157 @property +1158 def imag(self): +1159 def return_imag(obs_OR_cobs): +1160 if isinstance(obs_OR_cobs, CObs): +1161 return obs_OR_cobs.imag +1162 else: +1163 return obs_OR_cobs * 0 # So it stays the right type +1164 +1165 return self._apply_func_to_corr(return_imag) +1166 +1167 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): +1168 r''' Project large correlation matrix to lowest states 1169 -1170 Parameters -1171 ---------- -1172 Ntrunc: int -1173 Rank of the target matrix. -1174 tproj: int -1175 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. -1176 The default value is 3. -1177 t0proj: int -1178 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly -1179 discouraged for O(a) improved theories, since the correctness of the procedure -1180 cannot be granted in this case. The default value is 2. -1181 basematrix : Corr -1182 Correlation matrix that is used to determine the eigenvectors of the -1183 lowest states based on a GEVP. basematrix is taken to be the Corr itself if -1184 is is not specified. -1185 -1186 Notes -1187 ----- -1188 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving -1189 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}$ -1190 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the -1191 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via -1192 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large -1193 correlation matrix and to remove some noise that is added by irrelevant operators. -1194 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated -1195 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. -1196 ''' -1197 -1198 if self.N == 1: -1199 raise Exception('Method cannot be applied to one-dimensional correlators.') -1200 if basematrix is None: -1201 basematrix = self -1202 if Ntrunc >= basematrix.N: -1203 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) -1204 if basematrix.N != self.N: -1205 raise Exception('basematrix and targetmatrix have to be of the same size.') -1206 -1207 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] -1208 -1209 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) -1210 rmat = [] -1211 for t in range(basematrix.T): -1212 for i in range(Ntrunc): -1213 for j in range(Ntrunc): -1214 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] -1215 rmat.append(np.copy(tmpmat)) -1216 -1217 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] -1218 return Corr(newcontent) +1170 This method can be used to reduce the size of an (N x N) correlation matrix +1171 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise +1172 is still small. +1173 +1174 Parameters +1175 ---------- +1176 Ntrunc: int +1177 Rank of the target matrix. +1178 tproj: int +1179 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. +1180 The default value is 3. +1181 t0proj: int +1182 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly +1183 discouraged for O(a) improved theories, since the correctness of the procedure +1184 cannot be granted in this case. The default value is 2. +1185 basematrix : Corr +1186 Correlation matrix that is used to determine the eigenvectors of the +1187 lowest states based on a GEVP. basematrix is taken to be the Corr itself if +1188 is is not specified. +1189 +1190 Notes +1191 ----- +1192 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving +1193 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}$ +1194 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the +1195 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via +1196 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large +1197 correlation matrix and to remove some noise that is added by irrelevant operators. +1198 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated +1199 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. +1200 ''' +1201 +1202 if self.N == 1: +1203 raise Exception('Method cannot be applied to one-dimensional correlators.') +1204 if basematrix is None: +1205 basematrix = self +1206 if Ntrunc >= basematrix.N: +1207 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) +1208 if basematrix.N != self.N: +1209 raise Exception('basematrix and targetmatrix have to be of the same size.') +1210 +1211 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] +1212 +1213 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) +1214 rmat = [] +1215 for t in range(basematrix.T): +1216 for i in range(Ntrunc): +1217 for j in range(Ntrunc): +1218 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] +1219 rmat.append(np.copy(tmpmat)) +1220 +1221 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] +1222 return Corr(newcontent) @@ -3775,49 +3783,53 @@ Available choice: symmetric, improved, default: symmetric 605 for t in range(self.T - 1): 606 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): 607 newcontent.append(None) -608 else: -609 newcontent.append(self.content[t] / self.content[t + 1]) -610 if(all([x is None for x in newcontent])): -611 raise Exception('m_eff is undefined at all timeslices') -612 -613 return np.log(Corr(newcontent, padding=[0, 1])) +608 elif self.content[t][0].value / self.content[t + 1][0].value < 0: +609 newcontent.append(None) +610 else: +611 newcontent.append(self.content[t] / self.content[t + 1]) +612 if(all([x is None for x in newcontent])): +613 raise Exception('m_eff is undefined at all timeslices') 614 -615 elif variant in ['periodic', 'cosh', 'sinh']: -616 if variant in ['periodic', 'cosh']: -617 func = anp.cosh -618 else: -619 func = anp.sinh -620 -621 def root_function(x, d): -622 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d -623 -624 newcontent = [] -625 for t in range(self.T - 1): -626 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): -627 newcontent.append(None) -628 # Fill the two timeslices in the middle of the lattice with their predecessors -629 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: -630 newcontent.append(newcontent[-1]) -631 else: -632 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) -633 if(all([x is None for x in newcontent])): -634 raise Exception('m_eff is undefined at all timeslices') -635 -636 return Corr(newcontent, padding=[0, 1]) -637 -638 elif variant == 'arccosh': -639 newcontent = [] -640 for t in range(1, self.T - 1): -641 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): -642 newcontent.append(None) -643 else: -644 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) -645 if(all([x is None for x in newcontent])): -646 raise Exception("m_eff is undefined at all timeslices") -647 return np.arccosh(Corr(newcontent, padding=[1, 1])) -648 -649 else: -650 raise Exception('Unknown variant.') +615 return np.log(Corr(newcontent, padding=[0, 1])) +616 +617 elif variant in ['periodic', 'cosh', 'sinh']: +618 if variant in ['periodic', 'cosh']: +619 func = anp.cosh +620 else: +621 func = anp.sinh +622 +623 def root_function(x, d): +624 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d +625 +626 newcontent = [] +627 for t in range(self.T - 1): +628 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): +629 newcontent.append(None) +630 # Fill the two timeslices in the middle of the lattice with their predecessors +631 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: +632 newcontent.append(newcontent[-1]) +633 elif self.content[t][0].value / self.content[t + 1][0].value < 0: +634 newcontent.append(None) +635 else: +636 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) +637 if(all([x is None for x in newcontent])): +638 raise Exception('m_eff is undefined at all timeslices') +639 +640 return Corr(newcontent, padding=[0, 1]) +641 +642 elif variant == 'arccosh': +643 newcontent = [] +644 for t in range(1, self.T - 1): +645 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): +646 newcontent.append(None) +647 else: +648 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) +649 if(all([x is None for x in newcontent])): +650 raise Exception("m_eff is undefined at all timeslices") +651 return np.arccosh(Corr(newcontent, padding=[1, 1])) +652 +653 else: +654 raise Exception('Unknown variant.') @@ -3850,39 +3862,39 @@ guess for the root finder, only relevant for the root variant - 
652    def fit(self, function, fitrange=None, silent=False, **kwargs):
-653        r'''Fits function to the data
-654
-655        Parameters
-656        ----------
-657        function : obj
-658            function to fit to the data. See fits.least_squares for details.
-659        fitrange : list
-660            Two element list containing the timeslices on which the fit is supposed to start and stop.
-661            Caution: This range is inclusive as opposed to standard python indexing.
-662            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
-663            If not specified, self.prange or all timeslices are used.
-664        silent : bool
-665            Decides whether output is printed to the standard output.
-666        '''
-667        if self.N != 1:
-668            raise Exception("Correlator must be projected before fitting")
-669
-670        if fitrange is None:
-671            if self.prange:
-672                fitrange = self.prange
-673            else:
-674                fitrange = [0, self.T - 1]
-675        else:
-676            if not isinstance(fitrange, list):
-677                raise Exception("fitrange has to be a list with two elements")
-678            if len(fitrange) != 2:
-679                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
-680
-681        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-682        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-683        result = least_squares(xs, ys, function, silent=silent, **kwargs)
-684        return result
+            
656    def fit(self, function, fitrange=None, silent=False, **kwargs):
+657        r'''Fits function to the data
+658
+659        Parameters
+660        ----------
+661        function : obj
+662            function to fit to the data. See fits.least_squares for details.
+663        fitrange : list
+664            Two element list containing the timeslices on which the fit is supposed to start and stop.
+665            Caution: This range is inclusive as opposed to standard python indexing.
+666            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+667            If not specified, self.prange or all timeslices are used.
+668        silent : bool
+669            Decides whether output is printed to the standard output.
+670        '''
+671        if self.N != 1:
+672            raise Exception("Correlator must be projected before fitting")
+673
+674        if fitrange is None:
+675            if self.prange:
+676                fitrange = self.prange
+677            else:
+678                fitrange = [0, self.T - 1]
+679        else:
+680            if not isinstance(fitrange, list):
+681                raise Exception("fitrange has to be a list with two elements")
+682            if len(fitrange) != 2:
+683                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
+684
+685        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+686        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+687        result = least_squares(xs, ys, function, silent=silent, **kwargs)
+688        return result
 

 
 
@@ -3916,42 +3928,42 @@ Decides whether output is printed to the standard output.
- 
686    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
-687        """ Extract a plateau value from a Corr object
-688
-689        Parameters
-690        ----------
-691        plateau_range : list
-692            list with two entries, indicating the first and the last timeslice
-693            of the plateau region.
-694        method : str
-695            method to extract the plateau.
-696                'fit' fits a constant to the plateau region
-697                'avg', 'average' or 'mean' just average over the given timeslices.
-698        auto_gamma : bool
-699            apply gamma_method with default parameters to the Corr. Defaults to None
-700        """
-701        if not plateau_range:
-702            if self.prange:
-703                plateau_range = self.prange
-704            else:
-705                raise Exception("no plateau range provided")
-706        if self.N != 1:
-707            raise Exception("Correlator must be projected before getting a plateau.")
-708        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
-709            raise Exception("plateau is undefined at all timeslices in plateaurange.")
-710        if auto_gamma:
-711            self.gamma_method()
-712        if method == "fit":
-713            def const_func(a, t):
-714                return a[0]
-715            return self.fit(const_func, plateau_range)[0]
-716        elif method in ["avg", "average", "mean"]:
-717            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
-718            return returnvalue
-719
-720        else:
-721            raise Exception("Unsupported plateau method: " + method)
+            
690    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+691        """ Extract a plateau value from a Corr object
+692
+693        Parameters
+694        ----------
+695        plateau_range : list
+696            list with two entries, indicating the first and the last timeslice
+697            of the plateau region.
+698        method : str
+699            method to extract the plateau.
+700                'fit' fits a constant to the plateau region
+701                'avg', 'average' or 'mean' just average over the given timeslices.
+702        auto_gamma : bool
+703            apply gamma_method with default parameters to the Corr. Defaults to None
+704        """
+705        if not plateau_range:
+706            if self.prange:
+707                plateau_range = self.prange
+708            else:
+709                raise Exception("no plateau range provided")
+710        if self.N != 1:
+711            raise Exception("Correlator must be projected before getting a plateau.")
+712        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+713            raise Exception("plateau is undefined at all timeslices in plateaurange.")
+714        if auto_gamma:
+715            self.gamma_method()
+716        if method == "fit":
+717            def const_func(a, t):
+718                return a[0]
+719            return self.fit(const_func, plateau_range)[0]
+720        elif method in ["avg", "average", "mean"]:
+721            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+722            return returnvalue
+723
+724        else:
+725            raise Exception("Unsupported plateau method: " + method)
 

 
 
@@ -3985,17 +3997,17 @@ apply gamma_method with default parameters to the Corr. Defaults to None
- 
723    def set_prange(self, prange):
-724        """Sets the attribute prange of the Corr object."""
-725        if not len(prange) == 2:
-726            raise Exception("prange must be a list or array with two values")
-727        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
-728            raise Exception("Start and end point must be integers")
-729        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
-730            raise Exception("Start and end point must define a range in the interval 0,T")
-731
-732        self.prange = prange
-733        return
+            
727    def set_prange(self, prange):
+728        """Sets the attribute prange of the Corr object."""
+729        if not len(prange) == 2:
+730            raise Exception("prange must be a list or array with two values")
+731        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+732            raise Exception("Start and end point must be integers")
+733        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
+734            raise Exception("Start and end point must define a range in the interval 0,T")
+735
+736        self.prange = prange
+737        return
 

 
 
@@ -4015,124 +4027,124 @@ apply gamma_method with default parameters to the Corr. Defaults to None
- 
735    def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
-736        """Plots the correlator using the tag of the correlator as label if available.
-737
-738        Parameters
-739        ----------
-740        x_range : list
-741            list of two values, determining the range of the x-axis e.g. [4, 8].
-742        comp : Corr or list of Corr
-743            Correlator or list of correlators which are plotted for comparison.
-744            The tags of these correlators are used as labels if available.
-745        logscale : bool
-746            Sets y-axis to logscale.
-747        plateau : Obs
-748            Plateau value to be visualized in the figure.
-749        fit_res : Fit_result
-750            Fit_result object to be visualized.
-751        ylabel : str
-752            Label for the y-axis.
-753        save : str
-754            path to file in which the figure should be saved.
-755        auto_gamma : bool
-756            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
-757        hide_sigma : float
-758            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
-759        references : list
-760            List of floating point values that are displayed as horizontal lines for reference.
-761        title : string
-762            Optional title of the figure.
-763        """
-764        if self.N != 1:
-765            raise Exception("Correlator must be projected before plotting")
-766
-767        if auto_gamma:
-768            self.gamma_method()
-769
-770        if x_range is None:
-771            x_range = [0, self.T - 1]
-772
-773        fig = plt.figure()
-774        ax1 = fig.add_subplot(111)
-775
-776        x, y, y_err = self.plottable()
-777        if hide_sigma:
-778            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
-779        else:
-780            hide_from = None
-781        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
-782        if logscale:
-783            ax1.set_yscale('log')
-784        else:
-785            if y_range is None:
-786                try:
-787                    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)])
-788                    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)])
-789                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
-790                except Exception:
-791                    pass
-792            else:
-793                ax1.set_ylim(y_range)
-794        if comp:
-795            if isinstance(comp, (Corr, list)):
-796                for corr in comp if isinstance(comp, list) else [comp]:
-797                    if auto_gamma:
-798                        corr.gamma_method()
-799                    x, y, y_err = corr.plottable()
-800                    if hide_sigma:
-801                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
-802                    else:
-803                        hide_from = None
-804                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
-805            else:
-806                raise Exception("'comp' must be a correlator or a list of correlators.")
-807
-808        if plateau:
-809            if isinstance(plateau, Obs):
-810                if auto_gamma:
-811                    plateau.gamma_method()
-812                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
-813                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
-814            else:
-815                raise Exception("'plateau' must be an Obs")
-816
-817        if references:
-818            if isinstance(references, list):
-819                for ref in references:
-820                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
-821            else:
-822                raise Exception("'references' must be a list of floating pint values.")
-823
-824        if self.prange:
-825            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
-826            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+            
739    def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None, title=None):
+740        """Plots the correlator using the tag of the correlator as label if available.
+741
+742        Parameters
+743        ----------
+744        x_range : list
+745            list of two values, determining the range of the x-axis e.g. [4, 8].
+746        comp : Corr or list of Corr
+747            Correlator or list of correlators which are plotted for comparison.
+748            The tags of these correlators are used as labels if available.
+749        logscale : bool
+750            Sets y-axis to logscale.
+751        plateau : Obs
+752            Plateau value to be visualized in the figure.
+753        fit_res : Fit_result
+754            Fit_result object to be visualized.
+755        ylabel : str
+756            Label for the y-axis.
+757        save : str
+758            path to file in which the figure should be saved.
+759        auto_gamma : bool
+760            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+761        hide_sigma : float
+762            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+763        references : list
+764            List of floating point values that are displayed as horizontal lines for reference.
+765        title : string
+766            Optional title of the figure.
+767        """
+768        if self.N != 1:
+769            raise Exception("Correlator must be projected before plotting")
+770
+771        if auto_gamma:
+772            self.gamma_method()
+773
+774        if x_range is None:
+775            x_range = [0, self.T - 1]
+776
+777        fig = plt.figure()
+778        ax1 = fig.add_subplot(111)
+779
+780        x, y, y_err = self.plottable()
+781        if hide_sigma:
+782            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+783        else:
+784            hide_from = None
+785        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+786        if logscale:
+787            ax1.set_yscale('log')
+788        else:
+789            if y_range is None:
+790                try:
+791                    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)])
+792                    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)])
+793                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+794                except Exception:
+795                    pass
+796            else:
+797                ax1.set_ylim(y_range)
+798        if comp:
+799            if isinstance(comp, (Corr, list)):
+800                for corr in comp if isinstance(comp, list) else [comp]:
+801                    if auto_gamma:
+802                        corr.gamma_method()
+803                    x, y, y_err = corr.plottable()
+804                    if hide_sigma:
+805                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+806                    else:
+807                        hide_from = None
+808                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+809            else:
+810                raise Exception("'comp' must be a correlator or a list of correlators.")
+811
+812        if plateau:
+813            if isinstance(plateau, Obs):
+814                if auto_gamma:
+815                    plateau.gamma_method()
+816                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+817                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+818            else:
+819                raise Exception("'plateau' must be an Obs")
+820
+821        if references:
+822            if isinstance(references, list):
+823                for ref in references:
+824                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+825            else:
+826                raise Exception("'references' must be a list of floating pint values.")
 827
-828        if fit_res:
-829            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
-830            ax1.plot(x_samples,
-831                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
-832                     ls='-', marker=',', lw=2)
-833
-834        ax1.set_xlabel(r'$x_0 / a$')
-835        if ylabel:
-836            ax1.set_ylabel(ylabel)
-837        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
-838
-839        handles, labels = ax1.get_legend_handles_labels()
-840        if labels:
-841            ax1.legend()
+828        if self.prange:
+829            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
+830            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+831
+832        if fit_res:
+833            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+834            ax1.plot(x_samples,
+835                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
+836                     ls='-', marker=',', lw=2)
+837
+838        ax1.set_xlabel(r'$x_0 / a$')
+839        if ylabel:
+840            ax1.set_ylabel(ylabel)
+841        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
 842
-843        if title:
-844            plt.title(title)
-845
-846        plt.draw()
-847
-848        if save:
-849            if isinstance(save, str):
-850                fig.savefig(save, bbox_inches='tight')
-851            else:
-852                raise Exception("'save' has to be a string.")
+843        handles, labels = ax1.get_legend_handles_labels()
+844        if labels:
+845            ax1.legend()
+846
+847        if title:
+848            plt.title(title)
+849
+850        plt.draw()
+851
+852        if save:
+853            if isinstance(save, str):
+854                fig.savefig(save, bbox_inches='tight')
+855            else:
+856                raise Exception("'save' has to be a string.")
 

 
 
@@ -4180,34 +4192,34 @@ Optional title of the figure.
- 
854    def spaghetti_plot(self, logscale=True):
-855        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-856
-857        Parameters
-858        ----------
-859        logscale : bool
-860            Determines whether the scale of the y-axis is logarithmic or standard.
-861        """
-862        if self.N != 1:
-863            raise Exception("Correlator needs to be projected first.")
-864
-865        mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
-866        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
-867
-868        for name in mc_names:
-869            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
-870
-871            fig = plt.figure()
-872            ax = fig.add_subplot(111)
-873            for dat in data:
-874                ax.plot(x0_vals, dat, ls='-', marker='')
-875
-876            if logscale is True:
-877                ax.set_yscale('log')
-878
-879            ax.set_xlabel(r'$x_0 / a$')
-880            plt.title(name)
-881            plt.draw()
+            
858    def spaghetti_plot(self, logscale=True):
+859        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+860
+861        Parameters
+862        ----------
+863        logscale : bool
+864            Determines whether the scale of the y-axis is logarithmic or standard.
+865        """
+866        if self.N != 1:
+867            raise Exception("Correlator needs to be projected first.")
+868
+869        mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
+870        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
+871
+872        for name in mc_names:
+873            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+874
+875            fig = plt.figure()
+876            ax = fig.add_subplot(111)
+877            for dat in data:
+878                ax.plot(x0_vals, dat, ls='-', marker='')
+879
+880            if logscale is True:
+881                ax.set_yscale('log')
+882
+883            ax.set_xlabel(r'$x_0 / a$')
+884            plt.title(name)
+885            plt.draw()
 

 
 
@@ -4234,29 +4246,29 @@ Determines whether the scale of the y-axis is logarithmic or standard.
- 
883    def dump(self, filename, datatype="json.gz", **kwargs):
-884        """Dumps the Corr into a file of chosen type
-885        Parameters
-886        ----------
-887        filename : str
-888            Name of the file to be saved.
-889        datatype : str
-890            Format of the exported file. Supported formats include
-891            "json.gz" and "pickle"
-892        path : str
-893            specifies a custom path for the file (default '.')
-894        """
-895        if datatype == "json.gz":
-896            from .input.json import dump_to_json
-897            if 'path' in kwargs:
-898                file_name = kwargs.get('path') + '/' + filename
-899            else:
-900                file_name = filename
-901            dump_to_json(self, file_name)
-902        elif datatype == "pickle":
-903            dump_object(self, filename, **kwargs)
-904        else:
-905            raise Exception("Unknown datatype " + str(datatype))
+            
887    def dump(self, filename, datatype="json.gz", **kwargs):
+888        """Dumps the Corr into a file of chosen type
+889        Parameters
+890        ----------
+891        filename : str
+892            Name of the file to be saved.
+893        datatype : str
+894            Format of the exported file. Supported formats include
+895            "json.gz" and "pickle"
+896        path : str
+897            specifies a custom path for the file (default '.')
+898        """
+899        if datatype == "json.gz":
+900            from .input.json import dump_to_json
+901            if 'path' in kwargs:
+902                file_name = kwargs.get('path') + '/' + filename
+903            else:
+904                file_name = filename
+905            dump_to_json(self, file_name)
+906        elif datatype == "pickle":
+907            dump_object(self, filename, **kwargs)
+908        else:
+909            raise Exception("Unknown datatype " + str(datatype))
 

 
 
@@ -4288,8 +4300,8 @@ specifies a custom path for the file (default '.')
- 
907    def print(self, print_range=None):
-908        print(self.__repr__(print_range))
+            
911    def print(self, print_range=None):
+912        print(self.__repr__(print_range))
 

 
 
@@ -4307,8 +4319,8 @@ specifies a custom path for the file (default '.')
- 
1072    def sqrt(self):
-1073        return self ** 0.5
+            
1076    def sqrt(self):
+1077        return self ** 0.5
 

 
 
@@ -4326,9 +4338,9 @@ specifies a custom path for the file (default '.')
- 
1075    def log(self):
-1076        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
-1077        return Corr(newcontent, prange=self.prange)
+            
1079    def log(self):
+1080        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
+1081        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4346,9 +4358,9 @@ specifies a custom path for the file (default '.')
- 
1079    def exp(self):
-1080        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
-1081        return Corr(newcontent, prange=self.prange)
+            
1083    def exp(self):
+1084        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
+1085        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4366,8 +4378,8 @@ specifies a custom path for the file (default '.')
- 
1094    def sin(self):
-1095        return self._apply_func_to_corr(np.sin)
+            
1098    def sin(self):
+1099        return self._apply_func_to_corr(np.sin)
 

 
 
@@ -4385,8 +4397,8 @@ specifies a custom path for the file (default '.')
- 
1097    def cos(self):
-1098        return self._apply_func_to_corr(np.cos)
+            
1101    def cos(self):
+1102        return self._apply_func_to_corr(np.cos)
 

 
 
@@ -4404,8 +4416,8 @@ specifies a custom path for the file (default '.')
- 
1100    def tan(self):
-1101        return self._apply_func_to_corr(np.tan)
+            
1104    def tan(self):
+1105        return self._apply_func_to_corr(np.tan)
 

 
 
@@ -4423,8 +4435,8 @@ specifies a custom path for the file (default '.')
- 
1103    def sinh(self):
-1104        return self._apply_func_to_corr(np.sinh)
+            
1107    def sinh(self):
+1108        return self._apply_func_to_corr(np.sinh)
 

 
 
@@ -4442,8 +4454,8 @@ specifies a custom path for the file (default '.')
- 
1106    def cosh(self):
-1107        return self._apply_func_to_corr(np.cosh)
+            
1110    def cosh(self):
+1111        return self._apply_func_to_corr(np.cosh)
 

 
 
@@ -4461,8 +4473,8 @@ specifies a custom path for the file (default '.')
- 
1109    def tanh(self):
-1110        return self._apply_func_to_corr(np.tanh)
+            
1113    def tanh(self):
+1114        return self._apply_func_to_corr(np.tanh)
 

 
 
@@ -4480,8 +4492,8 @@ specifies a custom path for the file (default '.')
- 
1112    def arcsin(self):
-1113        return self._apply_func_to_corr(np.arcsin)
+            
1116    def arcsin(self):
+1117        return self._apply_func_to_corr(np.arcsin)
 

 
 
@@ -4499,8 +4511,8 @@ specifies a custom path for the file (default '.')
- 
1115    def arccos(self):
-1116        return self._apply_func_to_corr(np.arccos)
+            
1119    def arccos(self):
+1120        return self._apply_func_to_corr(np.arccos)
 

 
 
@@ -4518,8 +4530,8 @@ specifies a custom path for the file (default '.')
- 
1118    def arctan(self):
-1119        return self._apply_func_to_corr(np.arctan)
+            
1122    def arctan(self):
+1123        return self._apply_func_to_corr(np.arctan)
 

 
 
@@ -4537,8 +4549,8 @@ specifies a custom path for the file (default '.')
- 
1121    def arcsinh(self):
-1122        return self._apply_func_to_corr(np.arcsinh)
+            
1125    def arcsinh(self):
+1126        return self._apply_func_to_corr(np.arcsinh)
 

 
 
@@ -4556,8 +4568,8 @@ specifies a custom path for the file (default '.')
- 
1124    def arccosh(self):
-1125        return self._apply_func_to_corr(np.arccosh)
+            
1128    def arccosh(self):
+1129        return self._apply_func_to_corr(np.arccosh)
 

 
 
@@ -4575,8 +4587,8 @@ specifies a custom path for the file (default '.')
- 
1127    def arctanh(self):
-1128        return self._apply_func_to_corr(np.arctanh)
+            
1131    def arctanh(self):
+1132        return self._apply_func_to_corr(np.arctanh)
 

 
 
@@ -4616,62 +4628,62 @@ specifies a custom path for the file (default '.')
- 
1163    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1164        r''' Project large correlation matrix to lowest states
-1165
-1166        This method can be used to reduce the size of an (N x N) correlation matrix
-1167        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1168        is still small.
+            
1167    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1168        r''' Project large correlation matrix to lowest states
 1169
-1170        Parameters
-1171        ----------
-1172        Ntrunc: int
-1173            Rank of the target matrix.
-1174        tproj: int
-1175            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1176            The default value is 3.
-1177        t0proj: int
-1178            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1179            discouraged for O(a) improved theories, since the correctness of the procedure
-1180            cannot be granted in this case. The default value is 2.
-1181        basematrix : Corr
-1182            Correlation matrix that is used to determine the eigenvectors of the
-1183            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1184            is is not specified.
-1185
-1186        Notes
-1187        -----
-1188        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1189        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}$
-1190        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1191        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1192        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1193        correlation matrix and to remove some noise that is added by irrelevant operators.
-1194        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1195        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1196        '''
-1197
-1198        if self.N == 1:
-1199            raise Exception('Method cannot be applied to one-dimensional correlators.')
-1200        if basematrix is None:
-1201            basematrix = self
-1202        if Ntrunc >= basematrix.N:
-1203            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
-1204        if basematrix.N != self.N:
-1205            raise Exception('basematrix and targetmatrix have to be of the same size.')
-1206
-1207        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
-1208
-1209        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1210        rmat = []
-1211        for t in range(basematrix.T):
-1212            for i in range(Ntrunc):
-1213                for j in range(Ntrunc):
-1214                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1215            rmat.append(np.copy(tmpmat))
-1216
-1217        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
-1218        return Corr(newcontent)
+1170        This method can be used to reduce the size of an (N x N) correlation matrix
+1171        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1172        is still small.
+1173
+1174        Parameters
+1175        ----------
+1176        Ntrunc: int
+1177            Rank of the target matrix.
+1178        tproj: int
+1179            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1180            The default value is 3.
+1181        t0proj: int
+1182            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1183            discouraged for O(a) improved theories, since the correctness of the procedure
+1184            cannot be granted in this case. The default value is 2.
+1185        basematrix : Corr
+1186            Correlation matrix that is used to determine the eigenvectors of the
+1187            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1188            is is not specified.
+1189
+1190        Notes
+1191        -----
+1192        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1193        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}$
+1194        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1195        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1196        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1197        correlation matrix and to remove some noise that is added by irrelevant operators.
+1198        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1199        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1200        '''
+1201
+1202        if self.N == 1:
+1203            raise Exception('Method cannot be applied to one-dimensional correlators.')
+1204        if basematrix is None:
+1205            basematrix = self
+1206        if Ntrunc >= basematrix.N:
+1207            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
+1208        if basematrix.N != self.N:
+1209            raise Exception('basematrix and targetmatrix have to be of the same size.')
+1210
+1211        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
+1212
+1213        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1214        rmat = []
+1215        for t in range(basematrix.T):
+1216            for i in range(Ntrunc):
+1217                for j in range(Ntrunc):
+1218                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1219            rmat.append(np.copy(tmpmat))
+1220
+1221        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
+1222        return Corr(newcontent)