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
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)
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
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.")
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()
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))
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))
1072 def sqrt(self):
-1073 return self ** 0.5
+ 1076 def sqrt(self):
+1077 return self ** 0.5
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)