482 def deriv(self, variant="symmetric"):
-483 """Return the first derivative of the correlator with respect to x0.
-484
-485 Parameters
-486 ----------
-487 variant : str
-488 decides which definition of the finite differences derivative is used.
-489 Available choice: symmetric, forward, backward, improved, default: symmetric
-490 """
-491 if variant == "symmetric":
-492 newcontent = []
-493 for t in range(1, self.T - 1):
-494 if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-495 newcontent.append(None)
-496 else:
-497 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
-498 if(all([x is None for x in newcontent])):
-499 raise Exception('Derivative is undefined at all timeslices')
-500 return Corr(newcontent, padding=[1, 1])
-501 elif variant == "forward":
-502 newcontent = []
-503 for t in range(self.T - 1):
-504 if (self.content[t] is None) or (self.content[t + 1] is None):
-505 newcontent.append(None)
-506 else:
-507 newcontent.append(self.content[t + 1] - self.content[t])
-508 if(all([x is None for x in newcontent])):
-509 raise Exception("Derivative is undefined at all timeslices")
-510 return Corr(newcontent, padding=[0, 1])
-511 elif variant == "backward":
-512 newcontent = []
-513 for t in range(1, self.T):
-514 if (self.content[t - 1] is None) or (self.content[t] is None):
-515 newcontent.append(None)
-516 else:
-517 newcontent.append(self.content[t] - self.content[t - 1])
-518 if(all([x is None for x in newcontent])):
-519 raise Exception("Derivative is undefined at all timeslices")
-520 return Corr(newcontent, padding=[1, 0])
-521 elif variant == "improved":
-522 newcontent = []
-523 for t in range(2, self.T - 2):
-524 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
-525 newcontent.append(None)
-526 else:
-527 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
-528 if(all([x is None for x in newcontent])):
-529 raise Exception('Derivative is undefined at all timeslices')
-530 return Corr(newcontent, padding=[2, 2])
-531 else:
-532 raise Exception("Unknown variant.")
+ 484 def deriv(self, variant="symmetric"):
+485 """Return the first derivative of the correlator with respect to x0.
+486
+487 Parameters
+488 ----------
+489 variant : str
+490 decides which definition of the finite differences derivative is used.
+491 Available choice: symmetric, forward, backward, improved, default: symmetric
+492 """
+493 if variant == "symmetric":
+494 newcontent = []
+495 for t in range(1, self.T - 1):
+496 if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+497 newcontent.append(None)
+498 else:
+499 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
+500 if(all([x is None for x in newcontent])):
+501 raise Exception('Derivative is undefined at all timeslices')
+502 return Corr(newcontent, padding=[1, 1])
+503 elif variant == "forward":
+504 newcontent = []
+505 for t in range(self.T - 1):
+506 if (self.content[t] is None) or (self.content[t + 1] is None):
+507 newcontent.append(None)
+508 else:
+509 newcontent.append(self.content[t + 1] - self.content[t])
+510 if(all([x is None for x in newcontent])):
+511 raise Exception("Derivative is undefined at all timeslices")
+512 return Corr(newcontent, padding=[0, 1])
+513 elif variant == "backward":
+514 newcontent = []
+515 for t in range(1, self.T):
+516 if (self.content[t - 1] is None) or (self.content[t] is None):
+517 newcontent.append(None)
+518 else:
+519 newcontent.append(self.content[t] - self.content[t - 1])
+520 if(all([x is None for x in newcontent])):
+521 raise Exception("Derivative is undefined at all timeslices")
+522 return Corr(newcontent, padding=[1, 0])
+523 elif variant == "improved":
+524 newcontent = []
+525 for t in range(2, self.T - 2):
+526 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
+527 newcontent.append(None)
+528 else:
+529 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
+530 if(all([x is None for x in newcontent])):
+531 raise Exception('Derivative is undefined at all timeslices')
+532 return Corr(newcontent, padding=[2, 2])
+533 else:
+534 raise Exception("Unknown variant.")
534 def second_deriv(self, variant="symmetric"):
-535 """Return the second derivative of the correlator with respect to x0.
-536
-537 Parameters
-538 ----------
-539 variant : str
-540 decides which definition of the finite differences derivative is used.
-541 Available choice: symmetric, improved, default: symmetric
-542 """
-543 if variant == "symmetric":
-544 newcontent = []
-545 for t in range(1, self.T - 1):
-546 if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-547 newcontent.append(None)
-548 else:
-549 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
-550 if(all([x is None for x in newcontent])):
-551 raise Exception("Derivative is undefined at all timeslices")
-552 return Corr(newcontent, padding=[1, 1])
-553 elif variant == "improved":
-554 newcontent = []
-555 for t in range(2, self.T - 2):
-556 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
-557 newcontent.append(None)
-558 else:
-559 newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2]))
-560 if(all([x is None for x in newcontent])):
-561 raise Exception("Derivative is undefined at all timeslices")
-562 return Corr(newcontent, padding=[2, 2])
-563 else:
-564 raise Exception("Unknown variant.")
+ 536 def second_deriv(self, variant="symmetric"):
+537 """Return the second derivative of the correlator with respect to x0.
+538
+539 Parameters
+540 ----------
+541 variant : str
+542 decides which definition of the finite differences derivative is used.
+543 Available choice: symmetric, improved, default: symmetric
+544 """
+545 if variant == "symmetric":
+546 newcontent = []
+547 for t in range(1, self.T - 1):
+548 if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+549 newcontent.append(None)
+550 else:
+551 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
+552 if(all([x is None for x in newcontent])):
+553 raise Exception("Derivative is undefined at all timeslices")
+554 return Corr(newcontent, padding=[1, 1])
+555 elif variant == "improved":
+556 newcontent = []
+557 for t in range(2, self.T - 2):
+558 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
+559 newcontent.append(None)
+560 else:
+561 newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2]))
+562 if(all([x is None for x in newcontent])):
+563 raise Exception("Derivative is undefined at all timeslices")
+564 return Corr(newcontent, padding=[2, 2])
+565 else:
+566 raise Exception("Unknown variant.")
566 def m_eff(self, variant='log', guess=1.0):
-567 """Returns the effective mass of the correlator as correlator object
-568
-569 Parameters
-570 ----------
-571 variant : str
-572 log : uses the standard effective mass log(C(t) / C(t+1))
-573 cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
-574 sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
-575 See, e.g., arXiv:1205.5380
-576 arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
-577 guess : float
-578 guess for the root finder, only relevant for the root variant
-579 """
-580 if self.N != 1:
-581 raise Exception('Correlator must be projected before getting m_eff')
-582 if variant == 'log':
-583 newcontent = []
-584 for t in range(self.T - 1):
-585 if (self.content[t] is None) or (self.content[t + 1] is None):
-586 newcontent.append(None)
-587 else:
-588 newcontent.append(self.content[t] / self.content[t + 1])
-589 if(all([x is None for x in newcontent])):
-590 raise Exception('m_eff is undefined at all timeslices')
-591
-592 return np.log(Corr(newcontent, padding=[0, 1]))
+ 568 def m_eff(self, variant='log', guess=1.0):
+569 """Returns the effective mass of the correlator as correlator object
+570
+571 Parameters
+572 ----------
+573 variant : str
+574 log : uses the standard effective mass log(C(t) / C(t+1))
+575 cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
+576 sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
+577 See, e.g., arXiv:1205.5380
+578 arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
+579 guess : float
+580 guess for the root finder, only relevant for the root variant
+581 """
+582 if self.N != 1:
+583 raise Exception('Correlator must be projected before getting m_eff')
+584 if variant == 'log':
+585 newcontent = []
+586 for t in range(self.T - 1):
+587 if (self.content[t] is None) or (self.content[t + 1] is None):
+588 newcontent.append(None)
+589 else:
+590 newcontent.append(self.content[t] / self.content[t + 1])
+591 if(all([x is None for x in newcontent])):
+592 raise Exception('m_eff is undefined at all timeslices')
593
-594 elif variant in ['periodic', 'cosh', 'sinh']:
-595 if variant in ['periodic', 'cosh']:
-596 func = anp.cosh
-597 else:
-598 func = anp.sinh
-599
-600 def root_function(x, d):
-601 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
-602
-603 newcontent = []
-604 for t in range(self.T - 1):
-605 if (self.content[t] is None) or (self.content[t + 1] is None):
-606 newcontent.append(None)
-607 # Fill the two timeslices in the middle of the lattice with their predecessors
-608 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
-609 newcontent.append(newcontent[-1])
-610 else:
-611 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
-612 if(all([x is None for x in newcontent])):
-613 raise Exception('m_eff is undefined at all timeslices')
-614
-615 return Corr(newcontent, padding=[0, 1])
+594 return np.log(Corr(newcontent, padding=[0, 1]))
+595
+596 elif variant in ['periodic', 'cosh', 'sinh']:
+597 if variant in ['periodic', 'cosh']:
+598 func = anp.cosh
+599 else:
+600 func = anp.sinh
+601
+602 def root_function(x, d):
+603 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
+604
+605 newcontent = []
+606 for t in range(self.T - 1):
+607 if (self.content[t] is None) or (self.content[t + 1] is None):
+608 newcontent.append(None)
+609 # Fill the two timeslices in the middle of the lattice with their predecessors
+610 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+611 newcontent.append(newcontent[-1])
+612 else:
+613 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+614 if(all([x is None for x in newcontent])):
+615 raise Exception('m_eff is undefined at all timeslices')
616
-617 elif variant == 'arccosh':
-618 newcontent = []
-619 for t in range(1, self.T - 1):
-620 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
-621 newcontent.append(None)
-622 else:
-623 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
-624 if(all([x is None for x in newcontent])):
-625 raise Exception("m_eff is undefined at all timeslices")
-626 return np.arccosh(Corr(newcontent, padding=[1, 1]))
-627
-628 else:
-629 raise Exception('Unknown variant.')
+617 return Corr(newcontent, padding=[0, 1])
+618
+619 elif variant == 'arccosh':
+620 newcontent = []
+621 for t in range(1, self.T - 1):
+622 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
+623 newcontent.append(None)
+624 else:
+625 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+626 if(all([x is None for x in newcontent])):
+627 raise Exception("m_eff is undefined at all timeslices")
+628 return np.arccosh(Corr(newcontent, padding=[1, 1]))
+629
+630 else:
+631 raise Exception('Unknown variant.')
631 def fit(self, function, fitrange=None, silent=False, **kwargs):
-632 r'''Fits function to the data
-633
-634 Parameters
-635 ----------
-636 function : obj
-637 function to fit to the data. See fits.least_squares for details.
-638 fitrange : list
-639 Two element list containing the timeslices on which the fit is supposed to start and stop.
-640 Caution: This range is inclusive as opposed to standard python indexing.
-641 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
-642 If not specified, self.prange or all timeslices are used.
-643 silent : bool
-644 Decides whether output is printed to the standard output.
-645 '''
-646 if self.N != 1:
-647 raise Exception("Correlator must be projected before fitting")
-648
-649 if fitrange is None:
-650 if self.prange:
-651 fitrange = self.prange
-652 else:
-653 fitrange = [0, self.T - 1]
-654 else:
-655 if not isinstance(fitrange, list):
-656 raise Exception("fitrange has to be a list with two elements")
-657 if len(fitrange) != 2:
-658 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
-659
-660 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-661 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-662 result = least_squares(xs, ys, function, silent=silent, **kwargs)
-663 return result
+ 633 def fit(self, function, fitrange=None, silent=False, **kwargs):
+634 r'''Fits function to the data
+635
+636 Parameters
+637 ----------
+638 function : obj
+639 function to fit to the data. See fits.least_squares for details.
+640 fitrange : list
+641 Two element list containing the timeslices on which the fit is supposed to start and stop.
+642 Caution: This range is inclusive as opposed to standard python indexing.
+643 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+644 If not specified, self.prange or all timeslices are used.
+645 silent : bool
+646 Decides whether output is printed to the standard output.
+647 '''
+648 if self.N != 1:
+649 raise Exception("Correlator must be projected before fitting")
+650
+651 if fitrange is None:
+652 if self.prange:
+653 fitrange = self.prange
+654 else:
+655 fitrange = [0, self.T - 1]
+656 else:
+657 if not isinstance(fitrange, list):
+658 raise Exception("fitrange has to be a list with two elements")
+659 if len(fitrange) != 2:
+660 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
+661
+662 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+663 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+664 result = least_squares(xs, ys, function, silent=silent, **kwargs)
+665 return result
665 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
-666 """ Extract a plateau value from a Corr object
-667
-668 Parameters
-669 ----------
-670 plateau_range : list
-671 list with two entries, indicating the first and the last timeslice
-672 of the plateau region.
-673 method : str
-674 method to extract the plateau.
-675 'fit' fits a constant to the plateau region
-676 'avg', 'average' or 'mean' just average over the given timeslices.
-677 auto_gamma : bool
-678 apply gamma_method with default parameters to the Corr. Defaults to None
-679 """
-680 if not plateau_range:
-681 if self.prange:
-682 plateau_range = self.prange
-683 else:
-684 raise Exception("no plateau range provided")
-685 if self.N != 1:
-686 raise Exception("Correlator must be projected before getting a plateau.")
-687 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
-688 raise Exception("plateau is undefined at all timeslices in plateaurange.")
-689 if auto_gamma:
-690 self.gamma_method()
-691 if method == "fit":
-692 def const_func(a, t):
-693 return a[0]
-694 return self.fit(const_func, plateau_range)[0]
-695 elif method in ["avg", "average", "mean"]:
-696 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
-697 return returnvalue
-698
-699 else:
-700 raise Exception("Unsupported plateau method: " + method)
+ 667 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+668 """ Extract a plateau value from a Corr object
+669
+670 Parameters
+671 ----------
+672 plateau_range : list
+673 list with two entries, indicating the first and the last timeslice
+674 of the plateau region.
+675 method : str
+676 method to extract the plateau.
+677 'fit' fits a constant to the plateau region
+678 'avg', 'average' or 'mean' just average over the given timeslices.
+679 auto_gamma : bool
+680 apply gamma_method with default parameters to the Corr. Defaults to None
+681 """
+682 if not plateau_range:
+683 if self.prange:
+684 plateau_range = self.prange
+685 else:
+686 raise Exception("no plateau range provided")
+687 if self.N != 1:
+688 raise Exception("Correlator must be projected before getting a plateau.")
+689 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+690 raise Exception("plateau is undefined at all timeslices in plateaurange.")
+691 if auto_gamma:
+692 self.gamma_method()
+693 if method == "fit":
+694 def const_func(a, t):
+695 return a[0]
+696 return self.fit(const_func, plateau_range)[0]
+697 elif method in ["avg", "average", "mean"]:
+698 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+699 return returnvalue
+700
+701 else:
+702 raise Exception("Unsupported plateau method: " + method)
702 def set_prange(self, prange):
-703 """Sets the attribute prange of the Corr object."""
-704 if not len(prange) == 2:
-705 raise Exception("prange must be a list or array with two values")
-706 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
-707 raise Exception("Start and end point must be integers")
-708 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
-709 raise Exception("Start and end point must define a range in the interval 0,T")
-710
-711 self.prange = prange
-712 return
+ 704 def set_prange(self, prange):
+705 """Sets the attribute prange of the Corr object."""
+706 if not len(prange) == 2:
+707 raise Exception("prange must be a list or array with two values")
+708 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+709 raise Exception("Start and end point must be integers")
+710 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
+711 raise Exception("Start and end point must define a range in the interval 0,T")
+712
+713 self.prange = prange
+714 return
714 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):
-715 """Plots the correlator using the tag of the correlator as label if available.
-716
-717 Parameters
-718 ----------
-719 x_range : list
-720 list of two values, determining the range of the x-axis e.g. [4, 8]
-721 comp : Corr or list of Corr
-722 Correlator or list of correlators which are plotted for comparison.
-723 The tags of these correlators are used as labels if available.
-724 logscale : bool
-725 Sets y-axis to logscale
-726 plateau : Obs
-727 Plateau value to be visualized in the figure
-728 fit_res : Fit_result
-729 Fit_result object to be visualized
-730 ylabel : str
-731 Label for the y-axis
-732 save : str
-733 path to file in which the figure should be saved
-734 auto_gamma : bool
-735 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
-736 hide_sigma : float
-737 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
-738 references : list
-739 List of floating point values that are displayed as horizontal lines for reference.
-740 """
-741 if self.N != 1:
-742 raise Exception("Correlator must be projected before plotting")
-743
-744 if auto_gamma:
-745 self.gamma_method()
-746
-747 if x_range is None:
-748 x_range = [0, self.T - 1]
-749
-750 fig = plt.figure()
-751 ax1 = fig.add_subplot(111)
-752
-753 x, y, y_err = self.plottable()
-754 if hide_sigma:
-755 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
-756 else:
-757 hide_from = None
-758 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
-759 if logscale:
-760 ax1.set_yscale('log')
-761 else:
-762 if y_range is None:
-763 try:
-764 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)])
-765 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)])
-766 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
-767 except Exception:
-768 pass
-769 else:
-770 ax1.set_ylim(y_range)
-771 if comp:
-772 if isinstance(comp, (Corr, list)):
-773 for corr in comp if isinstance(comp, list) else [comp]:
-774 if auto_gamma:
-775 corr.gamma_method()
-776 x, y, y_err = corr.plottable()
-777 if hide_sigma:
-778 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
-779 else:
-780 hide_from = None
-781 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
-782 else:
-783 raise Exception("'comp' must be a correlator or a list of correlators.")
-784
-785 if plateau:
-786 if isinstance(plateau, Obs):
-787 if auto_gamma:
-788 plateau.gamma_method()
-789 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
-790 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
-791 else:
-792 raise Exception("'plateau' must be an Obs")
-793
-794 if references:
-795 if isinstance(references, list):
-796 for ref in references:
-797 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
-798 else:
-799 raise Exception("'references' must be a list of floating pint values.")
-800
-801 if self.prange:
-802 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
-803 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
-804
-805 if fit_res:
-806 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
-807 ax1.plot(x_samples,
-808 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
-809 ls='-', marker=',', lw=2)
-810
-811 ax1.set_xlabel(r'$x_0 / a$')
-812 if ylabel:
-813 ax1.set_ylabel(ylabel)
-814 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
-815
-816 handles, labels = ax1.get_legend_handles_labels()
-817 if labels:
-818 ax1.legend()
-819 plt.draw()
-820
-821 if save:
-822 if isinstance(save, str):
-823 fig.savefig(save)
-824 else:
-825 raise Exception("'save' has to be a string.")
+ 716 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):
+717 """Plots the correlator using the tag of the correlator as label if available.
+718
+719 Parameters
+720 ----------
+721 x_range : list
+722 list of two values, determining the range of the x-axis e.g. [4, 8]
+723 comp : Corr or list of Corr
+724 Correlator or list of correlators which are plotted for comparison.
+725 The tags of these correlators are used as labels if available.
+726 logscale : bool
+727 Sets y-axis to logscale
+728 plateau : Obs
+729 Plateau value to be visualized in the figure
+730 fit_res : Fit_result
+731 Fit_result object to be visualized
+732 ylabel : str
+733 Label for the y-axis
+734 save : str
+735 path to file in which the figure should be saved
+736 auto_gamma : bool
+737 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+738 hide_sigma : float
+739 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+740 references : list
+741 List of floating point values that are displayed as horizontal lines for reference.
+742 """
+743 if self.N != 1:
+744 raise Exception("Correlator must be projected before plotting")
+745
+746 if auto_gamma:
+747 self.gamma_method()
+748
+749 if x_range is None:
+750 x_range = [0, self.T - 1]
+751
+752 fig = plt.figure()
+753 ax1 = fig.add_subplot(111)
+754
+755 x, y, y_err = self.plottable()
+756 if hide_sigma:
+757 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
+758 else:
+759 hide_from = None
+760 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+761 if logscale:
+762 ax1.set_yscale('log')
+763 else:
+764 if y_range is None:
+765 try:
+766 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)])
+767 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)])
+768 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+769 except Exception:
+770 pass
+771 else:
+772 ax1.set_ylim(y_range)
+773 if comp:
+774 if isinstance(comp, (Corr, list)):
+775 for corr in comp if isinstance(comp, list) else [comp]:
+776 if auto_gamma:
+777 corr.gamma_method()
+778 x, y, y_err = corr.plottable()
+779 if hide_sigma:
+780 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
+781 else:
+782 hide_from = None
+783 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+784 else:
+785 raise Exception("'comp' must be a correlator or a list of correlators.")
+786
+787 if plateau:
+788 if isinstance(plateau, Obs):
+789 if auto_gamma:
+790 plateau.gamma_method()
+791 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+792 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+793 else:
+794 raise Exception("'plateau' must be an Obs")
+795
+796 if references:
+797 if isinstance(references, list):
+798 for ref in references:
+799 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+800 else:
+801 raise Exception("'references' must be a list of floating pint values.")
+802
+803 if self.prange:
+804 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
+805 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+806
+807 if fit_res:
+808 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+809 ax1.plot(x_samples,
+810 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
+811 ls='-', marker=',', lw=2)
+812
+813 ax1.set_xlabel(r'$x_0 / a$')
+814 if ylabel:
+815 ax1.set_ylabel(ylabel)
+816 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
+817
+818 handles, labels = ax1.get_legend_handles_labels()
+819 if labels:
+820 ax1.legend()
+821 plt.draw()
+822
+823 if save:
+824 if isinstance(save, str):
+825 fig.savefig(save)
+826 else:
+827 raise Exception("'save' has to be a string.")
827 def spaghetti_plot(self, logscale=True):
-828 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-829
-830 Parameters
-831 ----------
-832 logscale : bool
-833 Determines whether the scale of the y-axis is logarithmic or standard.
-834 """
-835 if self.N != 1:
-836 raise Exception("Correlator needs to be projected first.")
-837
-838 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]))
-839 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
-840
-841 for name in mc_names:
-842 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
-843
-844 fig = plt.figure()
-845 ax = fig.add_subplot(111)
-846 for dat in data:
-847 ax.plot(x0_vals, dat, ls='-', marker='')
-848
-849 if logscale is True:
-850 ax.set_yscale('log')
-851
-852 ax.set_xlabel(r'$x_0 / a$')
-853 plt.title(name)
-854 plt.draw()
+ 829 def spaghetti_plot(self, logscale=True):
+830 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+831
+832 Parameters
+833 ----------
+834 logscale : bool
+835 Determines whether the scale of the y-axis is logarithmic or standard.
+836 """
+837 if self.N != 1:
+838 raise Exception("Correlator needs to be projected first.")
+839
+840 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]))
+841 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
+842
+843 for name in mc_names:
+844 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+845
+846 fig = plt.figure()
+847 ax = fig.add_subplot(111)
+848 for dat in data:
+849 ax.plot(x0_vals, dat, ls='-', marker='')
+850
+851 if logscale is True:
+852 ax.set_yscale('log')
+853
+854 ax.set_xlabel(r'$x_0 / a$')
+855 plt.title(name)
+856 plt.draw()
856 def dump(self, filename, datatype="json.gz", **kwargs):
-857 """Dumps the Corr into a file of chosen type
-858 Parameters
-859 ----------
-860 filename : str
-861 Name of the file to be saved.
-862 datatype : str
-863 Format of the exported file. Supported formats include
-864 "json.gz" and "pickle"
-865 path : str
-866 specifies a custom path for the file (default '.')
-867 """
-868 if datatype == "json.gz":
-869 from .input.json import dump_to_json
-870 if 'path' in kwargs:
-871 file_name = kwargs.get('path') + '/' + filename
-872 else:
-873 file_name = filename
-874 dump_to_json(self, file_name)
-875 elif datatype == "pickle":
-876 dump_object(self, filename, **kwargs)
-877 else:
-878 raise Exception("Unknown datatype " + str(datatype))
+ 858 def dump(self, filename, datatype="json.gz", **kwargs):
+859 """Dumps the Corr into a file of chosen type
+860 Parameters
+861 ----------
+862 filename : str
+863 Name of the file to be saved.
+864 datatype : str
+865 Format of the exported file. Supported formats include
+866 "json.gz" and "pickle"
+867 path : str
+868 specifies a custom path for the file (default '.')
+869 """
+870 if datatype == "json.gz":
+871 from .input.json import dump_to_json
+872 if 'path' in kwargs:
+873 file_name = kwargs.get('path') + '/' + filename
+874 else:
+875 file_name = filename
+876 dump_to_json(self, file_name)
+877 elif datatype == "pickle":
+878 dump_object(self, filename, **kwargs)
+879 else:
+880 raise Exception("Unknown datatype " + str(datatype))
880 def print(self, range=[0, None]):
-881 print(self.__repr__(range))
+ 882 def print(self, range=[0, None]):
+883 print(self.__repr__(range))
1043 def sqrt(self):
-1044 return self**0.5
+ 1045 def sqrt(self):
+1046 return self**0.5
1046 def log(self):
-1047 newcontent = [None if (item is None) else np.log(item) for item in self.content]
-1048 return Corr(newcontent, prange=self.prange)
+ 1048 def log(self):
+1049 newcontent = [None if (item is None) else np.log(item) for item in self.content]
+1050 return Corr(newcontent, prange=self.prange)
1050 def exp(self):
-1051 newcontent = [None if (item is None) else np.exp(item) for item in self.content]
-1052 return Corr(newcontent, prange=self.prange)
+ 1052 def exp(self):
+1053 newcontent = [None if (item is None) else np.exp(item) for item in self.content]
+1054 return Corr(newcontent, prange=self.prange)
1065 def sin(self):
-1066 return self._apply_func_to_corr(np.sin)
+ 1067 def sin(self):
+1068 return self._apply_func_to_corr(np.sin)
1068 def cos(self):
-1069 return self._apply_func_to_corr(np.cos)
+ 1070 def cos(self):
+1071 return self._apply_func_to_corr(np.cos)
1071 def tan(self):
-1072 return self._apply_func_to_corr(np.tan)
+ 1073 def tan(self):
+1074 return self._apply_func_to_corr(np.tan)
1074 def sinh(self):
-1075 return self._apply_func_to_corr(np.sinh)
+ 1076 def sinh(self):
+1077 return self._apply_func_to_corr(np.sinh)
1077 def cosh(self):
-1078 return self._apply_func_to_corr(np.cosh)
+ 1079 def cosh(self):
+1080 return self._apply_func_to_corr(np.cosh)
1080 def tanh(self):
-1081 return self._apply_func_to_corr(np.tanh)
+ 1082 def tanh(self):
+1083 return self._apply_func_to_corr(np.tanh)
1083 def arcsin(self):
-1084 return self._apply_func_to_corr(np.arcsin)
+ 1085 def arcsin(self):
+1086 return self._apply_func_to_corr(np.arcsin)
1086 def arccos(self):
-1087 return self._apply_func_to_corr(np.arccos)
+ 1088 def arccos(self):
+1089 return self._apply_func_to_corr(np.arccos)
1089 def arctan(self):
-1090 return self._apply_func_to_corr(np.arctan)
+ 1091 def arctan(self):
+1092 return self._apply_func_to_corr(np.arctan)
1092 def arcsinh(self):
-1093 return self._apply_func_to_corr(np.arcsinh)
+ 1094 def arcsinh(self):
+1095 return self._apply_func_to_corr(np.arcsinh)
1095 def arccosh(self):
-1096 return self._apply_func_to_corr(np.arccosh)
+ 1097 def arccosh(self):
+1098 return self._apply_func_to_corr(np.arccosh)
1098 def arctanh(self):
-1099 return self._apply_func_to_corr(np.arctanh)
+ 1100 def arctanh(self):
+1101 return self._apply_func_to_corr(np.arctanh)
1134 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1135 r''' Project large correlation matrix to lowest states
-1136
-1137 This method can be used to reduce the size of an (N x N) correlation matrix
-1138 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1139 is still small.
-1140
-1141 Parameters
-1142 ----------
-1143 Ntrunc: int
-1144 Rank of the target matrix.
-1145 tproj: int
-1146 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1147 The default value is 3.
-1148 t0proj: int
-1149 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1150 discouraged for O(a) improved theories, since the correctness of the procedure
-1151 cannot be granted in this case. The default value is 2.
-1152 basematrix : Corr
-1153 Correlation matrix that is used to determine the eigenvectors of the
-1154 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1155 is is not specified.
-1156
-1157 Notes
-1158 -----
-1159 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1160 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}$
-1161 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1162 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1163 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1164 correlation matrix and to remove some noise that is added by irrelevant operators.
-1165 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1166 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1167 '''
-1168
-1169 if self.N == 1:
-1170 raise Exception('Method cannot be applied to one-dimensional correlators.')
-1171 if basematrix is None:
-1172 basematrix = self
-1173 if Ntrunc >= basematrix.N:
-1174 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
-1175 if basematrix.N != self.N:
-1176 raise Exception('basematrix and targetmatrix have to be of the same size.')
-1177
-1178 evecs = []
-1179 for i in range(Ntrunc):
-1180 evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
-1181
-1182 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1183 rmat = []
-1184 for t in range(basematrix.T):
-1185 for i in range(Ntrunc):
-1186 for j in range(Ntrunc):
-1187 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1188 rmat.append(np.copy(tmpmat))
-1189
-1190 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
-1191 return Corr(newcontent)
+ 1136 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1137 r''' Project large correlation matrix to lowest states
+1138
+1139 This method can be used to reduce the size of an (N x N) correlation matrix
+1140 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1141 is still small.
+1142
+1143 Parameters
+1144 ----------
+1145 Ntrunc: int
+1146 Rank of the target matrix.
+1147 tproj: int
+1148 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1149 The default value is 3.
+1150 t0proj: int
+1151 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1152 discouraged for O(a) improved theories, since the correctness of the procedure
+1153 cannot be granted in this case. The default value is 2.
+1154 basematrix : Corr
+1155 Correlation matrix that is used to determine the eigenvectors of the
+1156 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1157 is is not specified.
+1158
+1159 Notes
+1160 -----
+1161 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1162 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}$
+1163 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1164 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1165 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1166 correlation matrix and to remove some noise that is added by irrelevant operators.
+1167 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1168 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1169 '''
+1170
+1171 if self.N == 1:
+1172 raise Exception('Method cannot be applied to one-dimensional correlators.')
+1173 if basematrix is None:
+1174 basematrix = self
+1175 if Ntrunc >= basematrix.N:
+1176 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
+1177 if basematrix.N != self.N:
+1178 raise Exception('basematrix and targetmatrix have to be of the same size.')
+1179
+1180 evecs = []
+1181 for i in range(Ntrunc):
+1182 evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
+1183
+1184 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1185 rmat = []
+1186 for t in range(basematrix.T):
+1187 for i in range(Ntrunc):
+1188 for j in range(Ntrunc):
+1189 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1190 rmat.append(np.copy(tmpmat))
+1191
+1192 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
+1193 return Corr(newcontent)