From ea760795371398f74578dddfedb2d572d5fa1a71 Mon Sep 17 00:00:00 2001 From: fjosw Date: Wed, 30 Nov 2022 16:49:03 +0000 Subject: [PATCH] Documentation updated --- docs/pyerrors/correlators.html | 5377 ++++++++++++++++---------------- 1 file changed, 2690 insertions(+), 2687 deletions(-) diff --git a/docs/pyerrors/correlators.html b/docs/pyerrors/correlators.html index bf8c779b..753989dd 100644 --- a/docs/pyerrors/correlators.html +++ b/docs/pyerrors/correlators.html @@ -519,1010 +519,1011 @@ 306 else: 307 symmetric_corr = self.matrix_symmetric() 308 - 309 if sort is None: - 310 if (ts is None): - 311 raise Exception("ts is required if sort=None.") - 312 if (self.content[t0] is None) or (self.content[ts] is None): - 313 raise Exception("Corr not defined at t0/ts.") - 314 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) - 315 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts]) - 316 reordered_vecs = _GEVP_solver(Gt, G0) - 317 - 318 elif sort in ["Eigenvalue", "Eigenvector"]: - 319 if sort == "Eigenvalue" and ts is not None: - 320 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) - 321 all_vecs = [None] * (t0 + 1) - 322 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) - 323 for t in range(t0 + 1, self.T): - 324 try: - 325 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t]) - 326 all_vecs.append(_GEVP_solver(Gt, G0)) - 327 except Exception: - 328 all_vecs.append(None) - 329 if sort == "Eigenvector": - 330 if (ts is None): - 331 raise Exception("ts is required for the Eigenvector sorting method.") - 332 all_vecs = _sort_vectors(all_vecs, ts) - 333 - 334 reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] - 335 else: - 336 raise Exception("Unkown value for 'sort'.") - 337 - 338 if "state" in kwargs: - 339 return reordered_vecs[kwargs.get("state")] - 340 else: - 341 return reordered_vecs - 342 - 343 def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"): - 344 """Determines the eigenvalue of the GEVP by solving and projecting the correlator - 345 - 346 Parameters - 347 ---------- - 348 state : int - 349 The state one is interested in ordered by energy. The lowest state is zero. - 350 - 351 All other parameters are identical to the ones of Corr.GEVP. - 352 """ - 353 vec = self.GEVP(t0, ts=ts, sort=sort)[state] - 354 return self.projected(vec) - 355 - 356 def Hankel(self, N, periodic=False): - 357 """Constructs an NxN Hankel matrix - 358 - 359 C(t) c(t+1) ... c(t+n-1) - 360 C(t+1) c(t+2) ... c(t+n) - 361 ................. - 362 C(t+(n-1)) c(t+n) ... c(t+2(n-1)) - 363 - 364 Parameters - 365 ---------- - 366 N : int - 367 Dimension of the Hankel matrix - 368 periodic : bool, optional - 369 determines whether the matrix is extended periodically - 370 """ - 371 - 372 if self.N != 1: - 373 raise Exception("Multi-operator Prony not implemented!") - 374 - 375 array = np.empty([N, N], dtype="object") - 376 new_content = [] - 377 for t in range(self.T): - 378 new_content.append(array.copy()) - 379 - 380 def wrap(i): - 381 while i >= self.T: - 382 i -= self.T - 383 return i - 384 - 385 for t in range(self.T): - 386 for i in range(N): - 387 for j in range(N): - 388 if periodic: - 389 new_content[t][i, j] = self.content[wrap(t + i + j)][0] - 390 elif (t + i + j) >= self.T: - 391 new_content[t] = None - 392 else: - 393 new_content[t][i, j] = self.content[t + i + j][0] - 394 - 395 return Corr(new_content) - 396 - 397 def roll(self, dt): - 398 """Periodically shift the correlator by dt timeslices - 399 - 400 Parameters - 401 ---------- - 402 dt : int - 403 number of timeslices - 404 """ - 405 return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) - 406 - 407 def reverse(self): - 408 """Reverse the time ordering of the Corr""" - 409 return Corr(self.content[:: -1]) - 410 - 411 def thin(self, spacing=2, offset=0): - 412 """Thin out a correlator to suppress correlations - 413 - 414 Parameters - 415 ---------- - 416 spacing : int - 417 Keep only every 'spacing'th entry of the correlator - 418 offset : int - 419 Offset the equal spacing - 420 """ - 421 new_content = [] - 422 for t in range(self.T): - 423 if (offset + t) % spacing != 0: - 424 new_content.append(None) - 425 else: - 426 new_content.append(self.content[t]) - 427 return Corr(new_content) - 428 - 429 def correlate(self, partner): - 430 """Correlate the correlator with another correlator or Obs - 431 - 432 Parameters - 433 ---------- - 434 partner : Obs or Corr - 435 partner to correlate the correlator with. - 436 Can either be an Obs which is correlated with all entries of the - 437 correlator or a Corr of same length. - 438 """ - 439 if self.N != 1: - 440 raise Exception("Only one-dimensional correlators can be safely correlated.") - 441 new_content = [] - 442 for x0, t_slice in enumerate(self.content): - 443 if _check_for_none(self, t_slice): - 444 new_content.append(None) - 445 else: - 446 if isinstance(partner, Corr): - 447 if _check_for_none(partner, partner.content[x0]): - 448 new_content.append(None) - 449 else: - 450 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) - 451 elif isinstance(partner, Obs): # Should this include CObs? - 452 new_content.append(np.array([correlate(o, partner) for o in t_slice])) - 453 else: - 454 raise Exception("Can only correlate with an Obs or a Corr.") - 455 - 456 return Corr(new_content) - 457 - 458 def reweight(self, weight, **kwargs): - 459 """Reweight the correlator. - 460 - 461 Parameters - 462 ---------- - 463 weight : Obs - 464 Reweighting factor. An Observable that has to be defined on a superset of the - 465 configurations in obs[i].idl for all i. - 466 all_configs : bool - 467 if True, the reweighted observables are normalized by the average of - 468 the reweighting factor on all configurations in weight.idl and not - 469 on the configurations in obs[i].idl. - 470 """ - 471 if self.N != 1: - 472 raise Exception("Reweighting only implemented for one-dimensional correlators.") - 473 new_content = [] - 474 for t_slice in self.content: - 475 if _check_for_none(self, t_slice): - 476 new_content.append(None) - 477 else: - 478 new_content.append(np.array(reweight(weight, t_slice, **kwargs))) - 479 return Corr(new_content) - 480 - 481 def T_symmetry(self, partner, parity=+1): - 482 """Return the time symmetry average of the correlator and its partner - 483 - 484 Parameters - 485 ---------- - 486 partner : Corr - 487 Time symmetry partner of the Corr - 488 partity : int - 489 Parity quantum number of the correlator, can be +1 or -1 - 490 """ - 491 if self.N != 1: - 492 raise Exception("T_symmetry only implemented for one-dimensional correlators.") - 493 if not isinstance(partner, Corr): - 494 raise Exception("T partner has to be a Corr object.") - 495 if parity not in [+1, -1]: - 496 raise Exception("Parity has to be +1 or -1.") - 497 T_partner = parity * partner.reverse() - 498 - 499 t_slices = [] - 500 test = (self - T_partner) - 501 test.gamma_method() - 502 for x0, t_slice in enumerate(test.content): - 503 if t_slice is not None: - 504 if not t_slice[0].is_zero_within_error(5): - 505 t_slices.append(x0) - 506 if t_slices: - 507 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning) - 508 - 509 return (self + T_partner) / 2 - 510 - 511 def deriv(self, variant="symmetric"): - 512 """Return the first derivative of the correlator with respect to x0. - 513 - 514 Parameters - 515 ---------- - 516 variant : str - 517 decides which definition of the finite differences derivative is used. - 518 Available choice: symmetric, forward, backward, improved, log, default: symmetric - 519 """ - 520 if self.N != 1: - 521 raise Exception("deriv only implemented for one-dimensional correlators.") - 522 if variant == "symmetric": - 523 newcontent = [] - 524 for t in range(1, self.T - 1): - 525 if (self.content[t - 1] is None) or (self.content[t + 1] is None): - 526 newcontent.append(None) - 527 else: - 528 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) - 529 if (all([x is None for x in newcontent])): - 530 raise Exception('Derivative is undefined at all timeslices') - 531 return Corr(newcontent, padding=[1, 1]) - 532 elif variant == "forward": - 533 newcontent = [] - 534 for t in range(self.T - 1): - 535 if (self.content[t] is None) or (self.content[t + 1] is None): - 536 newcontent.append(None) - 537 else: - 538 newcontent.append(self.content[t + 1] - self.content[t]) - 539 if (all([x is None for x in newcontent])): - 540 raise Exception("Derivative is undefined at all timeslices") - 541 return Corr(newcontent, padding=[0, 1]) - 542 elif variant == "backward": - 543 newcontent = [] - 544 for t in range(1, self.T): - 545 if (self.content[t - 1] is None) or (self.content[t] is None): - 546 newcontent.append(None) - 547 else: - 548 newcontent.append(self.content[t] - self.content[t - 1]) - 549 if (all([x is None for x in newcontent])): - 550 raise Exception("Derivative is undefined at all timeslices") - 551 return Corr(newcontent, padding=[1, 0]) - 552 elif variant == "improved": - 553 newcontent = [] - 554 for t in range(2, self.T - 2): - 555 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): - 556 newcontent.append(None) - 557 else: - 558 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2])) - 559 if (all([x is None for x in newcontent])): - 560 raise Exception('Derivative is undefined at all timeslices') - 561 return Corr(newcontent, padding=[2, 2]) - 562 elif variant == 'log': - 563 newcontent = [] - 564 for t in range(self.T): - 565 if (self.content[t] is None) or (self.content[t] <= 0): - 566 newcontent.append(None) - 567 else: - 568 newcontent.append(np.log(self.content[t])) - 569 if (all([x is None for x in newcontent])): - 570 raise Exception("Log is undefined at all timeslices") - 571 logcorr = Corr(newcontent) - 572 return self * logcorr.deriv('symmetric') - 573 else: - 574 raise Exception("Unknown variant.") - 575 - 576 def second_deriv(self, variant="symmetric"): - 577 """Return the second derivative of the correlator with respect to x0. - 578 - 579 Parameters - 580 ---------- - 581 variant : str - 582 decides which definition of the finite differences derivative is used. - 583 Available choice: symmetric, improved, log, default: symmetric - 584 """ - 585 if self.N != 1: - 586 raise Exception("second_deriv only implemented for one-dimensional correlators.") - 587 if variant == "symmetric": - 588 newcontent = [] - 589 for t in range(1, self.T - 1): - 590 if (self.content[t - 1] is None) or (self.content[t + 1] is None): - 591 newcontent.append(None) - 592 else: - 593 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) - 594 if (all([x is None for x in newcontent])): - 595 raise Exception("Derivative is undefined at all timeslices") - 596 return Corr(newcontent, padding=[1, 1]) - 597 elif variant == "improved": - 598 newcontent = [] - 599 for t in range(2, self.T - 2): - 600 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): - 601 newcontent.append(None) - 602 else: - 603 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])) - 604 if (all([x is None for x in newcontent])): - 605 raise Exception("Derivative is undefined at all timeslices") - 606 return Corr(newcontent, padding=[2, 2]) - 607 elif variant == 'log': - 608 newcontent = [] - 609 for t in range(self.T): - 610 if (self.content[t] is None) or (self.content[t] <= 0): - 611 newcontent.append(None) - 612 else: - 613 newcontent.append(np.log(self.content[t])) - 614 if (all([x is None for x in newcontent])): - 615 raise Exception("Log is undefined at all timeslices") - 616 logcorr = Corr(newcontent) - 617 return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2) - 618 else: - 619 raise Exception("Unknown variant.") - 620 - 621 def m_eff(self, variant='log', guess=1.0): - 622 """Returns the effective mass of the correlator as correlator object - 623 - 624 Parameters - 625 ---------- - 626 variant : str - 627 log : uses the standard effective mass log(C(t) / C(t+1)) - 628 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. - 629 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. - 630 See, e.g., arXiv:1205.5380 - 631 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) - 632 logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2 - 633 guess : float - 634 guess for the root finder, only relevant for the root variant - 635 """ - 636 if self.N != 1: - 637 raise Exception('Correlator must be projected before getting m_eff') - 638 if variant == 'log': - 639 newcontent = [] - 640 for t in range(self.T - 1): - 641 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): - 642 newcontent.append(None) - 643 elif self.content[t][0].value / self.content[t + 1][0].value < 0: - 644 newcontent.append(None) - 645 else: - 646 newcontent.append(self.content[t] / self.content[t + 1]) - 647 if (all([x is None for x in newcontent])): - 648 raise Exception('m_eff is undefined at all timeslices') - 649 - 650 return np.log(Corr(newcontent, padding=[0, 1])) - 651 - 652 elif variant == 'logsym': - 653 newcontent = [] - 654 for t in range(1, self.T - 1): - 655 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): - 656 newcontent.append(None) - 657 elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0: - 658 newcontent.append(None) - 659 else: - 660 newcontent.append(self.content[t - 1] / self.content[t + 1]) - 661 if (all([x is None for x in newcontent])): - 662 raise Exception('m_eff is undefined at all timeslices') - 663 - 664 return np.log(Corr(newcontent, padding=[1, 1])) / 2 - 665 - 666 elif variant in ['periodic', 'cosh', 'sinh']: - 667 if variant in ['periodic', 'cosh']: - 668 func = anp.cosh - 669 else: - 670 func = anp.sinh - 671 - 672 def root_function(x, d): - 673 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d - 674 - 675 newcontent = [] - 676 for t in range(self.T - 1): - 677 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): - 678 newcontent.append(None) - 679 # Fill the two timeslices in the middle of the lattice with their predecessors - 680 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: - 681 newcontent.append(newcontent[-1]) - 682 elif self.content[t][0].value / self.content[t + 1][0].value < 0: - 683 newcontent.append(None) - 684 else: - 685 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) - 686 if (all([x is None for x in newcontent])): - 687 raise Exception('m_eff is undefined at all timeslices') - 688 - 689 return Corr(newcontent, padding=[0, 1]) - 690 - 691 elif variant == 'arccosh': - 692 newcontent = [] - 693 for t in range(1, self.T - 1): - 694 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): - 695 newcontent.append(None) - 696 else: - 697 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) - 698 if (all([x is None for x in newcontent])): - 699 raise Exception("m_eff is undefined at all timeslices") - 700 return np.arccosh(Corr(newcontent, padding=[1, 1])) - 701 - 702 else: - 703 raise Exception('Unknown variant.') - 704 - 705 def fit(self, function, fitrange=None, silent=False, **kwargs): - 706 r'''Fits function to the data - 707 - 708 Parameters - 709 ---------- - 710 function : obj - 711 function to fit to the data. See fits.least_squares for details. - 712 fitrange : list - 713 Two element list containing the timeslices on which the fit is supposed to start and stop. - 714 Caution: This range is inclusive as opposed to standard python indexing. - 715 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. - 716 If not specified, self.prange or all timeslices are used. - 717 silent : bool - 718 Decides whether output is printed to the standard output. - 719 ''' - 720 if self.N != 1: - 721 raise Exception("Correlator must be projected before fitting") - 722 - 723 if fitrange is None: - 724 if self.prange: - 725 fitrange = self.prange - 726 else: - 727 fitrange = [0, self.T - 1] - 728 else: - 729 if not isinstance(fitrange, list): - 730 raise Exception("fitrange has to be a list with two elements") - 731 if len(fitrange) != 2: - 732 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") - 733 - 734 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] - 735 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] - 736 result = least_squares(xs, ys, function, silent=silent, **kwargs) - 737 return result - 738 - 739 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): - 740 """ Extract a plateau value from a Corr object - 741 - 742 Parameters - 743 ---------- - 744 plateau_range : list - 745 list with two entries, indicating the first and the last timeslice - 746 of the plateau region. - 747 method : str - 748 method to extract the plateau. - 749 'fit' fits a constant to the plateau region - 750 'avg', 'average' or 'mean' just average over the given timeslices. - 751 auto_gamma : bool - 752 apply gamma_method with default parameters to the Corr. Defaults to None - 753 """ - 754 if not plateau_range: - 755 if self.prange: - 756 plateau_range = self.prange - 757 else: - 758 raise Exception("no plateau range provided") - 759 if self.N != 1: - 760 raise Exception("Correlator must be projected before getting a plateau.") - 761 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): - 762 raise Exception("plateau is undefined at all timeslices in plateaurange.") - 763 if auto_gamma: - 764 self.gamma_method() - 765 if method == "fit": - 766 def const_func(a, t): - 767 return a[0] - 768 return self.fit(const_func, plateau_range)[0] - 769 elif method in ["avg", "average", "mean"]: - 770 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) - 771 return returnvalue - 772 - 773 else: - 774 raise Exception("Unsupported plateau method: " + method) - 775 - 776 def set_prange(self, prange): - 777 """Sets the attribute prange of the Corr object.""" - 778 if not len(prange) == 2: - 779 raise Exception("prange must be a list or array with two values") - 780 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): - 781 raise Exception("Start and end point must be integers") - 782 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): - 783 raise Exception("Start and end point must define a range in the interval 0,T") - 784 - 785 self.prange = prange - 786 return - 787 - 788 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): - 789 """Plots the correlator using the tag of the correlator as label if available. - 790 - 791 Parameters - 792 ---------- - 793 x_range : list - 794 list of two values, determining the range of the x-axis e.g. [4, 8]. - 795 comp : Corr or list of Corr - 796 Correlator or list of correlators which are plotted for comparison. - 797 The tags of these correlators are used as labels if available. - 798 logscale : bool - 799 Sets y-axis to logscale. - 800 plateau : Obs - 801 Plateau value to be visualized in the figure. - 802 fit_res : Fit_result - 803 Fit_result object to be visualized. - 804 ylabel : str - 805 Label for the y-axis. - 806 save : str - 807 path to file in which the figure should be saved. - 808 auto_gamma : bool - 809 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. - 810 hide_sigma : float - 811 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. - 812 references : list - 813 List of floating point values that are displayed as horizontal lines for reference. - 814 title : string - 815 Optional title of the figure. - 816 """ - 817 if self.N != 1: - 818 raise Exception("Correlator must be projected before plotting") - 819 - 820 if auto_gamma: - 821 self.gamma_method() - 822 - 823 if x_range is None: - 824 x_range = [0, self.T - 1] - 825 - 826 fig = plt.figure() - 827 ax1 = fig.add_subplot(111) - 828 - 829 x, y, y_err = self.plottable() - 830 if hide_sigma: - 831 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 - 832 else: - 833 hide_from = None - 834 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) - 835 if logscale: - 836 ax1.set_yscale('log') - 837 else: - 838 if y_range is None: - 839 try: - 840 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)]) - 841 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)]) - 842 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) - 843 except Exception: - 844 pass - 845 else: - 846 ax1.set_ylim(y_range) - 847 if comp: - 848 if isinstance(comp, (Corr, list)): - 849 for corr in comp if isinstance(comp, list) else [comp]: - 850 if auto_gamma: - 851 corr.gamma_method() - 852 x, y, y_err = corr.plottable() - 853 if hide_sigma: - 854 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 - 855 else: - 856 hide_from = None - 857 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) - 858 else: - 859 raise Exception("'comp' must be a correlator or a list of correlators.") - 860 - 861 if plateau: - 862 if isinstance(plateau, Obs): - 863 if auto_gamma: - 864 plateau.gamma_method() - 865 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) - 866 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') - 867 else: - 868 raise Exception("'plateau' must be an Obs") - 869 - 870 if references: - 871 if isinstance(references, list): - 872 for ref in references: - 873 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') - 874 else: - 875 raise Exception("'references' must be a list of floating pint values.") - 876 - 877 if self.prange: - 878 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') - 879 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') - 880 - 881 if fit_res: - 882 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) - 883 ax1.plot(x_samples, - 884 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), - 885 ls='-', marker=',', lw=2) - 886 - 887 ax1.set_xlabel(r'$x_0 / a$') - 888 if ylabel: - 889 ax1.set_ylabel(ylabel) - 890 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) - 891 - 892 handles, labels = ax1.get_legend_handles_labels() - 893 if labels: - 894 ax1.legend() - 895 - 896 if title: - 897 plt.title(title) - 898 - 899 plt.draw() - 900 - 901 if save: - 902 if isinstance(save, str): - 903 fig.savefig(save, bbox_inches='tight') - 904 else: - 905 raise Exception("'save' has to be a string.") - 906 - 907 def spaghetti_plot(self, logscale=True): - 908 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. - 909 - 910 Parameters - 911 ---------- - 912 logscale : bool - 913 Determines whether the scale of the y-axis is logarithmic or standard. - 914 """ - 915 if self.N != 1: - 916 raise Exception("Correlator needs to be projected first.") - 917 - 918 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])) - 919 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] - 920 - 921 for name in mc_names: - 922 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T - 923 - 924 fig = plt.figure() - 925 ax = fig.add_subplot(111) - 926 for dat in data: - 927 ax.plot(x0_vals, dat, ls='-', marker='') - 928 - 929 if logscale is True: - 930 ax.set_yscale('log') - 931 - 932 ax.set_xlabel(r'$x_0 / a$') - 933 plt.title(name) - 934 plt.draw() - 935 - 936 def dump(self, filename, datatype="json.gz", **kwargs): - 937 """Dumps the Corr into a file of chosen type - 938 Parameters - 939 ---------- - 940 filename : str - 941 Name of the file to be saved. - 942 datatype : str - 943 Format of the exported file. Supported formats include - 944 "json.gz" and "pickle" - 945 path : str - 946 specifies a custom path for the file (default '.') - 947 """ - 948 if datatype == "json.gz": - 949 from .input.json import dump_to_json - 950 if 'path' in kwargs: - 951 file_name = kwargs.get('path') + '/' + filename - 952 else: - 953 file_name = filename - 954 dump_to_json(self, file_name) - 955 elif datatype == "pickle": - 956 dump_object(self, filename, **kwargs) - 957 else: - 958 raise Exception("Unknown datatype " + str(datatype)) - 959 - 960 def print(self, print_range=None): - 961 print(self.__repr__(print_range)) - 962 - 963 def __repr__(self, print_range=None): - 964 if print_range is None: - 965 print_range = [0, None] - 966 - 967 content_string = "" - 968 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 - 969 - 970 if self.tag is not None: - 971 content_string += "Description: " + self.tag + "\n" - 972 if self.N != 1: - 973 return content_string - 974 - 975 if print_range[1]: - 976 print_range[1] += 1 - 977 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' - 978 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): - 979 if sub_corr is None: - 980 content_string += str(i + print_range[0]) + '\n' - 981 else: - 982 content_string += str(i + print_range[0]) - 983 for element in sub_corr: - 984 content_string += '\t' + ' ' * int(element >= 0) + str(element) - 985 content_string += '\n' - 986 return content_string - 987 - 988 def __str__(self): - 989 return self.__repr__() - 990 - 991 # We define the basic operations, that can be performed with correlators. - 992 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. - 993 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. - 994 # One could try and tell Obs to check if the y in __mul__ is a Corr and - 995 - 996 def __add__(self, y): - 997 if isinstance(y, Corr): - 998 if ((self.N != y.N) or (self.T != y.T)): - 999 raise Exception("Addition of Corrs with different shape") -1000 newcontent = [] -1001 for t in range(self.T): -1002 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1003 newcontent.append(None) -1004 else: -1005 newcontent.append(self.content[t] + y.content[t]) -1006 return Corr(newcontent) -1007 -1008 elif isinstance(y, (Obs, int, float, CObs)): -1009 newcontent = [] -1010 for t in range(self.T): -1011 if _check_for_none(self, self.content[t]): -1012 newcontent.append(None) -1013 else: -1014 newcontent.append(self.content[t] + y) -1015 return Corr(newcontent, prange=self.prange) -1016 elif isinstance(y, np.ndarray): -1017 if y.shape == (self.T,): -1018 return Corr(list((np.array(self.content).T + y).T)) -1019 else: -1020 raise ValueError("operands could not be broadcast together") -1021 else: -1022 raise TypeError("Corr + wrong type") -1023 -1024 def __mul__(self, y): -1025 if isinstance(y, Corr): -1026 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1027 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1028 newcontent = [] -1029 for t in range(self.T): -1030 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1031 newcontent.append(None) -1032 else: -1033 newcontent.append(self.content[t] * y.content[t]) -1034 return Corr(newcontent) -1035 -1036 elif isinstance(y, (Obs, int, float, CObs)): -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 __truediv__(self, y): -1053 if isinstance(y, Corr): -1054 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1055 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1056 newcontent = [] -1057 for t in range(self.T): -1058 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1059 newcontent.append(None) -1060 else: -1061 newcontent.append(self.content[t] / y.content[t]) -1062 for t in range(self.T): -1063 if _check_for_none(self, newcontent[t]): -1064 continue -1065 if np.isnan(np.sum(newcontent[t]).value): -1066 newcontent[t] = None -1067 -1068 if all([item is None for item in newcontent]): -1069 raise Exception("Division returns completely undefined correlator") -1070 return Corr(newcontent) -1071 -1072 elif isinstance(y, (Obs, CObs)): -1073 if isinstance(y, Obs): -1074 if y.value == 0: -1075 raise Exception('Division by zero will return undefined correlator') -1076 if isinstance(y, CObs): -1077 if y.is_zero(): -1078 raise Exception('Division by zero will return undefined correlator') -1079 -1080 newcontent = [] -1081 for t in range(self.T): -1082 if _check_for_none(self, self.content[t]): -1083 newcontent.append(None) -1084 else: -1085 newcontent.append(self.content[t] / y) -1086 return Corr(newcontent, prange=self.prange) -1087 -1088 elif isinstance(y, (int, float)): -1089 if y == 0: -1090 raise Exception('Division by zero will return undefined correlator') -1091 newcontent = [] -1092 for t in range(self.T): -1093 if _check_for_none(self, self.content[t]): -1094 newcontent.append(None) -1095 else: -1096 newcontent.append(self.content[t] / y) -1097 return Corr(newcontent, prange=self.prange) -1098 elif isinstance(y, np.ndarray): -1099 if y.shape == (self.T,): -1100 return Corr(list((np.array(self.content).T / y).T)) -1101 else: -1102 raise ValueError("operands could not be broadcast together") -1103 else: -1104 raise TypeError('Corr / wrong type') -1105 -1106 def __neg__(self): -1107 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] -1108 return Corr(newcontent, prange=self.prange) -1109 -1110 def __sub__(self, y): -1111 return self + (-y) -1112 -1113 def __pow__(self, y): -1114 if isinstance(y, (Obs, int, float, CObs)): -1115 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] -1116 return Corr(newcontent, prange=self.prange) -1117 else: -1118 raise TypeError('Type of exponent not supported') -1119 -1120 def __abs__(self): -1121 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] -1122 return Corr(newcontent, prange=self.prange) -1123 -1124 # The numpy functions: -1125 def sqrt(self): -1126 return self ** 0.5 -1127 -1128 def log(self): -1129 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] -1130 return Corr(newcontent, prange=self.prange) -1131 -1132 def exp(self): -1133 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] -1134 return Corr(newcontent, prange=self.prange) -1135 -1136 def _apply_func_to_corr(self, func): -1137 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] -1138 for t in range(self.T): -1139 if _check_for_none(self, newcontent[t]): -1140 continue -1141 if np.isnan(np.sum(newcontent[t]).value): -1142 newcontent[t] = None -1143 if all([item is None for item in newcontent]): -1144 raise Exception('Operation returns undefined correlator') -1145 return Corr(newcontent) -1146 -1147 def sin(self): -1148 return self._apply_func_to_corr(np.sin) -1149 -1150 def cos(self): -1151 return self._apply_func_to_corr(np.cos) -1152 -1153 def tan(self): -1154 return self._apply_func_to_corr(np.tan) -1155 -1156 def sinh(self): -1157 return self._apply_func_to_corr(np.sinh) -1158 -1159 def cosh(self): -1160 return self._apply_func_to_corr(np.cosh) -1161 -1162 def tanh(self): -1163 return self._apply_func_to_corr(np.tanh) -1164 -1165 def arcsin(self): -1166 return self._apply_func_to_corr(np.arcsin) -1167 -1168 def arccos(self): -1169 return self._apply_func_to_corr(np.arccos) -1170 -1171 def arctan(self): -1172 return self._apply_func_to_corr(np.arctan) -1173 -1174 def arcsinh(self): -1175 return self._apply_func_to_corr(np.arcsinh) -1176 -1177 def arccosh(self): -1178 return self._apply_func_to_corr(np.arccosh) -1179 -1180 def arctanh(self): -1181 return self._apply_func_to_corr(np.arctanh) -1182 -1183 # Right hand side operations (require tweak in main module to work) -1184 def __radd__(self, y): -1185 return self + y -1186 -1187 def __rsub__(self, y): -1188 return -self + y -1189 -1190 def __rmul__(self, y): -1191 return self * y -1192 -1193 def __rtruediv__(self, y): -1194 return (self / y) ** (-1) -1195 -1196 @property -1197 def real(self): -1198 def return_real(obs_OR_cobs): -1199 if isinstance(obs_OR_cobs, CObs): -1200 return obs_OR_cobs.real -1201 else: -1202 return obs_OR_cobs -1203 -1204 return self._apply_func_to_corr(return_real) -1205 -1206 @property -1207 def imag(self): -1208 def return_imag(obs_OR_cobs): -1209 if isinstance(obs_OR_cobs, CObs): -1210 return obs_OR_cobs.imag -1211 else: -1212 return obs_OR_cobs * 0 # So it stays the right type -1213 -1214 return self._apply_func_to_corr(return_imag) -1215 -1216 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): -1217 r''' Project large correlation matrix to lowest states -1218 -1219 This method can be used to reduce the size of an (N x N) correlation matrix -1220 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise -1221 is still small. -1222 -1223 Parameters -1224 ---------- -1225 Ntrunc: int -1226 Rank of the target matrix. -1227 tproj: int -1228 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. -1229 The default value is 3. -1230 t0proj: int -1231 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly -1232 discouraged for O(a) improved theories, since the correctness of the procedure -1233 cannot be granted in this case. The default value is 2. -1234 basematrix : Corr -1235 Correlation matrix that is used to determine the eigenvectors of the -1236 lowest states based on a GEVP. basematrix is taken to be the Corr itself if -1237 is is not specified. -1238 -1239 Notes -1240 ----- -1241 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving -1242 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}$ -1243 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the -1244 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via -1245 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large -1246 correlation matrix and to remove some noise that is added by irrelevant operators. -1247 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated -1248 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. -1249 ''' -1250 -1251 if self.N == 1: -1252 raise Exception('Method cannot be applied to one-dimensional correlators.') -1253 if basematrix is None: -1254 basematrix = self -1255 if Ntrunc >= basematrix.N: -1256 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) -1257 if basematrix.N != self.N: -1258 raise Exception('basematrix and targetmatrix have to be of the same size.') -1259 -1260 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] -1261 -1262 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) -1263 rmat = [] -1264 for t in range(basematrix.T): -1265 for i in range(Ntrunc): -1266 for j in range(Ntrunc): -1267 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] -1268 rmat.append(np.copy(tmpmat)) -1269 -1270 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] -1271 return Corr(newcontent) -1272 + 309 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) + 310 np.linalg.cholesky(G0) # Check if matrix G0 is positive-semidefinite. + 311 + 312 if sort is None: + 313 if (ts is None): + 314 raise Exception("ts is required if sort=None.") + 315 if (self.content[t0] is None) or (self.content[ts] is None): + 316 raise Exception("Corr not defined at t0/ts.") + 317 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts]) + 318 reordered_vecs = _GEVP_solver(Gt, G0) + 319 + 320 elif sort in ["Eigenvalue", "Eigenvector"]: + 321 if sort == "Eigenvalue" and ts is not None: + 322 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) + 323 all_vecs = [None] * (t0 + 1) + 324 for t in range(t0 + 1, self.T): + 325 try: + 326 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t]) + 327 all_vecs.append(_GEVP_solver(Gt, G0)) + 328 except Exception: + 329 all_vecs.append(None) + 330 if sort == "Eigenvector": + 331 if (ts is None): + 332 raise Exception("ts is required for the Eigenvector sorting method.") + 333 all_vecs = _sort_vectors(all_vecs, ts) + 334 + 335 reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] + 336 else: + 337 raise Exception("Unkown value for 'sort'.") + 338 + 339 if "state" in kwargs: + 340 return reordered_vecs[kwargs.get("state")] + 341 else: + 342 return reordered_vecs + 343 + 344 def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"): + 345 """Determines the eigenvalue of the GEVP by solving and projecting the correlator + 346 + 347 Parameters + 348 ---------- + 349 state : int + 350 The state one is interested in ordered by energy. The lowest state is zero. + 351 + 352 All other parameters are identical to the ones of Corr.GEVP. + 353 """ + 354 vec = self.GEVP(t0, ts=ts, sort=sort)[state] + 355 return self.projected(vec) + 356 + 357 def Hankel(self, N, periodic=False): + 358 """Constructs an NxN Hankel matrix + 359 + 360 C(t) c(t+1) ... c(t+n-1) + 361 C(t+1) c(t+2) ... c(t+n) + 362 ................. + 363 C(t+(n-1)) c(t+n) ... c(t+2(n-1)) + 364 + 365 Parameters + 366 ---------- + 367 N : int + 368 Dimension of the Hankel matrix + 369 periodic : bool, optional + 370 determines whether the matrix is extended periodically + 371 """ + 372 + 373 if self.N != 1: + 374 raise Exception("Multi-operator Prony not implemented!") + 375 + 376 array = np.empty([N, N], dtype="object") + 377 new_content = [] + 378 for t in range(self.T): + 379 new_content.append(array.copy()) + 380 + 381 def wrap(i): + 382 while i >= self.T: + 383 i -= self.T + 384 return i + 385 + 386 for t in range(self.T): + 387 for i in range(N): + 388 for j in range(N): + 389 if periodic: + 390 new_content[t][i, j] = self.content[wrap(t + i + j)][0] + 391 elif (t + i + j) >= self.T: + 392 new_content[t] = None + 393 else: + 394 new_content[t][i, j] = self.content[t + i + j][0] + 395 + 396 return Corr(new_content) + 397 + 398 def roll(self, dt): + 399 """Periodically shift the correlator by dt timeslices + 400 + 401 Parameters + 402 ---------- + 403 dt : int + 404 number of timeslices + 405 """ + 406 return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) + 407 + 408 def reverse(self): + 409 """Reverse the time ordering of the Corr""" + 410 return Corr(self.content[:: -1]) + 411 + 412 def thin(self, spacing=2, offset=0): + 413 """Thin out a correlator to suppress correlations + 414 + 415 Parameters + 416 ---------- + 417 spacing : int + 418 Keep only every 'spacing'th entry of the correlator + 419 offset : int + 420 Offset the equal spacing + 421 """ + 422 new_content = [] + 423 for t in range(self.T): + 424 if (offset + t) % spacing != 0: + 425 new_content.append(None) + 426 else: + 427 new_content.append(self.content[t]) + 428 return Corr(new_content) + 429 + 430 def correlate(self, partner): + 431 """Correlate the correlator with another correlator or Obs + 432 + 433 Parameters + 434 ---------- + 435 partner : Obs or Corr + 436 partner to correlate the correlator with. + 437 Can either be an Obs which is correlated with all entries of the + 438 correlator or a Corr of same length. + 439 """ + 440 if self.N != 1: + 441 raise Exception("Only one-dimensional correlators can be safely correlated.") + 442 new_content = [] + 443 for x0, t_slice in enumerate(self.content): + 444 if _check_for_none(self, t_slice): + 445 new_content.append(None) + 446 else: + 447 if isinstance(partner, Corr): + 448 if _check_for_none(partner, partner.content[x0]): + 449 new_content.append(None) + 450 else: + 451 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) + 452 elif isinstance(partner, Obs): # Should this include CObs? + 453 new_content.append(np.array([correlate(o, partner) for o in t_slice])) + 454 else: + 455 raise Exception("Can only correlate with an Obs or a Corr.") + 456 + 457 return Corr(new_content) + 458 + 459 def reweight(self, weight, **kwargs): + 460 """Reweight the correlator. + 461 + 462 Parameters + 463 ---------- + 464 weight : Obs + 465 Reweighting factor. An Observable that has to be defined on a superset of the + 466 configurations in obs[i].idl for all i. + 467 all_configs : bool + 468 if True, the reweighted observables are normalized by the average of + 469 the reweighting factor on all configurations in weight.idl and not + 470 on the configurations in obs[i].idl. + 471 """ + 472 if self.N != 1: + 473 raise Exception("Reweighting only implemented for one-dimensional correlators.") + 474 new_content = [] + 475 for t_slice in self.content: + 476 if _check_for_none(self, t_slice): + 477 new_content.append(None) + 478 else: + 479 new_content.append(np.array(reweight(weight, t_slice, **kwargs))) + 480 return Corr(new_content) + 481 + 482 def T_symmetry(self, partner, parity=+1): + 483 """Return the time symmetry average of the correlator and its partner + 484 + 485 Parameters + 486 ---------- + 487 partner : Corr + 488 Time symmetry partner of the Corr + 489 partity : int + 490 Parity quantum number of the correlator, can be +1 or -1 + 491 """ + 492 if self.N != 1: + 493 raise Exception("T_symmetry only implemented for one-dimensional correlators.") + 494 if not isinstance(partner, Corr): + 495 raise Exception("T partner has to be a Corr object.") + 496 if parity not in [+1, -1]: + 497 raise Exception("Parity has to be +1 or -1.") + 498 T_partner = parity * partner.reverse() + 499 + 500 t_slices = [] + 501 test = (self - T_partner) + 502 test.gamma_method() + 503 for x0, t_slice in enumerate(test.content): + 504 if t_slice is not None: + 505 if not t_slice[0].is_zero_within_error(5): + 506 t_slices.append(x0) + 507 if t_slices: + 508 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning) + 509 + 510 return (self + T_partner) / 2 + 511 + 512 def deriv(self, variant="symmetric"): + 513 """Return the first derivative of the correlator with respect to x0. + 514 + 515 Parameters + 516 ---------- + 517 variant : str + 518 decides which definition of the finite differences derivative is used. + 519 Available choice: symmetric, forward, backward, improved, log, default: symmetric + 520 """ + 521 if self.N != 1: + 522 raise Exception("deriv only implemented for one-dimensional correlators.") + 523 if variant == "symmetric": + 524 newcontent = [] + 525 for t in range(1, self.T - 1): + 526 if (self.content[t - 1] is None) or (self.content[t + 1] is None): + 527 newcontent.append(None) + 528 else: + 529 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) + 530 if (all([x is None for x in newcontent])): + 531 raise Exception('Derivative is undefined at all timeslices') + 532 return Corr(newcontent, padding=[1, 1]) + 533 elif variant == "forward": + 534 newcontent = [] + 535 for t in range(self.T - 1): + 536 if (self.content[t] is None) or (self.content[t + 1] is None): + 537 newcontent.append(None) + 538 else: + 539 newcontent.append(self.content[t + 1] - self.content[t]) + 540 if (all([x is None for x in newcontent])): + 541 raise Exception("Derivative is undefined at all timeslices") + 542 return Corr(newcontent, padding=[0, 1]) + 543 elif variant == "backward": + 544 newcontent = [] + 545 for t in range(1, self.T): + 546 if (self.content[t - 1] is None) or (self.content[t] is None): + 547 newcontent.append(None) + 548 else: + 549 newcontent.append(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, 0]) + 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 + 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] - 8 * self.content[t - 1] + 8 * 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 elif variant == 'log': + 564 newcontent = [] + 565 for t in range(self.T): + 566 if (self.content[t] is None) or (self.content[t] <= 0): + 567 newcontent.append(None) + 568 else: + 569 newcontent.append(np.log(self.content[t])) + 570 if (all([x is None for x in newcontent])): + 571 raise Exception("Log is undefined at all timeslices") + 572 logcorr = Corr(newcontent) + 573 return self * logcorr.deriv('symmetric') + 574 else: + 575 raise Exception("Unknown variant.") + 576 + 577 def second_deriv(self, variant="symmetric"): + 578 """Return the second derivative of the correlator with respect to x0. + 579 + 580 Parameters + 581 ---------- + 582 variant : str + 583 decides which definition of the finite differences derivative is used. + 584 Available choice: symmetric, improved, log, default: symmetric + 585 """ + 586 if self.N != 1: + 587 raise Exception("second_deriv only implemented for one-dimensional correlators.") + 588 if variant == "symmetric": + 589 newcontent = [] + 590 for t in range(1, self.T - 1): + 591 if (self.content[t - 1] is None) or (self.content[t + 1] is None): + 592 newcontent.append(None) + 593 else: + 594 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) + 595 if (all([x is None for x in newcontent])): + 596 raise Exception("Derivative is undefined at all timeslices") + 597 return Corr(newcontent, padding=[1, 1]) + 598 elif variant == "improved": + 599 newcontent = [] + 600 for t in range(2, self.T - 2): + 601 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): + 602 newcontent.append(None) + 603 else: + 604 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])) + 605 if (all([x is None for x in newcontent])): + 606 raise Exception("Derivative is undefined at all timeslices") + 607 return Corr(newcontent, padding=[2, 2]) + 608 elif variant == 'log': + 609 newcontent = [] + 610 for t in range(self.T): + 611 if (self.content[t] is None) or (self.content[t] <= 0): + 612 newcontent.append(None) + 613 else: + 614 newcontent.append(np.log(self.content[t])) + 615 if (all([x is None for x in newcontent])): + 616 raise Exception("Log is undefined at all timeslices") + 617 logcorr = Corr(newcontent) + 618 return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2) + 619 else: + 620 raise Exception("Unknown variant.") + 621 + 622 def m_eff(self, variant='log', guess=1.0): + 623 """Returns the effective mass of the correlator as correlator object + 624 + 625 Parameters + 626 ---------- + 627 variant : str + 628 log : uses the standard effective mass log(C(t) / C(t+1)) + 629 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. + 630 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. + 631 See, e.g., arXiv:1205.5380 + 632 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) + 633 logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2 + 634 guess : float + 635 guess for the root finder, only relevant for the root variant + 636 """ + 637 if self.N != 1: + 638 raise Exception('Correlator must be projected before getting m_eff') + 639 if variant == 'log': + 640 newcontent = [] + 641 for t in range(self.T - 1): + 642 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): + 643 newcontent.append(None) + 644 elif self.content[t][0].value / self.content[t + 1][0].value < 0: + 645 newcontent.append(None) + 646 else: + 647 newcontent.append(self.content[t] / self.content[t + 1]) + 648 if (all([x is None for x in newcontent])): + 649 raise Exception('m_eff is undefined at all timeslices') + 650 + 651 return np.log(Corr(newcontent, padding=[0, 1])) + 652 + 653 elif variant == 'logsym': + 654 newcontent = [] + 655 for t in range(1, self.T - 1): + 656 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): + 657 newcontent.append(None) + 658 elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0: + 659 newcontent.append(None) + 660 else: + 661 newcontent.append(self.content[t - 1] / self.content[t + 1]) + 662 if (all([x is None for x in newcontent])): + 663 raise Exception('m_eff is undefined at all timeslices') + 664 + 665 return np.log(Corr(newcontent, padding=[1, 1])) / 2 + 666 + 667 elif variant in ['periodic', 'cosh', 'sinh']: + 668 if variant in ['periodic', 'cosh']: + 669 func = anp.cosh + 670 else: + 671 func = anp.sinh + 672 + 673 def root_function(x, d): + 674 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d + 675 + 676 newcontent = [] + 677 for t in range(self.T - 1): + 678 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): + 679 newcontent.append(None) + 680 # Fill the two timeslices in the middle of the lattice with their predecessors + 681 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: + 682 newcontent.append(newcontent[-1]) + 683 elif self.content[t][0].value / self.content[t + 1][0].value < 0: + 684 newcontent.append(None) + 685 else: + 686 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) + 687 if (all([x is None for x in newcontent])): + 688 raise Exception('m_eff is undefined at all timeslices') + 689 + 690 return Corr(newcontent, padding=[0, 1]) + 691 + 692 elif variant == 'arccosh': + 693 newcontent = [] + 694 for t in range(1, self.T - 1): + 695 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): + 696 newcontent.append(None) + 697 else: + 698 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) + 699 if (all([x is None for x in newcontent])): + 700 raise Exception("m_eff is undefined at all timeslices") + 701 return np.arccosh(Corr(newcontent, padding=[1, 1])) + 702 + 703 else: + 704 raise Exception('Unknown variant.') + 705 + 706 def fit(self, function, fitrange=None, silent=False, **kwargs): + 707 r'''Fits function to the data + 708 + 709 Parameters + 710 ---------- + 711 function : obj + 712 function to fit to the data. See fits.least_squares for details. + 713 fitrange : list + 714 Two element list containing the timeslices on which the fit is supposed to start and stop. + 715 Caution: This range is inclusive as opposed to standard python indexing. + 716 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. + 717 If not specified, self.prange or all timeslices are used. + 718 silent : bool + 719 Decides whether output is printed to the standard output. + 720 ''' + 721 if self.N != 1: + 722 raise Exception("Correlator must be projected before fitting") + 723 + 724 if fitrange is None: + 725 if self.prange: + 726 fitrange = self.prange + 727 else: + 728 fitrange = [0, self.T - 1] + 729 else: + 730 if not isinstance(fitrange, list): + 731 raise Exception("fitrange has to be a list with two elements") + 732 if len(fitrange) != 2: + 733 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") + 734 + 735 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] + 736 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] + 737 result = least_squares(xs, ys, function, silent=silent, **kwargs) + 738 return result + 739 + 740 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): + 741 """ Extract a plateau value from a Corr object + 742 + 743 Parameters + 744 ---------- + 745 plateau_range : list + 746 list with two entries, indicating the first and the last timeslice + 747 of the plateau region. + 748 method : str + 749 method to extract the plateau. + 750 'fit' fits a constant to the plateau region + 751 'avg', 'average' or 'mean' just average over the given timeslices. + 752 auto_gamma : bool + 753 apply gamma_method with default parameters to the Corr. Defaults to None + 754 """ + 755 if not plateau_range: + 756 if self.prange: + 757 plateau_range = self.prange + 758 else: + 759 raise Exception("no plateau range provided") + 760 if self.N != 1: + 761 raise Exception("Correlator must be projected before getting a plateau.") + 762 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): + 763 raise Exception("plateau is undefined at all timeslices in plateaurange.") + 764 if auto_gamma: + 765 self.gamma_method() + 766 if method == "fit": + 767 def const_func(a, t): + 768 return a[0] + 769 return self.fit(const_func, plateau_range)[0] + 770 elif method in ["avg", "average", "mean"]: + 771 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) + 772 return returnvalue + 773 + 774 else: + 775 raise Exception("Unsupported plateau method: " + method) + 776 + 777 def set_prange(self, prange): + 778 """Sets the attribute prange of the Corr object.""" + 779 if not len(prange) == 2: + 780 raise Exception("prange must be a list or array with two values") + 781 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): + 782 raise Exception("Start and end point must be integers") + 783 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): + 784 raise Exception("Start and end point must define a range in the interval 0,T") + 785 + 786 self.prange = prange + 787 return + 788 + 789 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): + 790 """Plots the correlator using the tag of the correlator as label if available. + 791 + 792 Parameters + 793 ---------- + 794 x_range : list + 795 list of two values, determining the range of the x-axis e.g. [4, 8]. + 796 comp : Corr or list of Corr + 797 Correlator or list of correlators which are plotted for comparison. + 798 The tags of these correlators are used as labels if available. + 799 logscale : bool + 800 Sets y-axis to logscale. + 801 plateau : Obs + 802 Plateau value to be visualized in the figure. + 803 fit_res : Fit_result + 804 Fit_result object to be visualized. + 805 ylabel : str + 806 Label for the y-axis. + 807 save : str + 808 path to file in which the figure should be saved. + 809 auto_gamma : bool + 810 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. + 811 hide_sigma : float + 812 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. + 813 references : list + 814 List of floating point values that are displayed as horizontal lines for reference. + 815 title : string + 816 Optional title of the figure. + 817 """ + 818 if self.N != 1: + 819 raise Exception("Correlator must be projected before plotting") + 820 + 821 if auto_gamma: + 822 self.gamma_method() + 823 + 824 if x_range is None: + 825 x_range = [0, self.T - 1] + 826 + 827 fig = plt.figure() + 828 ax1 = fig.add_subplot(111) + 829 + 830 x, y, y_err = self.plottable() + 831 if hide_sigma: + 832 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 + 833 else: + 834 hide_from = None + 835 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) + 836 if logscale: + 837 ax1.set_yscale('log') + 838 else: + 839 if y_range is None: + 840 try: + 841 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)]) + 842 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)]) + 843 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) + 844 except Exception: + 845 pass + 846 else: + 847 ax1.set_ylim(y_range) + 848 if comp: + 849 if isinstance(comp, (Corr, list)): + 850 for corr in comp if isinstance(comp, list) else [comp]: + 851 if auto_gamma: + 852 corr.gamma_method() + 853 x, y, y_err = corr.plottable() + 854 if hide_sigma: + 855 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 + 856 else: + 857 hide_from = None + 858 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) + 859 else: + 860 raise Exception("'comp' must be a correlator or a list of correlators.") + 861 + 862 if plateau: + 863 if isinstance(plateau, Obs): + 864 if auto_gamma: + 865 plateau.gamma_method() + 866 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) + 867 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') + 868 else: + 869 raise Exception("'plateau' must be an Obs") + 870 + 871 if references: + 872 if isinstance(references, list): + 873 for ref in references: + 874 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') + 875 else: + 876 raise Exception("'references' must be a list of floating pint values.") + 877 + 878 if self.prange: + 879 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') + 880 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') + 881 + 882 if fit_res: + 883 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) + 884 ax1.plot(x_samples, + 885 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), + 886 ls='-', marker=',', lw=2) + 887 + 888 ax1.set_xlabel(r'$x_0 / a$') + 889 if ylabel: + 890 ax1.set_ylabel(ylabel) + 891 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) + 892 + 893 handles, labels = ax1.get_legend_handles_labels() + 894 if labels: + 895 ax1.legend() + 896 + 897 if title: + 898 plt.title(title) + 899 + 900 plt.draw() + 901 + 902 if save: + 903 if isinstance(save, str): + 904 fig.savefig(save, bbox_inches='tight') + 905 else: + 906 raise Exception("'save' has to be a string.") + 907 + 908 def spaghetti_plot(self, logscale=True): + 909 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. + 910 + 911 Parameters + 912 ---------- + 913 logscale : bool + 914 Determines whether the scale of the y-axis is logarithmic or standard. + 915 """ + 916 if self.N != 1: + 917 raise Exception("Correlator needs to be projected first.") + 918 + 919 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])) + 920 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] + 921 + 922 for name in mc_names: + 923 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T + 924 + 925 fig = plt.figure() + 926 ax = fig.add_subplot(111) + 927 for dat in data: + 928 ax.plot(x0_vals, dat, ls='-', marker='') + 929 + 930 if logscale is True: + 931 ax.set_yscale('log') + 932 + 933 ax.set_xlabel(r'$x_0 / a$') + 934 plt.title(name) + 935 plt.draw() + 936 + 937 def dump(self, filename, datatype="json.gz", **kwargs): + 938 """Dumps the Corr into a file of chosen type + 939 Parameters + 940 ---------- + 941 filename : str + 942 Name of the file to be saved. + 943 datatype : str + 944 Format of the exported file. Supported formats include + 945 "json.gz" and "pickle" + 946 path : str + 947 specifies a custom path for the file (default '.') + 948 """ + 949 if datatype == "json.gz": + 950 from .input.json import dump_to_json + 951 if 'path' in kwargs: + 952 file_name = kwargs.get('path') + '/' + filename + 953 else: + 954 file_name = filename + 955 dump_to_json(self, file_name) + 956 elif datatype == "pickle": + 957 dump_object(self, filename, **kwargs) + 958 else: + 959 raise Exception("Unknown datatype " + str(datatype)) + 960 + 961 def print(self, print_range=None): + 962 print(self.__repr__(print_range)) + 963 + 964 def __repr__(self, print_range=None): + 965 if print_range is None: + 966 print_range = [0, None] + 967 + 968 content_string = "" + 969 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 + 970 + 971 if self.tag is not None: + 972 content_string += "Description: " + self.tag + "\n" + 973 if self.N != 1: + 974 return content_string + 975 + 976 if print_range[1]: + 977 print_range[1] += 1 + 978 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' + 979 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): + 980 if sub_corr is None: + 981 content_string += str(i + print_range[0]) + '\n' + 982 else: + 983 content_string += str(i + print_range[0]) + 984 for element in sub_corr: + 985 content_string += '\t' + ' ' * int(element >= 0) + str(element) + 986 content_string += '\n' + 987 return content_string + 988 + 989 def __str__(self): + 990 return self.__repr__() + 991 + 992 # We define the basic operations, that can be performed with correlators. + 993 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. + 994 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. + 995 # One could try and tell Obs to check if the y in __mul__ is a Corr and + 996 + 997 def __add__(self, y): + 998 if isinstance(y, Corr): + 999 if ((self.N != y.N) or (self.T != y.T)): +1000 raise Exception("Addition of Corrs with different shape") +1001 newcontent = [] +1002 for t in range(self.T): +1003 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1004 newcontent.append(None) +1005 else: +1006 newcontent.append(self.content[t] + y.content[t]) +1007 return Corr(newcontent) +1008 +1009 elif isinstance(y, (Obs, int, float, CObs)): +1010 newcontent = [] +1011 for t in range(self.T): +1012 if _check_for_none(self, self.content[t]): +1013 newcontent.append(None) +1014 else: +1015 newcontent.append(self.content[t] + y) +1016 return Corr(newcontent, prange=self.prange) +1017 elif isinstance(y, np.ndarray): +1018 if y.shape == (self.T,): +1019 return Corr(list((np.array(self.content).T + y).T)) +1020 else: +1021 raise ValueError("operands could not be broadcast together") +1022 else: +1023 raise TypeError("Corr + wrong type") +1024 +1025 def __mul__(self, y): +1026 if isinstance(y, Corr): +1027 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1028 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1029 newcontent = [] +1030 for t in range(self.T): +1031 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1032 newcontent.append(None) +1033 else: +1034 newcontent.append(self.content[t] * y.content[t]) +1035 return Corr(newcontent) +1036 +1037 elif isinstance(y, (Obs, int, float, CObs)): +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 __truediv__(self, y): +1054 if isinstance(y, Corr): +1055 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1056 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1057 newcontent = [] +1058 for t in range(self.T): +1059 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1060 newcontent.append(None) +1061 else: +1062 newcontent.append(self.content[t] / y.content[t]) +1063 for t in range(self.T): +1064 if _check_for_none(self, newcontent[t]): +1065 continue +1066 if np.isnan(np.sum(newcontent[t]).value): +1067 newcontent[t] = None +1068 +1069 if all([item is None for item in newcontent]): +1070 raise Exception("Division returns completely undefined correlator") +1071 return Corr(newcontent) +1072 +1073 elif isinstance(y, (Obs, CObs)): +1074 if isinstance(y, Obs): +1075 if y.value == 0: +1076 raise Exception('Division by zero will return undefined correlator') +1077 if isinstance(y, CObs): +1078 if y.is_zero(): +1079 raise Exception('Division by zero will return undefined correlator') +1080 +1081 newcontent = [] +1082 for t in range(self.T): +1083 if _check_for_none(self, self.content[t]): +1084 newcontent.append(None) +1085 else: +1086 newcontent.append(self.content[t] / y) +1087 return Corr(newcontent, prange=self.prange) +1088 +1089 elif isinstance(y, (int, float)): +1090 if y == 0: +1091 raise Exception('Division by zero will return undefined correlator') +1092 newcontent = [] +1093 for t in range(self.T): +1094 if _check_for_none(self, self.content[t]): +1095 newcontent.append(None) +1096 else: +1097 newcontent.append(self.content[t] / y) +1098 return Corr(newcontent, prange=self.prange) +1099 elif isinstance(y, np.ndarray): +1100 if y.shape == (self.T,): +1101 return Corr(list((np.array(self.content).T / y).T)) +1102 else: +1103 raise ValueError("operands could not be broadcast together") +1104 else: +1105 raise TypeError('Corr / wrong type') +1106 +1107 def __neg__(self): +1108 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] +1109 return Corr(newcontent, prange=self.prange) +1110 +1111 def __sub__(self, y): +1112 return self + (-y) +1113 +1114 def __pow__(self, y): +1115 if isinstance(y, (Obs, int, float, CObs)): +1116 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] +1117 return Corr(newcontent, prange=self.prange) +1118 else: +1119 raise TypeError('Type of exponent not supported') +1120 +1121 def __abs__(self): +1122 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] +1123 return Corr(newcontent, prange=self.prange) +1124 +1125 # The numpy functions: +1126 def sqrt(self): +1127 return self ** 0.5 +1128 +1129 def log(self): +1130 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] +1131 return Corr(newcontent, prange=self.prange) +1132 +1133 def exp(self): +1134 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] +1135 return Corr(newcontent, prange=self.prange) +1136 +1137 def _apply_func_to_corr(self, func): +1138 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] +1139 for t in range(self.T): +1140 if _check_for_none(self, newcontent[t]): +1141 continue +1142 if np.isnan(np.sum(newcontent[t]).value): +1143 newcontent[t] = None +1144 if all([item is None for item in newcontent]): +1145 raise Exception('Operation returns undefined correlator') +1146 return Corr(newcontent) +1147 +1148 def sin(self): +1149 return self._apply_func_to_corr(np.sin) +1150 +1151 def cos(self): +1152 return self._apply_func_to_corr(np.cos) +1153 +1154 def tan(self): +1155 return self._apply_func_to_corr(np.tan) +1156 +1157 def sinh(self): +1158 return self._apply_func_to_corr(np.sinh) +1159 +1160 def cosh(self): +1161 return self._apply_func_to_corr(np.cosh) +1162 +1163 def tanh(self): +1164 return self._apply_func_to_corr(np.tanh) +1165 +1166 def arcsin(self): +1167 return self._apply_func_to_corr(np.arcsin) +1168 +1169 def arccos(self): +1170 return self._apply_func_to_corr(np.arccos) +1171 +1172 def arctan(self): +1173 return self._apply_func_to_corr(np.arctan) +1174 +1175 def arcsinh(self): +1176 return self._apply_func_to_corr(np.arcsinh) +1177 +1178 def arccosh(self): +1179 return self._apply_func_to_corr(np.arccosh) +1180 +1181 def arctanh(self): +1182 return self._apply_func_to_corr(np.arctanh) +1183 +1184 # Right hand side operations (require tweak in main module to work) +1185 def __radd__(self, y): +1186 return self + y +1187 +1188 def __rsub__(self, y): +1189 return -self + y +1190 +1191 def __rmul__(self, y): +1192 return self * y +1193 +1194 def __rtruediv__(self, y): +1195 return (self / y) ** (-1) +1196 +1197 @property +1198 def real(self): +1199 def return_real(obs_OR_cobs): +1200 if isinstance(obs_OR_cobs, CObs): +1201 return obs_OR_cobs.real +1202 else: +1203 return obs_OR_cobs +1204 +1205 return self._apply_func_to_corr(return_real) +1206 +1207 @property +1208 def imag(self): +1209 def return_imag(obs_OR_cobs): +1210 if isinstance(obs_OR_cobs, CObs): +1211 return obs_OR_cobs.imag +1212 else: +1213 return obs_OR_cobs * 0 # So it stays the right type +1214 +1215 return self._apply_func_to_corr(return_imag) +1216 +1217 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): +1218 r''' Project large correlation matrix to lowest states +1219 +1220 This method can be used to reduce the size of an (N x N) correlation matrix +1221 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise +1222 is still small. +1223 +1224 Parameters +1225 ---------- +1226 Ntrunc: int +1227 Rank of the target matrix. +1228 tproj: int +1229 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. +1230 The default value is 3. +1231 t0proj: int +1232 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly +1233 discouraged for O(a) improved theories, since the correctness of the procedure +1234 cannot be granted in this case. The default value is 2. +1235 basematrix : Corr +1236 Correlation matrix that is used to determine the eigenvectors of the +1237 lowest states based on a GEVP. basematrix is taken to be the Corr itself if +1238 is is not specified. +1239 +1240 Notes +1241 ----- +1242 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving +1243 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}$ +1244 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the +1245 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via +1246 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large +1247 correlation matrix and to remove some noise that is added by irrelevant operators. +1248 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated +1249 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. +1250 ''' +1251 +1252 if self.N == 1: +1253 raise Exception('Method cannot be applied to one-dimensional correlators.') +1254 if basematrix is None: +1255 basematrix = self +1256 if Ntrunc >= basematrix.N: +1257 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) +1258 if basematrix.N != self.N: +1259 raise Exception('basematrix and targetmatrix have to be of the same size.') +1260 +1261 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] +1262 +1263 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) +1264 rmat = [] +1265 for t in range(basematrix.T): +1266 for i in range(Ntrunc): +1267 for j in range(Ntrunc): +1268 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] +1269 rmat.append(np.copy(tmpmat)) +1270 +1271 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] +1272 return Corr(newcontent) 1273 -1274def _sort_vectors(vec_set, ts): -1275 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" -1276 reference_sorting = np.array(vec_set[ts]) -1277 N = reference_sorting.shape[0] -1278 sorted_vec_set = [] -1279 for t in range(len(vec_set)): -1280 if vec_set[t] is None: -1281 sorted_vec_set.append(None) -1282 elif not t == ts: -1283 perms = [list(o) for o in permutations([i for i in range(N)], N)] -1284 best_score = 0 -1285 for perm in perms: -1286 current_score = 1 -1287 for k in range(N): -1288 new_sorting = reference_sorting.copy() -1289 new_sorting[perm[k], :] = vec_set[t][k] -1290 current_score *= abs(np.linalg.det(new_sorting)) -1291 if current_score > best_score: -1292 best_score = current_score -1293 best_perm = perm -1294 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) -1295 else: -1296 sorted_vec_set.append(vec_set[t]) -1297 -1298 return sorted_vec_set -1299 +1274 +1275def _sort_vectors(vec_set, ts): +1276 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" +1277 reference_sorting = np.array(vec_set[ts]) +1278 N = reference_sorting.shape[0] +1279 sorted_vec_set = [] +1280 for t in range(len(vec_set)): +1281 if vec_set[t] is None: +1282 sorted_vec_set.append(None) +1283 elif not t == ts: +1284 perms = [list(o) for o in permutations([i for i in range(N)], N)] +1285 best_score = 0 +1286 for perm in perms: +1287 current_score = 1 +1288 for k in range(N): +1289 new_sorting = reference_sorting.copy() +1290 new_sorting[perm[k], :] = vec_set[t][k] +1291 current_score *= abs(np.linalg.det(new_sorting)) +1292 if current_score > best_score: +1293 best_score = current_score +1294 best_perm = perm +1295 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) +1296 else: +1297 sorted_vec_set.append(vec_set[t]) +1298 +1299 return sorted_vec_set 1300 -1301def _check_for_none(corr, entry): -1302 """Checks if entry for correlator corr is None""" -1303 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2 -1304 +1301 +1302def _check_for_none(corr, entry): +1303 """Checks if entry for correlator corr is None""" +1304 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2 1305 -1306def _GEVP_solver(Gt, G0): -1307 """Helper function for solving the GEVP and sorting the eigenvectors. -1308 -1309 The helper function assumes that both provided matrices are symmetric and -1310 only processes the lower triangular part of both matrices. In case the matrices -1311 are not symmetric the upper triangular parts are effectively discarded.""" -1312 return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1] +1306 +1307def _GEVP_solver(Gt, G0): +1308 """Helper function for solving the GEVP and sorting the eigenvectors. +1309 +1310 The helper function assumes that both provided matrices are symmetric and +1311 only processes the lower triangular part of both matrices. In case the matrices +1312 are not symmetric the upper triangular parts are effectively discarded.""" +1313 return scipy.linalg.eigh(Gt, G0, lower=True)[1].T[::-1] @@ -1834,969 +1835,970 @@ 307 else: 308 symmetric_corr = self.matrix_symmetric() 309 - 310 if sort is None: - 311 if (ts is None): - 312 raise Exception("ts is required if sort=None.") - 313 if (self.content[t0] is None) or (self.content[ts] is None): - 314 raise Exception("Corr not defined at t0/ts.") - 315 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) - 316 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts]) - 317 reordered_vecs = _GEVP_solver(Gt, G0) - 318 - 319 elif sort in ["Eigenvalue", "Eigenvector"]: - 320 if sort == "Eigenvalue" and ts is not None: - 321 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) - 322 all_vecs = [None] * (t0 + 1) - 323 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) - 324 for t in range(t0 + 1, self.T): - 325 try: - 326 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t]) - 327 all_vecs.append(_GEVP_solver(Gt, G0)) - 328 except Exception: - 329 all_vecs.append(None) - 330 if sort == "Eigenvector": - 331 if (ts is None): - 332 raise Exception("ts is required for the Eigenvector sorting method.") - 333 all_vecs = _sort_vectors(all_vecs, ts) - 334 - 335 reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] - 336 else: - 337 raise Exception("Unkown value for 'sort'.") - 338 - 339 if "state" in kwargs: - 340 return reordered_vecs[kwargs.get("state")] - 341 else: - 342 return reordered_vecs - 343 - 344 def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"): - 345 """Determines the eigenvalue of the GEVP by solving and projecting the correlator - 346 - 347 Parameters - 348 ---------- - 349 state : int - 350 The state one is interested in ordered by energy. The lowest state is zero. - 351 - 352 All other parameters are identical to the ones of Corr.GEVP. - 353 """ - 354 vec = self.GEVP(t0, ts=ts, sort=sort)[state] - 355 return self.projected(vec) - 356 - 357 def Hankel(self, N, periodic=False): - 358 """Constructs an NxN Hankel matrix - 359 - 360 C(t) c(t+1) ... c(t+n-1) - 361 C(t+1) c(t+2) ... c(t+n) - 362 ................. - 363 C(t+(n-1)) c(t+n) ... c(t+2(n-1)) - 364 - 365 Parameters - 366 ---------- - 367 N : int - 368 Dimension of the Hankel matrix - 369 periodic : bool, optional - 370 determines whether the matrix is extended periodically - 371 """ - 372 - 373 if self.N != 1: - 374 raise Exception("Multi-operator Prony not implemented!") - 375 - 376 array = np.empty([N, N], dtype="object") - 377 new_content = [] - 378 for t in range(self.T): - 379 new_content.append(array.copy()) - 380 - 381 def wrap(i): - 382 while i >= self.T: - 383 i -= self.T - 384 return i - 385 - 386 for t in range(self.T): - 387 for i in range(N): - 388 for j in range(N): - 389 if periodic: - 390 new_content[t][i, j] = self.content[wrap(t + i + j)][0] - 391 elif (t + i + j) >= self.T: - 392 new_content[t] = None - 393 else: - 394 new_content[t][i, j] = self.content[t + i + j][0] - 395 - 396 return Corr(new_content) - 397 - 398 def roll(self, dt): - 399 """Periodically shift the correlator by dt timeslices - 400 - 401 Parameters - 402 ---------- - 403 dt : int - 404 number of timeslices - 405 """ - 406 return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) - 407 - 408 def reverse(self): - 409 """Reverse the time ordering of the Corr""" - 410 return Corr(self.content[:: -1]) - 411 - 412 def thin(self, spacing=2, offset=0): - 413 """Thin out a correlator to suppress correlations - 414 - 415 Parameters - 416 ---------- - 417 spacing : int - 418 Keep only every 'spacing'th entry of the correlator - 419 offset : int - 420 Offset the equal spacing - 421 """ - 422 new_content = [] - 423 for t in range(self.T): - 424 if (offset + t) % spacing != 0: - 425 new_content.append(None) - 426 else: - 427 new_content.append(self.content[t]) - 428 return Corr(new_content) - 429 - 430 def correlate(self, partner): - 431 """Correlate the correlator with another correlator or Obs - 432 - 433 Parameters - 434 ---------- - 435 partner : Obs or Corr - 436 partner to correlate the correlator with. - 437 Can either be an Obs which is correlated with all entries of the - 438 correlator or a Corr of same length. - 439 """ - 440 if self.N != 1: - 441 raise Exception("Only one-dimensional correlators can be safely correlated.") - 442 new_content = [] - 443 for x0, t_slice in enumerate(self.content): - 444 if _check_for_none(self, t_slice): - 445 new_content.append(None) - 446 else: - 447 if isinstance(partner, Corr): - 448 if _check_for_none(partner, partner.content[x0]): - 449 new_content.append(None) - 450 else: - 451 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) - 452 elif isinstance(partner, Obs): # Should this include CObs? - 453 new_content.append(np.array([correlate(o, partner) for o in t_slice])) - 454 else: - 455 raise Exception("Can only correlate with an Obs or a Corr.") - 456 - 457 return Corr(new_content) - 458 - 459 def reweight(self, weight, **kwargs): - 460 """Reweight the correlator. - 461 - 462 Parameters - 463 ---------- - 464 weight : Obs - 465 Reweighting factor. An Observable that has to be defined on a superset of the - 466 configurations in obs[i].idl for all i. - 467 all_configs : bool - 468 if True, the reweighted observables are normalized by the average of - 469 the reweighting factor on all configurations in weight.idl and not - 470 on the configurations in obs[i].idl. - 471 """ - 472 if self.N != 1: - 473 raise Exception("Reweighting only implemented for one-dimensional correlators.") - 474 new_content = [] - 475 for t_slice in self.content: - 476 if _check_for_none(self, t_slice): - 477 new_content.append(None) - 478 else: - 479 new_content.append(np.array(reweight(weight, t_slice, **kwargs))) - 480 return Corr(new_content) - 481 - 482 def T_symmetry(self, partner, parity=+1): - 483 """Return the time symmetry average of the correlator and its partner - 484 - 485 Parameters - 486 ---------- - 487 partner : Corr - 488 Time symmetry partner of the Corr - 489 partity : int - 490 Parity quantum number of the correlator, can be +1 or -1 - 491 """ - 492 if self.N != 1: - 493 raise Exception("T_symmetry only implemented for one-dimensional correlators.") - 494 if not isinstance(partner, Corr): - 495 raise Exception("T partner has to be a Corr object.") - 496 if parity not in [+1, -1]: - 497 raise Exception("Parity has to be +1 or -1.") - 498 T_partner = parity * partner.reverse() - 499 - 500 t_slices = [] - 501 test = (self - T_partner) - 502 test.gamma_method() - 503 for x0, t_slice in enumerate(test.content): - 504 if t_slice is not None: - 505 if not t_slice[0].is_zero_within_error(5): - 506 t_slices.append(x0) - 507 if t_slices: - 508 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning) - 509 - 510 return (self + T_partner) / 2 - 511 - 512 def deriv(self, variant="symmetric"): - 513 """Return the first derivative of the correlator with respect to x0. - 514 - 515 Parameters - 516 ---------- - 517 variant : str - 518 decides which definition of the finite differences derivative is used. - 519 Available choice: symmetric, forward, backward, improved, log, default: symmetric - 520 """ - 521 if self.N != 1: - 522 raise Exception("deriv only implemented for one-dimensional correlators.") - 523 if variant == "symmetric": - 524 newcontent = [] - 525 for t in range(1, self.T - 1): - 526 if (self.content[t - 1] is None) or (self.content[t + 1] is None): - 527 newcontent.append(None) - 528 else: - 529 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) - 530 if (all([x is None for x in newcontent])): - 531 raise Exception('Derivative is undefined at all timeslices') - 532 return Corr(newcontent, padding=[1, 1]) - 533 elif variant == "forward": - 534 newcontent = [] - 535 for t in range(self.T - 1): - 536 if (self.content[t] is None) or (self.content[t + 1] is None): - 537 newcontent.append(None) - 538 else: - 539 newcontent.append(self.content[t + 1] - self.content[t]) - 540 if (all([x is None for x in newcontent])): - 541 raise Exception("Derivative is undefined at all timeslices") - 542 return Corr(newcontent, padding=[0, 1]) - 543 elif variant == "backward": - 544 newcontent = [] - 545 for t in range(1, self.T): - 546 if (self.content[t - 1] is None) or (self.content[t] is None): - 547 newcontent.append(None) - 548 else: - 549 newcontent.append(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, 0]) - 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 + 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] - 8 * self.content[t - 1] + 8 * 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 elif variant == 'log': - 564 newcontent = [] - 565 for t in range(self.T): - 566 if (self.content[t] is None) or (self.content[t] <= 0): - 567 newcontent.append(None) - 568 else: - 569 newcontent.append(np.log(self.content[t])) - 570 if (all([x is None for x in newcontent])): - 571 raise Exception("Log is undefined at all timeslices") - 572 logcorr = Corr(newcontent) - 573 return self * logcorr.deriv('symmetric') - 574 else: - 575 raise Exception("Unknown variant.") - 576 - 577 def second_deriv(self, variant="symmetric"): - 578 """Return the second derivative of the correlator with respect to x0. - 579 - 580 Parameters - 581 ---------- - 582 variant : str - 583 decides which definition of the finite differences derivative is used. - 584 Available choice: symmetric, improved, log, default: symmetric - 585 """ - 586 if self.N != 1: - 587 raise Exception("second_deriv only implemented for one-dimensional correlators.") - 588 if variant == "symmetric": - 589 newcontent = [] - 590 for t in range(1, self.T - 1): - 591 if (self.content[t - 1] is None) or (self.content[t + 1] is None): - 592 newcontent.append(None) - 593 else: - 594 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) - 595 if (all([x is None for x in newcontent])): - 596 raise Exception("Derivative is undefined at all timeslices") - 597 return Corr(newcontent, padding=[1, 1]) - 598 elif variant == "improved": - 599 newcontent = [] - 600 for t in range(2, self.T - 2): - 601 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): - 602 newcontent.append(None) - 603 else: - 604 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])) - 605 if (all([x is None for x in newcontent])): - 606 raise Exception("Derivative is undefined at all timeslices") - 607 return Corr(newcontent, padding=[2, 2]) - 608 elif variant == 'log': - 609 newcontent = [] - 610 for t in range(self.T): - 611 if (self.content[t] is None) or (self.content[t] <= 0): - 612 newcontent.append(None) - 613 else: - 614 newcontent.append(np.log(self.content[t])) - 615 if (all([x is None for x in newcontent])): - 616 raise Exception("Log is undefined at all timeslices") - 617 logcorr = Corr(newcontent) - 618 return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2) - 619 else: - 620 raise Exception("Unknown variant.") - 621 - 622 def m_eff(self, variant='log', guess=1.0): - 623 """Returns the effective mass of the correlator as correlator object - 624 - 625 Parameters - 626 ---------- - 627 variant : str - 628 log : uses the standard effective mass log(C(t) / C(t+1)) - 629 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. - 630 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. - 631 See, e.g., arXiv:1205.5380 - 632 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) - 633 logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2 - 634 guess : float - 635 guess for the root finder, only relevant for the root variant - 636 """ - 637 if self.N != 1: - 638 raise Exception('Correlator must be projected before getting m_eff') - 639 if variant == 'log': - 640 newcontent = [] - 641 for t in range(self.T - 1): - 642 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): - 643 newcontent.append(None) - 644 elif self.content[t][0].value / self.content[t + 1][0].value < 0: - 645 newcontent.append(None) - 646 else: - 647 newcontent.append(self.content[t] / self.content[t + 1]) - 648 if (all([x is None for x in newcontent])): - 649 raise Exception('m_eff is undefined at all timeslices') - 650 - 651 return np.log(Corr(newcontent, padding=[0, 1])) - 652 - 653 elif variant == 'logsym': - 654 newcontent = [] - 655 for t in range(1, self.T - 1): - 656 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): - 657 newcontent.append(None) - 658 elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0: - 659 newcontent.append(None) - 660 else: - 661 newcontent.append(self.content[t - 1] / self.content[t + 1]) - 662 if (all([x is None for x in newcontent])): - 663 raise Exception('m_eff is undefined at all timeslices') - 664 - 665 return np.log(Corr(newcontent, padding=[1, 1])) / 2 - 666 - 667 elif variant in ['periodic', 'cosh', 'sinh']: - 668 if variant in ['periodic', 'cosh']: - 669 func = anp.cosh - 670 else: - 671 func = anp.sinh - 672 - 673 def root_function(x, d): - 674 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d - 675 - 676 newcontent = [] - 677 for t in range(self.T - 1): - 678 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): - 679 newcontent.append(None) - 680 # Fill the two timeslices in the middle of the lattice with their predecessors - 681 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: - 682 newcontent.append(newcontent[-1]) - 683 elif self.content[t][0].value / self.content[t + 1][0].value < 0: - 684 newcontent.append(None) - 685 else: - 686 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) - 687 if (all([x is None for x in newcontent])): - 688 raise Exception('m_eff is undefined at all timeslices') - 689 - 690 return Corr(newcontent, padding=[0, 1]) - 691 - 692 elif variant == 'arccosh': - 693 newcontent = [] - 694 for t in range(1, self.T - 1): - 695 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): - 696 newcontent.append(None) - 697 else: - 698 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) - 699 if (all([x is None for x in newcontent])): - 700 raise Exception("m_eff is undefined at all timeslices") - 701 return np.arccosh(Corr(newcontent, padding=[1, 1])) - 702 - 703 else: - 704 raise Exception('Unknown variant.') - 705 - 706 def fit(self, function, fitrange=None, silent=False, **kwargs): - 707 r'''Fits function to the data - 708 - 709 Parameters - 710 ---------- - 711 function : obj - 712 function to fit to the data. See fits.least_squares for details. - 713 fitrange : list - 714 Two element list containing the timeslices on which the fit is supposed to start and stop. - 715 Caution: This range is inclusive as opposed to standard python indexing. - 716 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. - 717 If not specified, self.prange or all timeslices are used. - 718 silent : bool - 719 Decides whether output is printed to the standard output. - 720 ''' - 721 if self.N != 1: - 722 raise Exception("Correlator must be projected before fitting") - 723 - 724 if fitrange is None: - 725 if self.prange: - 726 fitrange = self.prange - 727 else: - 728 fitrange = [0, self.T - 1] - 729 else: - 730 if not isinstance(fitrange, list): - 731 raise Exception("fitrange has to be a list with two elements") - 732 if len(fitrange) != 2: - 733 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") - 734 - 735 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] - 736 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] - 737 result = least_squares(xs, ys, function, silent=silent, **kwargs) - 738 return result - 739 - 740 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): - 741 """ Extract a plateau value from a Corr object - 742 - 743 Parameters - 744 ---------- - 745 plateau_range : list - 746 list with two entries, indicating the first and the last timeslice - 747 of the plateau region. - 748 method : str - 749 method to extract the plateau. - 750 'fit' fits a constant to the plateau region - 751 'avg', 'average' or 'mean' just average over the given timeslices. - 752 auto_gamma : bool - 753 apply gamma_method with default parameters to the Corr. Defaults to None - 754 """ - 755 if not plateau_range: - 756 if self.prange: - 757 plateau_range = self.prange - 758 else: - 759 raise Exception("no plateau range provided") - 760 if self.N != 1: - 761 raise Exception("Correlator must be projected before getting a plateau.") - 762 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): - 763 raise Exception("plateau is undefined at all timeslices in plateaurange.") - 764 if auto_gamma: - 765 self.gamma_method() - 766 if method == "fit": - 767 def const_func(a, t): - 768 return a[0] - 769 return self.fit(const_func, plateau_range)[0] - 770 elif method in ["avg", "average", "mean"]: - 771 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) - 772 return returnvalue - 773 - 774 else: - 775 raise Exception("Unsupported plateau method: " + method) - 776 - 777 def set_prange(self, prange): - 778 """Sets the attribute prange of the Corr object.""" - 779 if not len(prange) == 2: - 780 raise Exception("prange must be a list or array with two values") - 781 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): - 782 raise Exception("Start and end point must be integers") - 783 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): - 784 raise Exception("Start and end point must define a range in the interval 0,T") - 785 - 786 self.prange = prange - 787 return - 788 - 789 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): - 790 """Plots the correlator using the tag of the correlator as label if available. - 791 - 792 Parameters - 793 ---------- - 794 x_range : list - 795 list of two values, determining the range of the x-axis e.g. [4, 8]. - 796 comp : Corr or list of Corr - 797 Correlator or list of correlators which are plotted for comparison. - 798 The tags of these correlators are used as labels if available. - 799 logscale : bool - 800 Sets y-axis to logscale. - 801 plateau : Obs - 802 Plateau value to be visualized in the figure. - 803 fit_res : Fit_result - 804 Fit_result object to be visualized. - 805 ylabel : str - 806 Label for the y-axis. - 807 save : str - 808 path to file in which the figure should be saved. - 809 auto_gamma : bool - 810 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. - 811 hide_sigma : float - 812 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. - 813 references : list - 814 List of floating point values that are displayed as horizontal lines for reference. - 815 title : string - 816 Optional title of the figure. - 817 """ - 818 if self.N != 1: - 819 raise Exception("Correlator must be projected before plotting") - 820 - 821 if auto_gamma: - 822 self.gamma_method() - 823 - 824 if x_range is None: - 825 x_range = [0, self.T - 1] - 826 - 827 fig = plt.figure() - 828 ax1 = fig.add_subplot(111) - 829 - 830 x, y, y_err = self.plottable() - 831 if hide_sigma: - 832 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 - 833 else: - 834 hide_from = None - 835 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) - 836 if logscale: - 837 ax1.set_yscale('log') - 838 else: - 839 if y_range is None: - 840 try: - 841 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)]) - 842 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)]) - 843 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) - 844 except Exception: - 845 pass - 846 else: - 847 ax1.set_ylim(y_range) - 848 if comp: - 849 if isinstance(comp, (Corr, list)): - 850 for corr in comp if isinstance(comp, list) else [comp]: - 851 if auto_gamma: - 852 corr.gamma_method() - 853 x, y, y_err = corr.plottable() - 854 if hide_sigma: - 855 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 - 856 else: - 857 hide_from = None - 858 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) - 859 else: - 860 raise Exception("'comp' must be a correlator or a list of correlators.") - 861 - 862 if plateau: - 863 if isinstance(plateau, Obs): - 864 if auto_gamma: - 865 plateau.gamma_method() - 866 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) - 867 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') - 868 else: - 869 raise Exception("'plateau' must be an Obs") - 870 - 871 if references: - 872 if isinstance(references, list): - 873 for ref in references: - 874 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') - 875 else: - 876 raise Exception("'references' must be a list of floating pint values.") - 877 - 878 if self.prange: - 879 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') - 880 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') - 881 - 882 if fit_res: - 883 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) - 884 ax1.plot(x_samples, - 885 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), - 886 ls='-', marker=',', lw=2) - 887 - 888 ax1.set_xlabel(r'$x_0 / a$') - 889 if ylabel: - 890 ax1.set_ylabel(ylabel) - 891 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) - 892 - 893 handles, labels = ax1.get_legend_handles_labels() - 894 if labels: - 895 ax1.legend() - 896 - 897 if title: - 898 plt.title(title) - 899 - 900 plt.draw() - 901 - 902 if save: - 903 if isinstance(save, str): - 904 fig.savefig(save, bbox_inches='tight') - 905 else: - 906 raise Exception("'save' has to be a string.") - 907 - 908 def spaghetti_plot(self, logscale=True): - 909 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. - 910 - 911 Parameters - 912 ---------- - 913 logscale : bool - 914 Determines whether the scale of the y-axis is logarithmic or standard. - 915 """ - 916 if self.N != 1: - 917 raise Exception("Correlator needs to be projected first.") - 918 - 919 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])) - 920 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] - 921 - 922 for name in mc_names: - 923 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T - 924 - 925 fig = plt.figure() - 926 ax = fig.add_subplot(111) - 927 for dat in data: - 928 ax.plot(x0_vals, dat, ls='-', marker='') - 929 - 930 if logscale is True: - 931 ax.set_yscale('log') - 932 - 933 ax.set_xlabel(r'$x_0 / a$') - 934 plt.title(name) - 935 plt.draw() - 936 - 937 def dump(self, filename, datatype="json.gz", **kwargs): - 938 """Dumps the Corr into a file of chosen type - 939 Parameters - 940 ---------- - 941 filename : str - 942 Name of the file to be saved. - 943 datatype : str - 944 Format of the exported file. Supported formats include - 945 "json.gz" and "pickle" - 946 path : str - 947 specifies a custom path for the file (default '.') - 948 """ - 949 if datatype == "json.gz": - 950 from .input.json import dump_to_json - 951 if 'path' in kwargs: - 952 file_name = kwargs.get('path') + '/' + filename - 953 else: - 954 file_name = filename - 955 dump_to_json(self, file_name) - 956 elif datatype == "pickle": - 957 dump_object(self, filename, **kwargs) - 958 else: - 959 raise Exception("Unknown datatype " + str(datatype)) - 960 - 961 def print(self, print_range=None): - 962 print(self.__repr__(print_range)) - 963 - 964 def __repr__(self, print_range=None): - 965 if print_range is None: - 966 print_range = [0, None] - 967 - 968 content_string = "" - 969 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 - 970 - 971 if self.tag is not None: - 972 content_string += "Description: " + self.tag + "\n" - 973 if self.N != 1: - 974 return content_string - 975 - 976 if print_range[1]: - 977 print_range[1] += 1 - 978 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' - 979 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): - 980 if sub_corr is None: - 981 content_string += str(i + print_range[0]) + '\n' - 982 else: - 983 content_string += str(i + print_range[0]) - 984 for element in sub_corr: - 985 content_string += '\t' + ' ' * int(element >= 0) + str(element) - 986 content_string += '\n' - 987 return content_string - 988 - 989 def __str__(self): - 990 return self.__repr__() - 991 - 992 # We define the basic operations, that can be performed with correlators. - 993 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. - 994 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. - 995 # One could try and tell Obs to check if the y in __mul__ is a Corr and - 996 - 997 def __add__(self, y): - 998 if isinstance(y, Corr): - 999 if ((self.N != y.N) or (self.T != y.T)): -1000 raise Exception("Addition of Corrs with different shape") -1001 newcontent = [] -1002 for t in range(self.T): -1003 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1004 newcontent.append(None) -1005 else: -1006 newcontent.append(self.content[t] + y.content[t]) -1007 return Corr(newcontent) -1008 -1009 elif isinstance(y, (Obs, int, float, CObs)): -1010 newcontent = [] -1011 for t in range(self.T): -1012 if _check_for_none(self, self.content[t]): -1013 newcontent.append(None) -1014 else: -1015 newcontent.append(self.content[t] + y) -1016 return Corr(newcontent, prange=self.prange) -1017 elif isinstance(y, np.ndarray): -1018 if y.shape == (self.T,): -1019 return Corr(list((np.array(self.content).T + y).T)) -1020 else: -1021 raise ValueError("operands could not be broadcast together") -1022 else: -1023 raise TypeError("Corr + wrong type") -1024 -1025 def __mul__(self, y): -1026 if isinstance(y, Corr): -1027 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1028 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1029 newcontent = [] -1030 for t in range(self.T): -1031 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1032 newcontent.append(None) -1033 else: -1034 newcontent.append(self.content[t] * y.content[t]) -1035 return Corr(newcontent) -1036 -1037 elif isinstance(y, (Obs, int, float, CObs)): -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 __truediv__(self, y): -1054 if isinstance(y, Corr): -1055 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): -1056 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") -1057 newcontent = [] -1058 for t in range(self.T): -1059 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): -1060 newcontent.append(None) -1061 else: -1062 newcontent.append(self.content[t] / y.content[t]) -1063 for t in range(self.T): -1064 if _check_for_none(self, newcontent[t]): -1065 continue -1066 if np.isnan(np.sum(newcontent[t]).value): -1067 newcontent[t] = None -1068 -1069 if all([item is None for item in newcontent]): -1070 raise Exception("Division returns completely undefined correlator") -1071 return Corr(newcontent) -1072 -1073 elif isinstance(y, (Obs, CObs)): -1074 if isinstance(y, Obs): -1075 if y.value == 0: -1076 raise Exception('Division by zero will return undefined correlator') -1077 if isinstance(y, CObs): -1078 if y.is_zero(): -1079 raise Exception('Division by zero will return undefined correlator') -1080 -1081 newcontent = [] -1082 for t in range(self.T): -1083 if _check_for_none(self, self.content[t]): -1084 newcontent.append(None) -1085 else: -1086 newcontent.append(self.content[t] / y) -1087 return Corr(newcontent, prange=self.prange) -1088 -1089 elif isinstance(y, (int, float)): -1090 if y == 0: -1091 raise Exception('Division by zero will return undefined correlator') -1092 newcontent = [] -1093 for t in range(self.T): -1094 if _check_for_none(self, self.content[t]): -1095 newcontent.append(None) -1096 else: -1097 newcontent.append(self.content[t] / y) -1098 return Corr(newcontent, prange=self.prange) -1099 elif isinstance(y, np.ndarray): -1100 if y.shape == (self.T,): -1101 return Corr(list((np.array(self.content).T / y).T)) -1102 else: -1103 raise ValueError("operands could not be broadcast together") -1104 else: -1105 raise TypeError('Corr / wrong type') -1106 -1107 def __neg__(self): -1108 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] -1109 return Corr(newcontent, prange=self.prange) -1110 -1111 def __sub__(self, y): -1112 return self + (-y) -1113 -1114 def __pow__(self, y): -1115 if isinstance(y, (Obs, int, float, CObs)): -1116 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] -1117 return Corr(newcontent, prange=self.prange) -1118 else: -1119 raise TypeError('Type of exponent not supported') -1120 -1121 def __abs__(self): -1122 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] -1123 return Corr(newcontent, prange=self.prange) -1124 -1125 # The numpy functions: -1126 def sqrt(self): -1127 return self ** 0.5 -1128 -1129 def log(self): -1130 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] -1131 return Corr(newcontent, prange=self.prange) -1132 -1133 def exp(self): -1134 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] -1135 return Corr(newcontent, prange=self.prange) -1136 -1137 def _apply_func_to_corr(self, func): -1138 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] -1139 for t in range(self.T): -1140 if _check_for_none(self, newcontent[t]): -1141 continue -1142 if np.isnan(np.sum(newcontent[t]).value): -1143 newcontent[t] = None -1144 if all([item is None for item in newcontent]): -1145 raise Exception('Operation returns undefined correlator') -1146 return Corr(newcontent) -1147 -1148 def sin(self): -1149 return self._apply_func_to_corr(np.sin) -1150 -1151 def cos(self): -1152 return self._apply_func_to_corr(np.cos) -1153 -1154 def tan(self): -1155 return self._apply_func_to_corr(np.tan) -1156 -1157 def sinh(self): -1158 return self._apply_func_to_corr(np.sinh) -1159 -1160 def cosh(self): -1161 return self._apply_func_to_corr(np.cosh) -1162 -1163 def tanh(self): -1164 return self._apply_func_to_corr(np.tanh) -1165 -1166 def arcsin(self): -1167 return self._apply_func_to_corr(np.arcsin) -1168 -1169 def arccos(self): -1170 return self._apply_func_to_corr(np.arccos) -1171 -1172 def arctan(self): -1173 return self._apply_func_to_corr(np.arctan) -1174 -1175 def arcsinh(self): -1176 return self._apply_func_to_corr(np.arcsinh) -1177 -1178 def arccosh(self): -1179 return self._apply_func_to_corr(np.arccosh) -1180 -1181 def arctanh(self): -1182 return self._apply_func_to_corr(np.arctanh) -1183 -1184 # Right hand side operations (require tweak in main module to work) -1185 def __radd__(self, y): -1186 return self + y -1187 -1188 def __rsub__(self, y): -1189 return -self + y -1190 -1191 def __rmul__(self, y): -1192 return self * y -1193 -1194 def __rtruediv__(self, y): -1195 return (self / y) ** (-1) -1196 -1197 @property -1198 def real(self): -1199 def return_real(obs_OR_cobs): -1200 if isinstance(obs_OR_cobs, CObs): -1201 return obs_OR_cobs.real -1202 else: -1203 return obs_OR_cobs -1204 -1205 return self._apply_func_to_corr(return_real) -1206 -1207 @property -1208 def imag(self): -1209 def return_imag(obs_OR_cobs): -1210 if isinstance(obs_OR_cobs, CObs): -1211 return obs_OR_cobs.imag -1212 else: -1213 return obs_OR_cobs * 0 # So it stays the right type -1214 -1215 return self._apply_func_to_corr(return_imag) -1216 -1217 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): -1218 r''' Project large correlation matrix to lowest states -1219 -1220 This method can be used to reduce the size of an (N x N) correlation matrix -1221 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise -1222 is still small. -1223 -1224 Parameters -1225 ---------- -1226 Ntrunc: int -1227 Rank of the target matrix. -1228 tproj: int -1229 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. -1230 The default value is 3. -1231 t0proj: int -1232 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly -1233 discouraged for O(a) improved theories, since the correctness of the procedure -1234 cannot be granted in this case. The default value is 2. -1235 basematrix : Corr -1236 Correlation matrix that is used to determine the eigenvectors of the -1237 lowest states based on a GEVP. basematrix is taken to be the Corr itself if -1238 is is not specified. -1239 -1240 Notes -1241 ----- -1242 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving -1243 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}$ -1244 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the -1245 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via -1246 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large -1247 correlation matrix and to remove some noise that is added by irrelevant operators. -1248 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated -1249 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. -1250 ''' -1251 -1252 if self.N == 1: -1253 raise Exception('Method cannot be applied to one-dimensional correlators.') -1254 if basematrix is None: -1255 basematrix = self -1256 if Ntrunc >= basematrix.N: -1257 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) -1258 if basematrix.N != self.N: -1259 raise Exception('basematrix and targetmatrix have to be of the same size.') -1260 -1261 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] -1262 -1263 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) -1264 rmat = [] -1265 for t in range(basematrix.T): -1266 for i in range(Ntrunc): -1267 for j in range(Ntrunc): -1268 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] -1269 rmat.append(np.copy(tmpmat)) -1270 -1271 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] -1272 return Corr(newcontent) + 310 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) + 311 np.linalg.cholesky(G0) # Check if matrix G0 is positive-semidefinite. + 312 + 313 if sort is None: + 314 if (ts is None): + 315 raise Exception("ts is required if sort=None.") + 316 if (self.content[t0] is None) or (self.content[ts] is None): + 317 raise Exception("Corr not defined at t0/ts.") + 318 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts]) + 319 reordered_vecs = _GEVP_solver(Gt, G0) + 320 + 321 elif sort in ["Eigenvalue", "Eigenvector"]: + 322 if sort == "Eigenvalue" and ts is not None: + 323 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) + 324 all_vecs = [None] * (t0 + 1) + 325 for t in range(t0 + 1, self.T): + 326 try: + 327 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t]) + 328 all_vecs.append(_GEVP_solver(Gt, G0)) + 329 except Exception: + 330 all_vecs.append(None) + 331 if sort == "Eigenvector": + 332 if (ts is None): + 333 raise Exception("ts is required for the Eigenvector sorting method.") + 334 all_vecs = _sort_vectors(all_vecs, ts) + 335 + 336 reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] + 337 else: + 338 raise Exception("Unkown value for 'sort'.") + 339 + 340 if "state" in kwargs: + 341 return reordered_vecs[kwargs.get("state")] + 342 else: + 343 return reordered_vecs + 344 + 345 def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"): + 346 """Determines the eigenvalue of the GEVP by solving and projecting the correlator + 347 + 348 Parameters + 349 ---------- + 350 state : int + 351 The state one is interested in ordered by energy. The lowest state is zero. + 352 + 353 All other parameters are identical to the ones of Corr.GEVP. + 354 """ + 355 vec = self.GEVP(t0, ts=ts, sort=sort)[state] + 356 return self.projected(vec) + 357 + 358 def Hankel(self, N, periodic=False): + 359 """Constructs an NxN Hankel matrix + 360 + 361 C(t) c(t+1) ... c(t+n-1) + 362 C(t+1) c(t+2) ... c(t+n) + 363 ................. + 364 C(t+(n-1)) c(t+n) ... c(t+2(n-1)) + 365 + 366 Parameters + 367 ---------- + 368 N : int + 369 Dimension of the Hankel matrix + 370 periodic : bool, optional + 371 determines whether the matrix is extended periodically + 372 """ + 373 + 374 if self.N != 1: + 375 raise Exception("Multi-operator Prony not implemented!") + 376 + 377 array = np.empty([N, N], dtype="object") + 378 new_content = [] + 379 for t in range(self.T): + 380 new_content.append(array.copy()) + 381 + 382 def wrap(i): + 383 while i >= self.T: + 384 i -= self.T + 385 return i + 386 + 387 for t in range(self.T): + 388 for i in range(N): + 389 for j in range(N): + 390 if periodic: + 391 new_content[t][i, j] = self.content[wrap(t + i + j)][0] + 392 elif (t + i + j) >= self.T: + 393 new_content[t] = None + 394 else: + 395 new_content[t][i, j] = self.content[t + i + j][0] + 396 + 397 return Corr(new_content) + 398 + 399 def roll(self, dt): + 400 """Periodically shift the correlator by dt timeslices + 401 + 402 Parameters + 403 ---------- + 404 dt : int + 405 number of timeslices + 406 """ + 407 return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) + 408 + 409 def reverse(self): + 410 """Reverse the time ordering of the Corr""" + 411 return Corr(self.content[:: -1]) + 412 + 413 def thin(self, spacing=2, offset=0): + 414 """Thin out a correlator to suppress correlations + 415 + 416 Parameters + 417 ---------- + 418 spacing : int + 419 Keep only every 'spacing'th entry of the correlator + 420 offset : int + 421 Offset the equal spacing + 422 """ + 423 new_content = [] + 424 for t in range(self.T): + 425 if (offset + t) % spacing != 0: + 426 new_content.append(None) + 427 else: + 428 new_content.append(self.content[t]) + 429 return Corr(new_content) + 430 + 431 def correlate(self, partner): + 432 """Correlate the correlator with another correlator or Obs + 433 + 434 Parameters + 435 ---------- + 436 partner : Obs or Corr + 437 partner to correlate the correlator with. + 438 Can either be an Obs which is correlated with all entries of the + 439 correlator or a Corr of same length. + 440 """ + 441 if self.N != 1: + 442 raise Exception("Only one-dimensional correlators can be safely correlated.") + 443 new_content = [] + 444 for x0, t_slice in enumerate(self.content): + 445 if _check_for_none(self, t_slice): + 446 new_content.append(None) + 447 else: + 448 if isinstance(partner, Corr): + 449 if _check_for_none(partner, partner.content[x0]): + 450 new_content.append(None) + 451 else: + 452 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) + 453 elif isinstance(partner, Obs): # Should this include CObs? + 454 new_content.append(np.array([correlate(o, partner) for o in t_slice])) + 455 else: + 456 raise Exception("Can only correlate with an Obs or a Corr.") + 457 + 458 return Corr(new_content) + 459 + 460 def reweight(self, weight, **kwargs): + 461 """Reweight the correlator. + 462 + 463 Parameters + 464 ---------- + 465 weight : Obs + 466 Reweighting factor. An Observable that has to be defined on a superset of the + 467 configurations in obs[i].idl for all i. + 468 all_configs : bool + 469 if True, the reweighted observables are normalized by the average of + 470 the reweighting factor on all configurations in weight.idl and not + 471 on the configurations in obs[i].idl. + 472 """ + 473 if self.N != 1: + 474 raise Exception("Reweighting only implemented for one-dimensional correlators.") + 475 new_content = [] + 476 for t_slice in self.content: + 477 if _check_for_none(self, t_slice): + 478 new_content.append(None) + 479 else: + 480 new_content.append(np.array(reweight(weight, t_slice, **kwargs))) + 481 return Corr(new_content) + 482 + 483 def T_symmetry(self, partner, parity=+1): + 484 """Return the time symmetry average of the correlator and its partner + 485 + 486 Parameters + 487 ---------- + 488 partner : Corr + 489 Time symmetry partner of the Corr + 490 partity : int + 491 Parity quantum number of the correlator, can be +1 or -1 + 492 """ + 493 if self.N != 1: + 494 raise Exception("T_symmetry only implemented for one-dimensional correlators.") + 495 if not isinstance(partner, Corr): + 496 raise Exception("T partner has to be a Corr object.") + 497 if parity not in [+1, -1]: + 498 raise Exception("Parity has to be +1 or -1.") + 499 T_partner = parity * partner.reverse() + 500 + 501 t_slices = [] + 502 test = (self - T_partner) + 503 test.gamma_method() + 504 for x0, t_slice in enumerate(test.content): + 505 if t_slice is not None: + 506 if not t_slice[0].is_zero_within_error(5): + 507 t_slices.append(x0) + 508 if t_slices: + 509 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning) + 510 + 511 return (self + T_partner) / 2 + 512 + 513 def deriv(self, variant="symmetric"): + 514 """Return the first derivative of the correlator with respect to x0. + 515 + 516 Parameters + 517 ---------- + 518 variant : str + 519 decides which definition of the finite differences derivative is used. + 520 Available choice: symmetric, forward, backward, improved, log, default: symmetric + 521 """ + 522 if self.N != 1: + 523 raise Exception("deriv only implemented for one-dimensional correlators.") + 524 if variant == "symmetric": + 525 newcontent = [] + 526 for t in range(1, self.T - 1): + 527 if (self.content[t - 1] is None) or (self.content[t + 1] is None): + 528 newcontent.append(None) + 529 else: + 530 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) + 531 if (all([x is None for x in newcontent])): + 532 raise Exception('Derivative is undefined at all timeslices') + 533 return Corr(newcontent, padding=[1, 1]) + 534 elif variant == "forward": + 535 newcontent = [] + 536 for t in range(self.T - 1): + 537 if (self.content[t] is None) or (self.content[t + 1] is None): + 538 newcontent.append(None) + 539 else: + 540 newcontent.append(self.content[t + 1] - self.content[t]) + 541 if (all([x is None for x in newcontent])): + 542 raise Exception("Derivative is undefined at all timeslices") + 543 return Corr(newcontent, padding=[0, 1]) + 544 elif variant == "backward": + 545 newcontent = [] + 546 for t in range(1, self.T): + 547 if (self.content[t - 1] is None) or (self.content[t] is None): + 548 newcontent.append(None) + 549 else: + 550 newcontent.append(self.content[t] - self.content[t - 1]) + 551 if (all([x is None for x in newcontent])): + 552 raise Exception("Derivative is undefined at all timeslices") + 553 return Corr(newcontent, padding=[1, 0]) + 554 elif variant == "improved": + 555 newcontent = [] + 556 for t in range(2, self.T - 2): + 557 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): + 558 newcontent.append(None) + 559 else: + 560 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2])) + 561 if (all([x is None for x in newcontent])): + 562 raise Exception('Derivative is undefined at all timeslices') + 563 return Corr(newcontent, padding=[2, 2]) + 564 elif variant == 'log': + 565 newcontent = [] + 566 for t in range(self.T): + 567 if (self.content[t] is None) or (self.content[t] <= 0): + 568 newcontent.append(None) + 569 else: + 570 newcontent.append(np.log(self.content[t])) + 571 if (all([x is None for x in newcontent])): + 572 raise Exception("Log is undefined at all timeslices") + 573 logcorr = Corr(newcontent) + 574 return self * logcorr.deriv('symmetric') + 575 else: + 576 raise Exception("Unknown variant.") + 577 + 578 def second_deriv(self, variant="symmetric"): + 579 """Return the second derivative of the correlator with respect to x0. + 580 + 581 Parameters + 582 ---------- + 583 variant : str + 584 decides which definition of the finite differences derivative is used. + 585 Available choice: symmetric, improved, log, default: symmetric + 586 """ + 587 if self.N != 1: + 588 raise Exception("second_deriv only implemented for one-dimensional correlators.") + 589 if variant == "symmetric": + 590 newcontent = [] + 591 for t in range(1, self.T - 1): + 592 if (self.content[t - 1] is None) or (self.content[t + 1] is None): + 593 newcontent.append(None) + 594 else: + 595 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) + 596 if (all([x is None for x in newcontent])): + 597 raise Exception("Derivative is undefined at all timeslices") + 598 return Corr(newcontent, padding=[1, 1]) + 599 elif variant == "improved": + 600 newcontent = [] + 601 for t in range(2, self.T - 2): + 602 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): + 603 newcontent.append(None) + 604 else: + 605 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])) + 606 if (all([x is None for x in newcontent])): + 607 raise Exception("Derivative is undefined at all timeslices") + 608 return Corr(newcontent, padding=[2, 2]) + 609 elif variant == 'log': + 610 newcontent = [] + 611 for t in range(self.T): + 612 if (self.content[t] is None) or (self.content[t] <= 0): + 613 newcontent.append(None) + 614 else: + 615 newcontent.append(np.log(self.content[t])) + 616 if (all([x is None for x in newcontent])): + 617 raise Exception("Log is undefined at all timeslices") + 618 logcorr = Corr(newcontent) + 619 return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2) + 620 else: + 621 raise Exception("Unknown variant.") + 622 + 623 def m_eff(self, variant='log', guess=1.0): + 624 """Returns the effective mass of the correlator as correlator object + 625 + 626 Parameters + 627 ---------- + 628 variant : str + 629 log : uses the standard effective mass log(C(t) / C(t+1)) + 630 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. + 631 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. + 632 See, e.g., arXiv:1205.5380 + 633 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) + 634 logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2 + 635 guess : float + 636 guess for the root finder, only relevant for the root variant + 637 """ + 638 if self.N != 1: + 639 raise Exception('Correlator must be projected before getting m_eff') + 640 if variant == 'log': + 641 newcontent = [] + 642 for t in range(self.T - 1): + 643 if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): + 644 newcontent.append(None) + 645 elif self.content[t][0].value / self.content[t + 1][0].value < 0: + 646 newcontent.append(None) + 647 else: + 648 newcontent.append(self.content[t] / self.content[t + 1]) + 649 if (all([x is None for x in newcontent])): + 650 raise Exception('m_eff is undefined at all timeslices') + 651 + 652 return np.log(Corr(newcontent, padding=[0, 1])) + 653 + 654 elif variant == 'logsym': + 655 newcontent = [] + 656 for t in range(1, self.T - 1): + 657 if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0): + 658 newcontent.append(None) + 659 elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0: + 660 newcontent.append(None) + 661 else: + 662 newcontent.append(self.content[t - 1] / self.content[t + 1]) + 663 if (all([x is None for x in newcontent])): + 664 raise Exception('m_eff is undefined at all timeslices') + 665 + 666 return np.log(Corr(newcontent, padding=[1, 1])) / 2 + 667 + 668 elif variant in ['periodic', 'cosh', 'sinh']: + 669 if variant in ['periodic', 'cosh']: + 670 func = anp.cosh + 671 else: + 672 func = anp.sinh + 673 + 674 def root_function(x, d): + 675 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d + 676 + 677 newcontent = [] + 678 for t in range(self.T - 1): + 679 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0): + 680 newcontent.append(None) + 681 # Fill the two timeslices in the middle of the lattice with their predecessors + 682 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: + 683 newcontent.append(newcontent[-1]) + 684 elif self.content[t][0].value / self.content[t + 1][0].value < 0: + 685 newcontent.append(None) + 686 else: + 687 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) + 688 if (all([x is None for x in newcontent])): + 689 raise Exception('m_eff is undefined at all timeslices') + 690 + 691 return Corr(newcontent, padding=[0, 1]) + 692 + 693 elif variant == 'arccosh': + 694 newcontent = [] + 695 for t in range(1, self.T - 1): + 696 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): + 697 newcontent.append(None) + 698 else: + 699 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) + 700 if (all([x is None for x in newcontent])): + 701 raise Exception("m_eff is undefined at all timeslices") + 702 return np.arccosh(Corr(newcontent, padding=[1, 1])) + 703 + 704 else: + 705 raise Exception('Unknown variant.') + 706 + 707 def fit(self, function, fitrange=None, silent=False, **kwargs): + 708 r'''Fits function to the data + 709 + 710 Parameters + 711 ---------- + 712 function : obj + 713 function to fit to the data. See fits.least_squares for details. + 714 fitrange : list + 715 Two element list containing the timeslices on which the fit is supposed to start and stop. + 716 Caution: This range is inclusive as opposed to standard python indexing. + 717 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. + 718 If not specified, self.prange or all timeslices are used. + 719 silent : bool + 720 Decides whether output is printed to the standard output. + 721 ''' + 722 if self.N != 1: + 723 raise Exception("Correlator must be projected before fitting") + 724 + 725 if fitrange is None: + 726 if self.prange: + 727 fitrange = self.prange + 728 else: + 729 fitrange = [0, self.T - 1] + 730 else: + 731 if not isinstance(fitrange, list): + 732 raise Exception("fitrange has to be a list with two elements") + 733 if len(fitrange) != 2: + 734 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") + 735 + 736 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] + 737 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] + 738 result = least_squares(xs, ys, function, silent=silent, **kwargs) + 739 return result + 740 + 741 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): + 742 """ Extract a plateau value from a Corr object + 743 + 744 Parameters + 745 ---------- + 746 plateau_range : list + 747 list with two entries, indicating the first and the last timeslice + 748 of the plateau region. + 749 method : str + 750 method to extract the plateau. + 751 'fit' fits a constant to the plateau region + 752 'avg', 'average' or 'mean' just average over the given timeslices. + 753 auto_gamma : bool + 754 apply gamma_method with default parameters to the Corr. Defaults to None + 755 """ + 756 if not plateau_range: + 757 if self.prange: + 758 plateau_range = self.prange + 759 else: + 760 raise Exception("no plateau range provided") + 761 if self.N != 1: + 762 raise Exception("Correlator must be projected before getting a plateau.") + 763 if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): + 764 raise Exception("plateau is undefined at all timeslices in plateaurange.") + 765 if auto_gamma: + 766 self.gamma_method() + 767 if method == "fit": + 768 def const_func(a, t): + 769 return a[0] + 770 return self.fit(const_func, plateau_range)[0] + 771 elif method in ["avg", "average", "mean"]: + 772 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) + 773 return returnvalue + 774 + 775 else: + 776 raise Exception("Unsupported plateau method: " + method) + 777 + 778 def set_prange(self, prange): + 779 """Sets the attribute prange of the Corr object.""" + 780 if not len(prange) == 2: + 781 raise Exception("prange must be a list or array with two values") + 782 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): + 783 raise Exception("Start and end point must be integers") + 784 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): + 785 raise Exception("Start and end point must define a range in the interval 0,T") + 786 + 787 self.prange = prange + 788 return + 789 + 790 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): + 791 """Plots the correlator using the tag of the correlator as label if available. + 792 + 793 Parameters + 794 ---------- + 795 x_range : list + 796 list of two values, determining the range of the x-axis e.g. [4, 8]. + 797 comp : Corr or list of Corr + 798 Correlator or list of correlators which are plotted for comparison. + 799 The tags of these correlators are used as labels if available. + 800 logscale : bool + 801 Sets y-axis to logscale. + 802 plateau : Obs + 803 Plateau value to be visualized in the figure. + 804 fit_res : Fit_result + 805 Fit_result object to be visualized. + 806 ylabel : str + 807 Label for the y-axis. + 808 save : str + 809 path to file in which the figure should be saved. + 810 auto_gamma : bool + 811 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. + 812 hide_sigma : float + 813 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. + 814 references : list + 815 List of floating point values that are displayed as horizontal lines for reference. + 816 title : string + 817 Optional title of the figure. + 818 """ + 819 if self.N != 1: + 820 raise Exception("Correlator must be projected before plotting") + 821 + 822 if auto_gamma: + 823 self.gamma_method() + 824 + 825 if x_range is None: + 826 x_range = [0, self.T - 1] + 827 + 828 fig = plt.figure() + 829 ax1 = fig.add_subplot(111) + 830 + 831 x, y, y_err = self.plottable() + 832 if hide_sigma: + 833 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 + 834 else: + 835 hide_from = None + 836 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) + 837 if logscale: + 838 ax1.set_yscale('log') + 839 else: + 840 if y_range is None: + 841 try: + 842 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)]) + 843 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)]) + 844 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) + 845 except Exception: + 846 pass + 847 else: + 848 ax1.set_ylim(y_range) + 849 if comp: + 850 if isinstance(comp, (Corr, list)): + 851 for corr in comp if isinstance(comp, list) else [comp]: + 852 if auto_gamma: + 853 corr.gamma_method() + 854 x, y, y_err = corr.plottable() + 855 if hide_sigma: + 856 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 + 857 else: + 858 hide_from = None + 859 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) + 860 else: + 861 raise Exception("'comp' must be a correlator or a list of correlators.") + 862 + 863 if plateau: + 864 if isinstance(plateau, Obs): + 865 if auto_gamma: + 866 plateau.gamma_method() + 867 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) + 868 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') + 869 else: + 870 raise Exception("'plateau' must be an Obs") + 871 + 872 if references: + 873 if isinstance(references, list): + 874 for ref in references: + 875 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') + 876 else: + 877 raise Exception("'references' must be a list of floating pint values.") + 878 + 879 if self.prange: + 880 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') + 881 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') + 882 + 883 if fit_res: + 884 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) + 885 ax1.plot(x_samples, + 886 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), + 887 ls='-', marker=',', lw=2) + 888 + 889 ax1.set_xlabel(r'$x_0 / a$') + 890 if ylabel: + 891 ax1.set_ylabel(ylabel) + 892 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) + 893 + 894 handles, labels = ax1.get_legend_handles_labels() + 895 if labels: + 896 ax1.legend() + 897 + 898 if title: + 899 plt.title(title) + 900 + 901 plt.draw() + 902 + 903 if save: + 904 if isinstance(save, str): + 905 fig.savefig(save, bbox_inches='tight') + 906 else: + 907 raise Exception("'save' has to be a string.") + 908 + 909 def spaghetti_plot(self, logscale=True): + 910 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. + 911 + 912 Parameters + 913 ---------- + 914 logscale : bool + 915 Determines whether the scale of the y-axis is logarithmic or standard. + 916 """ + 917 if self.N != 1: + 918 raise Exception("Correlator needs to be projected first.") + 919 + 920 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])) + 921 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] + 922 + 923 for name in mc_names: + 924 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T + 925 + 926 fig = plt.figure() + 927 ax = fig.add_subplot(111) + 928 for dat in data: + 929 ax.plot(x0_vals, dat, ls='-', marker='') + 930 + 931 if logscale is True: + 932 ax.set_yscale('log') + 933 + 934 ax.set_xlabel(r'$x_0 / a$') + 935 plt.title(name) + 936 plt.draw() + 937 + 938 def dump(self, filename, datatype="json.gz", **kwargs): + 939 """Dumps the Corr into a file of chosen type + 940 Parameters + 941 ---------- + 942 filename : str + 943 Name of the file to be saved. + 944 datatype : str + 945 Format of the exported file. Supported formats include + 946 "json.gz" and "pickle" + 947 path : str + 948 specifies a custom path for the file (default '.') + 949 """ + 950 if datatype == "json.gz": + 951 from .input.json import dump_to_json + 952 if 'path' in kwargs: + 953 file_name = kwargs.get('path') + '/' + filename + 954 else: + 955 file_name = filename + 956 dump_to_json(self, file_name) + 957 elif datatype == "pickle": + 958 dump_object(self, filename, **kwargs) + 959 else: + 960 raise Exception("Unknown datatype " + str(datatype)) + 961 + 962 def print(self, print_range=None): + 963 print(self.__repr__(print_range)) + 964 + 965 def __repr__(self, print_range=None): + 966 if print_range is None: + 967 print_range = [0, None] + 968 + 969 content_string = "" + 970 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 + 971 + 972 if self.tag is not None: + 973 content_string += "Description: " + self.tag + "\n" + 974 if self.N != 1: + 975 return content_string + 976 + 977 if print_range[1]: + 978 print_range[1] += 1 + 979 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' + 980 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): + 981 if sub_corr is None: + 982 content_string += str(i + print_range[0]) + '\n' + 983 else: + 984 content_string += str(i + print_range[0]) + 985 for element in sub_corr: + 986 content_string += '\t' + ' ' * int(element >= 0) + str(element) + 987 content_string += '\n' + 988 return content_string + 989 + 990 def __str__(self): + 991 return self.__repr__() + 992 + 993 # We define the basic operations, that can be performed with correlators. + 994 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. + 995 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. + 996 # One could try and tell Obs to check if the y in __mul__ is a Corr and + 997 + 998 def __add__(self, y): + 999 if isinstance(y, Corr): +1000 if ((self.N != y.N) or (self.T != y.T)): +1001 raise Exception("Addition of Corrs with different shape") +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 return Corr(newcontent) +1009 +1010 elif isinstance(y, (Obs, int, float, CObs)): +1011 newcontent = [] +1012 for t in range(self.T): +1013 if _check_for_none(self, self.content[t]): +1014 newcontent.append(None) +1015 else: +1016 newcontent.append(self.content[t] + y) +1017 return Corr(newcontent, prange=self.prange) +1018 elif isinstance(y, np.ndarray): +1019 if y.shape == (self.T,): +1020 return Corr(list((np.array(self.content).T + y).T)) +1021 else: +1022 raise ValueError("operands could not be broadcast together") +1023 else: +1024 raise TypeError("Corr + wrong type") +1025 +1026 def __mul__(self, y): +1027 if isinstance(y, Corr): +1028 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1029 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1030 newcontent = [] +1031 for t in range(self.T): +1032 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1033 newcontent.append(None) +1034 else: +1035 newcontent.append(self.content[t] * y.content[t]) +1036 return Corr(newcontent) +1037 +1038 elif isinstance(y, (Obs, int, float, CObs)): +1039 newcontent = [] +1040 for t in range(self.T): +1041 if _check_for_none(self, self.content[t]): +1042 newcontent.append(None) +1043 else: +1044 newcontent.append(self.content[t] * y) +1045 return Corr(newcontent, prange=self.prange) +1046 elif isinstance(y, np.ndarray): +1047 if y.shape == (self.T,): +1048 return Corr(list((np.array(self.content).T * y).T)) +1049 else: +1050 raise ValueError("operands could not be broadcast together") +1051 else: +1052 raise TypeError("Corr * wrong type") +1053 +1054 def __truediv__(self, y): +1055 if isinstance(y, Corr): +1056 if not ((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): +1057 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") +1058 newcontent = [] +1059 for t in range(self.T): +1060 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): +1061 newcontent.append(None) +1062 else: +1063 newcontent.append(self.content[t] / y.content[t]) +1064 for t in range(self.T): +1065 if _check_for_none(self, newcontent[t]): +1066 continue +1067 if np.isnan(np.sum(newcontent[t]).value): +1068 newcontent[t] = None +1069 +1070 if all([item is None for item in newcontent]): +1071 raise Exception("Division returns completely undefined correlator") +1072 return Corr(newcontent) +1073 +1074 elif isinstance(y, (Obs, CObs)): +1075 if isinstance(y, Obs): +1076 if y.value == 0: +1077 raise Exception('Division by zero will return undefined correlator') +1078 if isinstance(y, CObs): +1079 if y.is_zero(): +1080 raise Exception('Division by zero will return undefined correlator') +1081 +1082 newcontent = [] +1083 for t in range(self.T): +1084 if _check_for_none(self, self.content[t]): +1085 newcontent.append(None) +1086 else: +1087 newcontent.append(self.content[t] / y) +1088 return Corr(newcontent, prange=self.prange) +1089 +1090 elif isinstance(y, (int, float)): +1091 if y == 0: +1092 raise Exception('Division by zero will return undefined correlator') +1093 newcontent = [] +1094 for t in range(self.T): +1095 if _check_for_none(self, self.content[t]): +1096 newcontent.append(None) +1097 else: +1098 newcontent.append(self.content[t] / y) +1099 return Corr(newcontent, prange=self.prange) +1100 elif isinstance(y, np.ndarray): +1101 if y.shape == (self.T,): +1102 return Corr(list((np.array(self.content).T / y).T)) +1103 else: +1104 raise ValueError("operands could not be broadcast together") +1105 else: +1106 raise TypeError('Corr / wrong type') +1107 +1108 def __neg__(self): +1109 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] +1110 return Corr(newcontent, prange=self.prange) +1111 +1112 def __sub__(self, y): +1113 return self + (-y) +1114 +1115 def __pow__(self, y): +1116 if isinstance(y, (Obs, int, float, CObs)): +1117 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] +1118 return Corr(newcontent, prange=self.prange) +1119 else: +1120 raise TypeError('Type of exponent not supported') +1121 +1122 def __abs__(self): +1123 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] +1124 return Corr(newcontent, prange=self.prange) +1125 +1126 # The numpy functions: +1127 def sqrt(self): +1128 return self ** 0.5 +1129 +1130 def log(self): +1131 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] +1132 return Corr(newcontent, prange=self.prange) +1133 +1134 def exp(self): +1135 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] +1136 return Corr(newcontent, prange=self.prange) +1137 +1138 def _apply_func_to_corr(self, func): +1139 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] +1140 for t in range(self.T): +1141 if _check_for_none(self, newcontent[t]): +1142 continue +1143 if np.isnan(np.sum(newcontent[t]).value): +1144 newcontent[t] = None +1145 if all([item is None for item in newcontent]): +1146 raise Exception('Operation returns undefined correlator') +1147 return Corr(newcontent) +1148 +1149 def sin(self): +1150 return self._apply_func_to_corr(np.sin) +1151 +1152 def cos(self): +1153 return self._apply_func_to_corr(np.cos) +1154 +1155 def tan(self): +1156 return self._apply_func_to_corr(np.tan) +1157 +1158 def sinh(self): +1159 return self._apply_func_to_corr(np.sinh) +1160 +1161 def cosh(self): +1162 return self._apply_func_to_corr(np.cosh) +1163 +1164 def tanh(self): +1165 return self._apply_func_to_corr(np.tanh) +1166 +1167 def arcsin(self): +1168 return self._apply_func_to_corr(np.arcsin) +1169 +1170 def arccos(self): +1171 return self._apply_func_to_corr(np.arccos) +1172 +1173 def arctan(self): +1174 return self._apply_func_to_corr(np.arctan) +1175 +1176 def arcsinh(self): +1177 return self._apply_func_to_corr(np.arcsinh) +1178 +1179 def arccosh(self): +1180 return self._apply_func_to_corr(np.arccosh) +1181 +1182 def arctanh(self): +1183 return self._apply_func_to_corr(np.arctanh) +1184 +1185 # Right hand side operations (require tweak in main module to work) +1186 def __radd__(self, y): +1187 return self + y +1188 +1189 def __rsub__(self, y): +1190 return -self + y +1191 +1192 def __rmul__(self, y): +1193 return self * y +1194 +1195 def __rtruediv__(self, y): +1196 return (self / y) ** (-1) +1197 +1198 @property +1199 def real(self): +1200 def return_real(obs_OR_cobs): +1201 if isinstance(obs_OR_cobs, CObs): +1202 return obs_OR_cobs.real +1203 else: +1204 return obs_OR_cobs +1205 +1206 return self._apply_func_to_corr(return_real) +1207 +1208 @property +1209 def imag(self): +1210 def return_imag(obs_OR_cobs): +1211 if isinstance(obs_OR_cobs, CObs): +1212 return obs_OR_cobs.imag +1213 else: +1214 return obs_OR_cobs * 0 # So it stays the right type +1215 +1216 return self._apply_func_to_corr(return_imag) +1217 +1218 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): +1219 r''' Project large correlation matrix to lowest states +1220 +1221 This method can be used to reduce the size of an (N x N) correlation matrix +1222 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise +1223 is still small. +1224 +1225 Parameters +1226 ---------- +1227 Ntrunc: int +1228 Rank of the target matrix. +1229 tproj: int +1230 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. +1231 The default value is 3. +1232 t0proj: int +1233 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly +1234 discouraged for O(a) improved theories, since the correctness of the procedure +1235 cannot be granted in this case. The default value is 2. +1236 basematrix : Corr +1237 Correlation matrix that is used to determine the eigenvectors of the +1238 lowest states based on a GEVP. basematrix is taken to be the Corr itself if +1239 is is not specified. +1240 +1241 Notes +1242 ----- +1243 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving +1244 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}$ +1245 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the +1246 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via +1247 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large +1248 correlation matrix and to remove some noise that is added by irrelevant operators. +1249 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated +1250 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. +1251 ''' +1252 +1253 if self.N == 1: +1254 raise Exception('Method cannot be applied to one-dimensional correlators.') +1255 if basematrix is None: +1256 basematrix = self +1257 if Ntrunc >= basematrix.N: +1258 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) +1259 if basematrix.N != self.N: +1260 raise Exception('basematrix and targetmatrix have to be of the same size.') +1261 +1262 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] +1263 +1264 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) +1265 rmat = [] +1266 for t in range(basematrix.T): +1267 for i in range(Ntrunc): +1268 for j in range(Ntrunc): +1269 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] +1270 rmat.append(np.copy(tmpmat)) +1271 +1272 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] +1273 return Corr(newcontent) @@ -3278,39 +3280,40 @@ timeslice and the error on each timeslice.

307 else: 308 symmetric_corr = self.matrix_symmetric() 309 -310 if sort is None: -311 if (ts is None): -312 raise Exception("ts is required if sort=None.") -313 if (self.content[t0] is None) or (self.content[ts] is None): -314 raise Exception("Corr not defined at t0/ts.") -315 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) -316 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts]) -317 reordered_vecs = _GEVP_solver(Gt, G0) -318 -319 elif sort in ["Eigenvalue", "Eigenvector"]: -320 if sort == "Eigenvalue" and ts is not None: -321 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) -322 all_vecs = [None] * (t0 + 1) -323 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) -324 for t in range(t0 + 1, self.T): -325 try: -326 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t]) -327 all_vecs.append(_GEVP_solver(Gt, G0)) -328 except Exception: -329 all_vecs.append(None) -330 if sort == "Eigenvector": -331 if (ts is None): -332 raise Exception("ts is required for the Eigenvector sorting method.") -333 all_vecs = _sort_vectors(all_vecs, ts) -334 -335 reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] -336 else: -337 raise Exception("Unkown value for 'sort'.") -338 -339 if "state" in kwargs: -340 return reordered_vecs[kwargs.get("state")] -341 else: -342 return reordered_vecs +310 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) +311 np.linalg.cholesky(G0) # Check if matrix G0 is positive-semidefinite. +312 +313 if sort is None: +314 if (ts is None): +315 raise Exception("ts is required if sort=None.") +316 if (self.content[t0] is None) or (self.content[ts] is None): +317 raise Exception("Corr not defined at t0/ts.") +318 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts]) +319 reordered_vecs = _GEVP_solver(Gt, G0) +320 +321 elif sort in ["Eigenvalue", "Eigenvector"]: +322 if sort == "Eigenvalue" and ts is not None: +323 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) +324 all_vecs = [None] * (t0 + 1) +325 for t in range(t0 + 1, self.T): +326 try: +327 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t]) +328 all_vecs.append(_GEVP_solver(Gt, G0)) +329 except Exception: +330 all_vecs.append(None) +331 if sort == "Eigenvector": +332 if (ts is None): +333 raise Exception("ts is required for the Eigenvector sorting method.") +334 all_vecs = _sort_vectors(all_vecs, ts) +335 +336 reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] +337 else: +338 raise Exception("Unkown value for 'sort'.") +339 +340 if "state" in kwargs: +341 return reordered_vecs[kwargs.get("state")] +342 else: +343 return reordered_vecs @@ -3363,18 +3366,18 @@ Returns only the vector(s) for a specified state. The lowest state is zero. - 
344    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
-345        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
-346
-347        Parameters
-348        ----------
-349        state : int
-350            The state one is interested in ordered by energy. The lowest state is zero.
-351
-352        All other parameters are identical to the ones of Corr.GEVP.
-353        """
-354        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
-355        return self.projected(vec)
+            
345    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
+346        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
+347
+348        Parameters
+349        ----------
+350        state : int
+351            The state one is interested in ordered by energy. The lowest state is zero.
+352
+353        All other parameters are identical to the ones of Corr.GEVP.
+354        """
+355        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
+356        return self.projected(vec)
 

 
 
@@ -3402,46 +3405,46 @@ The state one is interested in ordered by energy. The lowest state is zero.
- 
357    def Hankel(self, N, periodic=False):
-358        """Constructs an NxN Hankel matrix
-359
-360        C(t) c(t+1) ... c(t+n-1)
-361        C(t+1) c(t+2) ... c(t+n)
-362        .................
-363        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
-364
-365        Parameters
-366        ----------
-367        N : int
-368            Dimension of the Hankel matrix
-369        periodic : bool, optional
-370            determines whether the matrix is extended periodically
-371        """
-372
-373        if self.N != 1:
-374            raise Exception("Multi-operator Prony not implemented!")
-375
-376        array = np.empty([N, N], dtype="object")
-377        new_content = []
-378        for t in range(self.T):
-379            new_content.append(array.copy())
-380
-381        def wrap(i):
-382            while i >= self.T:
-383                i -= self.T
-384            return i
-385
-386        for t in range(self.T):
-387            for i in range(N):
-388                for j in range(N):
-389                    if periodic:
-390                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
-391                    elif (t + i + j) >= self.T:
-392                        new_content[t] = None
-393                    else:
-394                        new_content[t][i, j] = self.content[t + i + j][0]
-395
-396        return Corr(new_content)
+            
358    def Hankel(self, N, periodic=False):
+359        """Constructs an NxN Hankel matrix
+360
+361        C(t) c(t+1) ... c(t+n-1)
+362        C(t+1) c(t+2) ... c(t+n)
+363        .................
+364        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
+365
+366        Parameters
+367        ----------
+368        N : int
+369            Dimension of the Hankel matrix
+370        periodic : bool, optional
+371            determines whether the matrix is extended periodically
+372        """
+373
+374        if self.N != 1:
+375            raise Exception("Multi-operator Prony not implemented!")
+376
+377        array = np.empty([N, N], dtype="object")
+378        new_content = []
+379        for t in range(self.T):
+380            new_content.append(array.copy())
+381
+382        def wrap(i):
+383            while i >= self.T:
+384                i -= self.T
+385            return i
+386
+387        for t in range(self.T):
+388            for i in range(N):
+389                for j in range(N):
+390                    if periodic:
+391                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
+392                    elif (t + i + j) >= self.T:
+393                        new_content[t] = None
+394                    else:
+395                        new_content[t][i, j] = self.content[t + i + j][0]
+396
+397        return Corr(new_content)
 

 
 
@@ -3475,15 +3478,15 @@ determines whether the matrix is extended periodically
- 
398    def roll(self, dt):
-399        """Periodically shift the correlator by dt timeslices
-400
-401        Parameters
-402        ----------
-403        dt : int
-404            number of timeslices
-405        """
-406        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
+            
399    def roll(self, dt):
+400        """Periodically shift the correlator by dt timeslices
+401
+402        Parameters
+403        ----------
+404        dt : int
+405            number of timeslices
+406        """
+407        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
 

 
 
@@ -3510,9 +3513,9 @@ number of timeslices
- 
408    def reverse(self):
-409        """Reverse the time ordering of the Corr"""
-410        return Corr(self.content[:: -1])
+            
409    def reverse(self):
+410        """Reverse the time ordering of the Corr"""
+411        return Corr(self.content[:: -1])
 

 
 
@@ -3532,23 +3535,23 @@ number of timeslices
- 
412    def thin(self, spacing=2, offset=0):
-413        """Thin out a correlator to suppress correlations
-414
-415        Parameters
-416        ----------
-417        spacing : int
-418            Keep only every 'spacing'th entry of the correlator
-419        offset : int
-420            Offset the equal spacing
-421        """
-422        new_content = []
-423        for t in range(self.T):
-424            if (offset + t) % spacing != 0:
-425                new_content.append(None)
-426            else:
-427                new_content.append(self.content[t])
-428        return Corr(new_content)
+            
413    def thin(self, spacing=2, offset=0):
+414        """Thin out a correlator to suppress correlations
+415
+416        Parameters
+417        ----------
+418        spacing : int
+419            Keep only every 'spacing'th entry of the correlator
+420        offset : int
+421            Offset the equal spacing
+422        """
+423        new_content = []
+424        for t in range(self.T):
+425            if (offset + t) % spacing != 0:
+426                new_content.append(None)
+427            else:
+428                new_content.append(self.content[t])
+429        return Corr(new_content)
 

 
 
@@ -3577,34 +3580,34 @@ Offset the equal spacing
- 
430    def correlate(self, partner):
-431        """Correlate the correlator with another correlator or Obs
-432
-433        Parameters
-434        ----------
-435        partner : Obs or Corr
-436            partner to correlate the correlator with.
-437            Can either be an Obs which is correlated with all entries of the
-438            correlator or a Corr of same length.
-439        """
-440        if self.N != 1:
-441            raise Exception("Only one-dimensional correlators can be safely correlated.")
-442        new_content = []
-443        for x0, t_slice in enumerate(self.content):
-444            if _check_for_none(self, t_slice):
-445                new_content.append(None)
-446            else:
-447                if isinstance(partner, Corr):
-448                    if _check_for_none(partner, partner.content[x0]):
-449                        new_content.append(None)
-450                    else:
-451                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
-452                elif isinstance(partner, Obs):  # Should this include CObs?
-453                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
-454                else:
-455                    raise Exception("Can only correlate with an Obs or a Corr.")
-456
-457        return Corr(new_content)
+            
431    def correlate(self, partner):
+432        """Correlate the correlator with another correlator or Obs
+433
+434        Parameters
+435        ----------
+436        partner : Obs or Corr
+437            partner to correlate the correlator with.
+438            Can either be an Obs which is correlated with all entries of the
+439            correlator or a Corr of same length.
+440        """
+441        if self.N != 1:
+442            raise Exception("Only one-dimensional correlators can be safely correlated.")
+443        new_content = []
+444        for x0, t_slice in enumerate(self.content):
+445            if _check_for_none(self, t_slice):
+446                new_content.append(None)
+447            else:
+448                if isinstance(partner, Corr):
+449                    if _check_for_none(partner, partner.content[x0]):
+450                        new_content.append(None)
+451                    else:
+452                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
+453                elif isinstance(partner, Obs):  # Should this include CObs?
+454                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
+455                else:
+456                    raise Exception("Can only correlate with an Obs or a Corr.")
+457
+458        return Corr(new_content)
 

 
 
@@ -3633,28 +3636,28 @@ correlator or a Corr of same length.
- 
459    def reweight(self, weight, **kwargs):
-460        """Reweight the correlator.
-461
-462        Parameters
-463        ----------
-464        weight : Obs
-465            Reweighting factor. An Observable that has to be defined on a superset of the
-466            configurations in obs[i].idl for all i.
-467        all_configs : bool
-468            if True, the reweighted observables are normalized by the average of
-469            the reweighting factor on all configurations in weight.idl and not
-470            on the configurations in obs[i].idl.
-471        """
-472        if self.N != 1:
-473            raise Exception("Reweighting only implemented for one-dimensional correlators.")
-474        new_content = []
-475        for t_slice in self.content:
-476            if _check_for_none(self, t_slice):
-477                new_content.append(None)
-478            else:
-479                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
-480        return Corr(new_content)
+            
460    def reweight(self, weight, **kwargs):
+461        """Reweight the correlator.
+462
+463        Parameters
+464        ----------
+465        weight : Obs
+466            Reweighting factor. An Observable that has to be defined on a superset of the
+467            configurations in obs[i].idl for all i.
+468        all_configs : bool
+469            if True, the reweighted observables are normalized by the average of
+470            the reweighting factor on all configurations in weight.idl and not
+471            on the configurations in obs[i].idl.
+472        """
+473        if self.N != 1:
+474            raise Exception("Reweighting only implemented for one-dimensional correlators.")
+475        new_content = []
+476        for t_slice in self.content:
+477            if _check_for_none(self, t_slice):
+478                new_content.append(None)
+479            else:
+480                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
+481        return Corr(new_content)
 

 
 
@@ -3686,35 +3689,35 @@ on the configurations in obs[i].idl.
- 
482    def T_symmetry(self, partner, parity=+1):
-483        """Return the time symmetry average of the correlator and its partner
-484
-485        Parameters
-486        ----------
-487        partner : Corr
-488            Time symmetry partner of the Corr
-489        partity : int
-490            Parity quantum number of the correlator, can be +1 or -1
-491        """
-492        if self.N != 1:
-493            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
-494        if not isinstance(partner, Corr):
-495            raise Exception("T partner has to be a Corr object.")
-496        if parity not in [+1, -1]:
-497            raise Exception("Parity has to be +1 or -1.")
-498        T_partner = parity * partner.reverse()
-499
-500        t_slices = []
-501        test = (self - T_partner)
-502        test.gamma_method()
-503        for x0, t_slice in enumerate(test.content):
-504            if t_slice is not None:
-505                if not t_slice[0].is_zero_within_error(5):
-506                    t_slices.append(x0)
-507        if t_slices:
-508            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
-509
-510        return (self + T_partner) / 2
+            
483    def T_symmetry(self, partner, parity=+1):
+484        """Return the time symmetry average of the correlator and its partner
+485
+486        Parameters
+487        ----------
+488        partner : Corr
+489            Time symmetry partner of the Corr
+490        partity : int
+491            Parity quantum number of the correlator, can be +1 or -1
+492        """
+493        if self.N != 1:
+494            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
+495        if not isinstance(partner, Corr):
+496            raise Exception("T partner has to be a Corr object.")
+497        if parity not in [+1, -1]:
+498            raise Exception("Parity has to be +1 or -1.")
+499        T_partner = parity * partner.reverse()
+500
+501        t_slices = []
+502        test = (self - T_partner)
+503        test.gamma_method()
+504        for x0, t_slice in enumerate(test.content):
+505            if t_slice is not None:
+506                if not t_slice[0].is_zero_within_error(5):
+507                    t_slices.append(x0)
+508        if t_slices:
+509            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
+510
+511        return (self + T_partner) / 2
 

 
 
@@ -3743,70 +3746,70 @@ Parity quantum number of the correlator, can be +1 or -1
- 
512    def deriv(self, variant="symmetric"):
-513        """Return the first derivative of the correlator with respect to x0.
-514
-515        Parameters
-516        ----------
-517        variant : str
-518            decides which definition of the finite differences derivative is used.
-519            Available choice: symmetric, forward, backward, improved, log, default: symmetric
-520        """
-521        if self.N != 1:
-522            raise Exception("deriv only implemented for one-dimensional correlators.")
-523        if variant == "symmetric":
-524            newcontent = []
-525            for t in range(1, self.T - 1):
-526                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-527                    newcontent.append(None)
-528                else:
-529                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
-530            if (all([x is None for x in newcontent])):
-531                raise Exception('Derivative is undefined at all timeslices')
-532            return Corr(newcontent, padding=[1, 1])
-533        elif variant == "forward":
-534            newcontent = []
-535            for t in range(self.T - 1):
-536                if (self.content[t] is None) or (self.content[t + 1] is None):
-537                    newcontent.append(None)
-538                else:
-539                    newcontent.append(self.content[t + 1] - self.content[t])
-540            if (all([x is None for x in newcontent])):
-541                raise Exception("Derivative is undefined at all timeslices")
-542            return Corr(newcontent, padding=[0, 1])
-543        elif variant == "backward":
-544            newcontent = []
-545            for t in range(1, self.T):
-546                if (self.content[t - 1] is None) or (self.content[t] is None):
-547                    newcontent.append(None)
-548                else:
-549                    newcontent.append(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, 0])
-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 + 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] - 8 * self.content[t - 1] + 8 * 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        elif variant == 'log':
-564            newcontent = []
-565            for t in range(self.T):
-566                if (self.content[t] is None) or (self.content[t] <= 0):
-567                    newcontent.append(None)
-568                else:
-569                    newcontent.append(np.log(self.content[t]))
-570            if (all([x is None for x in newcontent])):
-571                raise Exception("Log is undefined at all timeslices")
-572            logcorr = Corr(newcontent)
-573            return self * logcorr.deriv('symmetric')
-574        else:
-575            raise Exception("Unknown variant.")
+            
513    def deriv(self, variant="symmetric"):
+514        """Return the first derivative of the correlator with respect to x0.
+515
+516        Parameters
+517        ----------
+518        variant : str
+519            decides which definition of the finite differences derivative is used.
+520            Available choice: symmetric, forward, backward, improved, log, default: symmetric
+521        """
+522        if self.N != 1:
+523            raise Exception("deriv only implemented for one-dimensional correlators.")
+524        if variant == "symmetric":
+525            newcontent = []
+526            for t in range(1, self.T - 1):
+527                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+528                    newcontent.append(None)
+529                else:
+530                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
+531            if (all([x is None for x in newcontent])):
+532                raise Exception('Derivative is undefined at all timeslices')
+533            return Corr(newcontent, padding=[1, 1])
+534        elif variant == "forward":
+535            newcontent = []
+536            for t in range(self.T - 1):
+537                if (self.content[t] is None) or (self.content[t + 1] is None):
+538                    newcontent.append(None)
+539                else:
+540                    newcontent.append(self.content[t + 1] - self.content[t])
+541            if (all([x is None for x in newcontent])):
+542                raise Exception("Derivative is undefined at all timeslices")
+543            return Corr(newcontent, padding=[0, 1])
+544        elif variant == "backward":
+545            newcontent = []
+546            for t in range(1, self.T):
+547                if (self.content[t - 1] is None) or (self.content[t] is None):
+548                    newcontent.append(None)
+549                else:
+550                    newcontent.append(self.content[t] - self.content[t - 1])
+551            if (all([x is None for x in newcontent])):
+552                raise Exception("Derivative is undefined at all timeslices")
+553            return Corr(newcontent, padding=[1, 0])
+554        elif variant == "improved":
+555            newcontent = []
+556            for t in range(2, self.T - 2):
+557                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):
+558                    newcontent.append(None)
+559                else:
+560                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
+561            if (all([x is None for x in newcontent])):
+562                raise Exception('Derivative is undefined at all timeslices')
+563            return Corr(newcontent, padding=[2, 2])
+564        elif variant == 'log':
+565            newcontent = []
+566            for t in range(self.T):
+567                if (self.content[t] is None) or (self.content[t] <= 0):
+568                    newcontent.append(None)
+569                else:
+570                    newcontent.append(np.log(self.content[t]))
+571            if (all([x is None for x in newcontent])):
+572                raise Exception("Log is undefined at all timeslices")
+573            logcorr = Corr(newcontent)
+574            return self * logcorr.deriv('symmetric')
+575        else:
+576            raise Exception("Unknown variant.")
 

 
 
@@ -3834,50 +3837,50 @@ Available choice: symmetric, forward, backward, improved, log, default: symmetri
- 
577    def second_deriv(self, variant="symmetric"):
-578        """Return the second derivative of the correlator with respect to x0.
-579
-580        Parameters
-581        ----------
-582        variant : str
-583            decides which definition of the finite differences derivative is used.
-584            Available choice: symmetric, improved, log, default: symmetric
-585        """
-586        if self.N != 1:
-587            raise Exception("second_deriv only implemented for one-dimensional correlators.")
-588        if variant == "symmetric":
-589            newcontent = []
-590            for t in range(1, self.T - 1):
-591                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-592                    newcontent.append(None)
-593                else:
-594                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
-595            if (all([x is None for x in newcontent])):
-596                raise Exception("Derivative is undefined at all timeslices")
-597            return Corr(newcontent, padding=[1, 1])
-598        elif variant == "improved":
-599            newcontent = []
-600            for t in range(2, self.T - 2):
-601                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):
-602                    newcontent.append(None)
-603                else:
-604                    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]))
-605            if (all([x is None for x in newcontent])):
-606                raise Exception("Derivative is undefined at all timeslices")
-607            return Corr(newcontent, padding=[2, 2])
-608        elif variant == 'log':
-609            newcontent = []
-610            for t in range(self.T):
-611                if (self.content[t] is None) or (self.content[t] <= 0):
-612                    newcontent.append(None)
-613                else:
-614                    newcontent.append(np.log(self.content[t]))
-615            if (all([x is None for x in newcontent])):
-616                raise Exception("Log is undefined at all timeslices")
-617            logcorr = Corr(newcontent)
-618            return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2)
-619        else:
-620            raise Exception("Unknown variant.")
+            
578    def second_deriv(self, variant="symmetric"):
+579        """Return the second derivative of the correlator with respect to x0.
+580
+581        Parameters
+582        ----------
+583        variant : str
+584            decides which definition of the finite differences derivative is used.
+585            Available choice: symmetric, improved, log, default: symmetric
+586        """
+587        if self.N != 1:
+588            raise Exception("second_deriv only implemented for one-dimensional correlators.")
+589        if variant == "symmetric":
+590            newcontent = []
+591            for t in range(1, self.T - 1):
+592                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+593                    newcontent.append(None)
+594                else:
+595                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
+596            if (all([x is None for x in newcontent])):
+597                raise Exception("Derivative is undefined at all timeslices")
+598            return Corr(newcontent, padding=[1, 1])
+599        elif variant == "improved":
+600            newcontent = []
+601            for t in range(2, self.T - 2):
+602                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):
+603                    newcontent.append(None)
+604                else:
+605                    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]))
+606            if (all([x is None for x in newcontent])):
+607                raise Exception("Derivative is undefined at all timeslices")
+608            return Corr(newcontent, padding=[2, 2])
+609        elif variant == 'log':
+610            newcontent = []
+611            for t in range(self.T):
+612                if (self.content[t] is None) or (self.content[t] <= 0):
+613                    newcontent.append(None)
+614                else:
+615                    newcontent.append(np.log(self.content[t]))
+616            if (all([x is None for x in newcontent])):
+617                raise Exception("Log is undefined at all timeslices")
+618            logcorr = Corr(newcontent)
+619            return self * (logcorr.second_deriv('symmetric') + (logcorr.deriv('symmetric'))**2)
+620        else:
+621            raise Exception("Unknown variant.")
 

 
 
@@ -3905,89 +3908,89 @@ Available choice: symmetric, improved, log, default: symmetric
- 
622    def m_eff(self, variant='log', guess=1.0):
-623        """Returns the effective mass of the correlator as correlator object
-624
-625        Parameters
-626        ----------
-627        variant : str
-628            log : uses the standard effective mass log(C(t) / C(t+1))
-629            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.
-630            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.
-631            See, e.g., arXiv:1205.5380
-632            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
-633            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
-634        guess : float
-635            guess for the root finder, only relevant for the root variant
-636        """
-637        if self.N != 1:
-638            raise Exception('Correlator must be projected before getting m_eff')
-639        if variant == 'log':
-640            newcontent = []
-641            for t in range(self.T - 1):
-642                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
-643                    newcontent.append(None)
-644                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
-645                    newcontent.append(None)
-646                else:
-647                    newcontent.append(self.content[t] / self.content[t + 1])
-648            if (all([x is None for x in newcontent])):
-649                raise Exception('m_eff is undefined at all timeslices')
-650
-651            return np.log(Corr(newcontent, padding=[0, 1]))
-652
-653        elif variant == 'logsym':
-654            newcontent = []
-655            for t in range(1, self.T - 1):
-656                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
-657                    newcontent.append(None)
-658                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
-659                    newcontent.append(None)
-660                else:
-661                    newcontent.append(self.content[t - 1] / self.content[t + 1])
-662            if (all([x is None for x in newcontent])):
-663                raise Exception('m_eff is undefined at all timeslices')
-664
-665            return np.log(Corr(newcontent, padding=[1, 1])) / 2
-666
-667        elif variant in ['periodic', 'cosh', 'sinh']:
-668            if variant in ['periodic', 'cosh']:
-669                func = anp.cosh
-670            else:
-671                func = anp.sinh
-672
-673            def root_function(x, d):
-674                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
-675
-676            newcontent = []
-677            for t in range(self.T - 1):
-678                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
-679                    newcontent.append(None)
-680                # Fill the two timeslices in the middle of the lattice with their predecessors
-681                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
-682                    newcontent.append(newcontent[-1])
-683                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
-684                    newcontent.append(None)
-685                else:
-686                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
-687            if (all([x is None for x in newcontent])):
-688                raise Exception('m_eff is undefined at all timeslices')
-689
-690            return Corr(newcontent, padding=[0, 1])
-691
-692        elif variant == 'arccosh':
-693            newcontent = []
-694            for t in range(1, self.T - 1):
-695                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):
-696                    newcontent.append(None)
-697                else:
-698                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
-699            if (all([x is None for x in newcontent])):
-700                raise Exception("m_eff is undefined at all timeslices")
-701            return np.arccosh(Corr(newcontent, padding=[1, 1]))
-702
-703        else:
-704            raise Exception('Unknown variant.')
+            
623    def m_eff(self, variant='log', guess=1.0):
+624        """Returns the effective mass of the correlator as correlator object
+625
+626        Parameters
+627        ----------
+628        variant : str
+629            log : uses the standard effective mass log(C(t) / C(t+1))
+630            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.
+631            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.
+632            See, e.g., arXiv:1205.5380
+633            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
+634            logsym: uses the symmetric effective mass log(C(t-1) / C(t+1))/2
+635        guess : float
+636            guess for the root finder, only relevant for the root variant
+637        """
+638        if self.N != 1:
+639            raise Exception('Correlator must be projected before getting m_eff')
+640        if variant == 'log':
+641            newcontent = []
+642            for t in range(self.T - 1):
+643                if ((self.content[t] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
+644                    newcontent.append(None)
+645                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
+646                    newcontent.append(None)
+647                else:
+648                    newcontent.append(self.content[t] / self.content[t + 1])
+649            if (all([x is None for x in newcontent])):
+650                raise Exception('m_eff is undefined at all timeslices')
+651
+652            return np.log(Corr(newcontent, padding=[0, 1]))
+653
+654        elif variant == 'logsym':
+655            newcontent = []
+656            for t in range(1, self.T - 1):
+657                if ((self.content[t - 1] is None) or (self.content[t + 1] is None)) or (self.content[t + 1][0].value == 0):
+658                    newcontent.append(None)
+659                elif self.content[t - 1][0].value / self.content[t + 1][0].value < 0:
+660                    newcontent.append(None)
+661                else:
+662                    newcontent.append(self.content[t - 1] / self.content[t + 1])
+663            if (all([x is None for x in newcontent])):
+664                raise Exception('m_eff is undefined at all timeslices')
+665
+666            return np.log(Corr(newcontent, padding=[1, 1])) / 2
+667
+668        elif variant in ['periodic', 'cosh', 'sinh']:
+669            if variant in ['periodic', 'cosh']:
+670                func = anp.cosh
+671            else:
+672                func = anp.sinh
+673
+674            def root_function(x, d):
+675                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
+676
+677            newcontent = []
+678            for t in range(self.T - 1):
+679                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 1][0].value == 0):
+680                    newcontent.append(None)
+681                # Fill the two timeslices in the middle of the lattice with their predecessors
+682                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+683                    newcontent.append(newcontent[-1])
+684                elif self.content[t][0].value / self.content[t + 1][0].value < 0:
+685                    newcontent.append(None)
+686                else:
+687                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+688            if (all([x is None for x in newcontent])):
+689                raise Exception('m_eff is undefined at all timeslices')
+690
+691            return Corr(newcontent, padding=[0, 1])
+692
+693        elif variant == 'arccosh':
+694            newcontent = []
+695            for t in range(1, self.T - 1):
+696                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):
+697                    newcontent.append(None)
+698                else:
+699                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+700            if (all([x is None for x in newcontent])):
+701                raise Exception("m_eff is undefined at all timeslices")
+702            return np.arccosh(Corr(newcontent, padding=[1, 1]))
+703
+704        else:
+705            raise Exception('Unknown variant.')
 

 
 
@@ -4021,39 +4024,39 @@ guess for the root finder, only relevant for the root variant
- 
706    def fit(self, function, fitrange=None, silent=False, **kwargs):
-707        r'''Fits function to the data
-708
-709        Parameters
-710        ----------
-711        function : obj
-712            function to fit to the data. See fits.least_squares for details.
-713        fitrange : list
-714            Two element list containing the timeslices on which the fit is supposed to start and stop.
-715            Caution: This range is inclusive as opposed to standard python indexing.
-716            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
-717            If not specified, self.prange or all timeslices are used.
-718        silent : bool
-719            Decides whether output is printed to the standard output.
-720        '''
-721        if self.N != 1:
-722            raise Exception("Correlator must be projected before fitting")
-723
-724        if fitrange is None:
-725            if self.prange:
-726                fitrange = self.prange
-727            else:
-728                fitrange = [0, self.T - 1]
-729        else:
-730            if not isinstance(fitrange, list):
-731                raise Exception("fitrange has to be a list with two elements")
-732            if len(fitrange) != 2:
-733                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
-734
-735        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-736        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-737        result = least_squares(xs, ys, function, silent=silent, **kwargs)
-738        return result
+            
707    def fit(self, function, fitrange=None, silent=False, **kwargs):
+708        r'''Fits function to the data
+709
+710        Parameters
+711        ----------
+712        function : obj
+713            function to fit to the data. See fits.least_squares for details.
+714        fitrange : list
+715            Two element list containing the timeslices on which the fit is supposed to start and stop.
+716            Caution: This range is inclusive as opposed to standard python indexing.
+717            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+718            If not specified, self.prange or all timeslices are used.
+719        silent : bool
+720            Decides whether output is printed to the standard output.
+721        '''
+722        if self.N != 1:
+723            raise Exception("Correlator must be projected before fitting")
+724
+725        if fitrange is None:
+726            if self.prange:
+727                fitrange = self.prange
+728            else:
+729                fitrange = [0, self.T - 1]
+730        else:
+731            if not isinstance(fitrange, list):
+732                raise Exception("fitrange has to be a list with two elements")
+733            if len(fitrange) != 2:
+734                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
+735
+736        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+737        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+738        result = least_squares(xs, ys, function, silent=silent, **kwargs)
+739        return result
 

 
 
@@ -4087,42 +4090,42 @@ Decides whether output is printed to the standard output.
- 
740    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
-741        """ Extract a plateau value from a Corr object
-742
-743        Parameters
-744        ----------
-745        plateau_range : list
-746            list with two entries, indicating the first and the last timeslice
-747            of the plateau region.
-748        method : str
-749            method to extract the plateau.
-750                'fit' fits a constant to the plateau region
-751                'avg', 'average' or 'mean' just average over the given timeslices.
-752        auto_gamma : bool
-753            apply gamma_method with default parameters to the Corr. Defaults to None
-754        """
-755        if not plateau_range:
-756            if self.prange:
-757                plateau_range = self.prange
-758            else:
-759                raise Exception("no plateau range provided")
-760        if self.N != 1:
-761            raise Exception("Correlator must be projected before getting a plateau.")
-762        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
-763            raise Exception("plateau is undefined at all timeslices in plateaurange.")
-764        if auto_gamma:
-765            self.gamma_method()
-766        if method == "fit":
-767            def const_func(a, t):
-768                return a[0]
-769            return self.fit(const_func, plateau_range)[0]
-770        elif method in ["avg", "average", "mean"]:
-771            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
-772            return returnvalue
-773
-774        else:
-775            raise Exception("Unsupported plateau method: " + method)
+            
741    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+742        """ Extract a plateau value from a Corr object
+743
+744        Parameters
+745        ----------
+746        plateau_range : list
+747            list with two entries, indicating the first and the last timeslice
+748            of the plateau region.
+749        method : str
+750            method to extract the plateau.
+751                'fit' fits a constant to the plateau region
+752                'avg', 'average' or 'mean' just average over the given timeslices.
+753        auto_gamma : bool
+754            apply gamma_method with default parameters to the Corr. Defaults to None
+755        """
+756        if not plateau_range:
+757            if self.prange:
+758                plateau_range = self.prange
+759            else:
+760                raise Exception("no plateau range provided")
+761        if self.N != 1:
+762            raise Exception("Correlator must be projected before getting a plateau.")
+763        if (all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+764            raise Exception("plateau is undefined at all timeslices in plateaurange.")
+765        if auto_gamma:
+766            self.gamma_method()
+767        if method == "fit":
+768            def const_func(a, t):
+769                return a[0]
+770            return self.fit(const_func, plateau_range)[0]
+771        elif method in ["avg", "average", "mean"]:
+772            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+773            return returnvalue
+774
+775        else:
+776            raise Exception("Unsupported plateau method: " + method)
 

 
 
@@ -4156,17 +4159,17 @@ apply gamma_method with default parameters to the Corr. Defaults to None
- 
777    def set_prange(self, prange):
-778        """Sets the attribute prange of the Corr object."""
-779        if not len(prange) == 2:
-780            raise Exception("prange must be a list or array with two values")
-781        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
-782            raise Exception("Start and end point must be integers")
-783        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
-784            raise Exception("Start and end point must define a range in the interval 0,T")
-785
-786        self.prange = prange
-787        return
+            
778    def set_prange(self, prange):
+779        """Sets the attribute prange of the Corr object."""
+780        if not len(prange) == 2:
+781            raise Exception("prange must be a list or array with two values")
+782        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+783            raise Exception("Start and end point must be integers")
+784        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
+785            raise Exception("Start and end point must define a range in the interval 0,T")
+786
+787        self.prange = prange
+788        return
 

 
 
@@ -4186,124 +4189,124 @@ apply gamma_method with default parameters to the Corr. Defaults to None
- 
789    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):
-790        """Plots the correlator using the tag of the correlator as label if available.
-791
-792        Parameters
-793        ----------
-794        x_range : list
-795            list of two values, determining the range of the x-axis e.g. [4, 8].
-796        comp : Corr or list of Corr
-797            Correlator or list of correlators which are plotted for comparison.
-798            The tags of these correlators are used as labels if available.
-799        logscale : bool
-800            Sets y-axis to logscale.
-801        plateau : Obs
-802            Plateau value to be visualized in the figure.
-803        fit_res : Fit_result
-804            Fit_result object to be visualized.
-805        ylabel : str
-806            Label for the y-axis.
-807        save : str
-808            path to file in which the figure should be saved.
-809        auto_gamma : bool
-810            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
-811        hide_sigma : float
-812            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
-813        references : list
-814            List of floating point values that are displayed as horizontal lines for reference.
-815        title : string
-816            Optional title of the figure.
-817        """
-818        if self.N != 1:
-819            raise Exception("Correlator must be projected before plotting")
-820
-821        if auto_gamma:
-822            self.gamma_method()
-823
-824        if x_range is None:
-825            x_range = [0, self.T - 1]
-826
-827        fig = plt.figure()
-828        ax1 = fig.add_subplot(111)
-829
-830        x, y, y_err = self.plottable()
-831        if hide_sigma:
-832            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
-833        else:
-834            hide_from = None
-835        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
-836        if logscale:
-837            ax1.set_yscale('log')
-838        else:
-839            if y_range is None:
-840                try:
-841                    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)])
-842                    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)])
-843                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
-844                except Exception:
-845                    pass
-846            else:
-847                ax1.set_ylim(y_range)
-848        if comp:
-849            if isinstance(comp, (Corr, list)):
-850                for corr in comp if isinstance(comp, list) else [comp]:
-851                    if auto_gamma:
-852                        corr.gamma_method()
-853                    x, y, y_err = corr.plottable()
-854                    if hide_sigma:
-855                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
-856                    else:
-857                        hide_from = None
-858                    ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
-859            else:
-860                raise Exception("'comp' must be a correlator or a list of correlators.")
-861
-862        if plateau:
-863            if isinstance(plateau, Obs):
-864                if auto_gamma:
-865                    plateau.gamma_method()
-866                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
-867                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
-868            else:
-869                raise Exception("'plateau' must be an Obs")
-870
-871        if references:
-872            if isinstance(references, list):
-873                for ref in references:
-874                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
-875            else:
-876                raise Exception("'references' must be a list of floating pint values.")
-877
-878        if self.prange:
-879            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
-880            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
-881
-882        if fit_res:
-883            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
-884            ax1.plot(x_samples,
-885                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
-886                     ls='-', marker=',', lw=2)
-887
-888        ax1.set_xlabel(r'$x_0 / a$')
-889        if ylabel:
-890            ax1.set_ylabel(ylabel)
-891        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
-892
-893        handles, labels = ax1.get_legend_handles_labels()
-894        if labels:
-895            ax1.legend()
-896
-897        if title:
-898            plt.title(title)
-899
-900        plt.draw()
-901
-902        if save:
-903            if isinstance(save, str):
-904                fig.savefig(save, bbox_inches='tight')
-905            else:
-906                raise Exception("'save' has to be a string.")
+            
790    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):
+791        """Plots the correlator using the tag of the correlator as label if available.
+792
+793        Parameters
+794        ----------
+795        x_range : list
+796            list of two values, determining the range of the x-axis e.g. [4, 8].
+797        comp : Corr or list of Corr
+798            Correlator or list of correlators which are plotted for comparison.
+799            The tags of these correlators are used as labels if available.
+800        logscale : bool
+801            Sets y-axis to logscale.
+802        plateau : Obs
+803            Plateau value to be visualized in the figure.
+804        fit_res : Fit_result
+805            Fit_result object to be visualized.
+806        ylabel : str
+807            Label for the y-axis.
+808        save : str
+809            path to file in which the figure should be saved.
+810        auto_gamma : bool
+811            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+812        hide_sigma : float
+813            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+814        references : list
+815            List of floating point values that are displayed as horizontal lines for reference.
+816        title : string
+817            Optional title of the figure.
+818        """
+819        if self.N != 1:
+820            raise Exception("Correlator must be projected before plotting")
+821
+822        if auto_gamma:
+823            self.gamma_method()
+824
+825        if x_range is None:
+826            x_range = [0, self.T - 1]
+827
+828        fig = plt.figure()
+829        ax1 = fig.add_subplot(111)
+830
+831        x, y, y_err = self.plottable()
+832        if hide_sigma:
+833            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+834        else:
+835            hide_from = None
+836        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+837        if logscale:
+838            ax1.set_yscale('log')
+839        else:
+840            if y_range is None:
+841                try:
+842                    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)])
+843                    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)])
+844                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+845                except Exception:
+846                    pass
+847            else:
+848                ax1.set_ylim(y_range)
+849        if comp:
+850            if isinstance(comp, (Corr, list)):
+851                for corr in comp if isinstance(comp, list) else [comp]:
+852                    if auto_gamma:
+853                        corr.gamma_method()
+854                    x, y, y_err = corr.plottable()
+855                    if hide_sigma:
+856                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+857                    else:
+858                        hide_from = None
+859                    ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+860            else:
+861                raise Exception("'comp' must be a correlator or a list of correlators.")
+862
+863        if plateau:
+864            if isinstance(plateau, Obs):
+865                if auto_gamma:
+866                    plateau.gamma_method()
+867                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+868                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+869            else:
+870                raise Exception("'plateau' must be an Obs")
+871
+872        if references:
+873            if isinstance(references, list):
+874                for ref in references:
+875                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+876            else:
+877                raise Exception("'references' must be a list of floating pint values.")
+878
+879        if self.prange:
+880            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
+881            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+882
+883        if fit_res:
+884            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+885            ax1.plot(x_samples,
+886                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
+887                     ls='-', marker=',', lw=2)
+888
+889        ax1.set_xlabel(r'$x_0 / a$')
+890        if ylabel:
+891            ax1.set_ylabel(ylabel)
+892        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
+893
+894        handles, labels = ax1.get_legend_handles_labels()
+895        if labels:
+896            ax1.legend()
+897
+898        if title:
+899            plt.title(title)
+900
+901        plt.draw()
+902
+903        if save:
+904            if isinstance(save, str):
+905                fig.savefig(save, bbox_inches='tight')
+906            else:
+907                raise Exception("'save' has to be a string.")
 

 
 
@@ -4351,34 +4354,34 @@ Optional title of the figure.
- 
908    def spaghetti_plot(self, logscale=True):
-909        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-910
-911        Parameters
-912        ----------
-913        logscale : bool
-914            Determines whether the scale of the y-axis is logarithmic or standard.
-915        """
-916        if self.N != 1:
-917            raise Exception("Correlator needs to be projected first.")
-918
-919        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]))
-920        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
-921
-922        for name in mc_names:
-923            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
-924
-925            fig = plt.figure()
-926            ax = fig.add_subplot(111)
-927            for dat in data:
-928                ax.plot(x0_vals, dat, ls='-', marker='')
-929
-930            if logscale is True:
-931                ax.set_yscale('log')
-932
-933            ax.set_xlabel(r'$x_0 / a$')
-934            plt.title(name)
-935            plt.draw()
+            
909    def spaghetti_plot(self, logscale=True):
+910        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+911
+912        Parameters
+913        ----------
+914        logscale : bool
+915            Determines whether the scale of the y-axis is logarithmic or standard.
+916        """
+917        if self.N != 1:
+918            raise Exception("Correlator needs to be projected first.")
+919
+920        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]))
+921        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
+922
+923        for name in mc_names:
+924            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+925
+926            fig = plt.figure()
+927            ax = fig.add_subplot(111)
+928            for dat in data:
+929                ax.plot(x0_vals, dat, ls='-', marker='')
+930
+931            if logscale is True:
+932                ax.set_yscale('log')
+933
+934            ax.set_xlabel(r'$x_0 / a$')
+935            plt.title(name)
+936            plt.draw()
 

 
 
@@ -4405,29 +4408,29 @@ Determines whether the scale of the y-axis is logarithmic or standard.
- 
937    def dump(self, filename, datatype="json.gz", **kwargs):
-938        """Dumps the Corr into a file of chosen type
-939        Parameters
-940        ----------
-941        filename : str
-942            Name of the file to be saved.
-943        datatype : str
-944            Format of the exported file. Supported formats include
-945            "json.gz" and "pickle"
-946        path : str
-947            specifies a custom path for the file (default '.')
-948        """
-949        if datatype == "json.gz":
-950            from .input.json import dump_to_json
-951            if 'path' in kwargs:
-952                file_name = kwargs.get('path') + '/' + filename
-953            else:
-954                file_name = filename
-955            dump_to_json(self, file_name)
-956        elif datatype == "pickle":
-957            dump_object(self, filename, **kwargs)
-958        else:
-959            raise Exception("Unknown datatype " + str(datatype))
+            
938    def dump(self, filename, datatype="json.gz", **kwargs):
+939        """Dumps the Corr into a file of chosen type
+940        Parameters
+941        ----------
+942        filename : str
+943            Name of the file to be saved.
+944        datatype : str
+945            Format of the exported file. Supported formats include
+946            "json.gz" and "pickle"
+947        path : str
+948            specifies a custom path for the file (default '.')
+949        """
+950        if datatype == "json.gz":
+951            from .input.json import dump_to_json
+952            if 'path' in kwargs:
+953                file_name = kwargs.get('path') + '/' + filename
+954            else:
+955                file_name = filename
+956            dump_to_json(self, file_name)
+957        elif datatype == "pickle":
+958            dump_object(self, filename, **kwargs)
+959        else:
+960            raise Exception("Unknown datatype " + str(datatype))
 

 
 
@@ -4459,8 +4462,8 @@ specifies a custom path for the file (default '.')
- 
961    def print(self, print_range=None):
-962        print(self.__repr__(print_range))
+            
962    def print(self, print_range=None):
+963        print(self.__repr__(print_range))
 

 
 
@@ -4478,8 +4481,8 @@ specifies a custom path for the file (default '.')
- 
1126    def sqrt(self):
-1127        return self ** 0.5
+            
1127    def sqrt(self):
+1128        return self ** 0.5
 

 
 
@@ -4497,9 +4500,9 @@ specifies a custom path for the file (default '.')
- 
1129    def log(self):
-1130        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
-1131        return Corr(newcontent, prange=self.prange)
+            
1130    def log(self):
+1131        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
+1132        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4517,9 +4520,9 @@ specifies a custom path for the file (default '.')
- 
1133    def exp(self):
-1134        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
-1135        return Corr(newcontent, prange=self.prange)
+            
1134    def exp(self):
+1135        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
+1136        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4537,8 +4540,8 @@ specifies a custom path for the file (default '.')
- 
1148    def sin(self):
-1149        return self._apply_func_to_corr(np.sin)
+            
1149    def sin(self):
+1150        return self._apply_func_to_corr(np.sin)
 

 
 
@@ -4556,8 +4559,8 @@ specifies a custom path for the file (default '.')
- 
1151    def cos(self):
-1152        return self._apply_func_to_corr(np.cos)
+            
1152    def cos(self):
+1153        return self._apply_func_to_corr(np.cos)
 

 
 
@@ -4575,8 +4578,8 @@ specifies a custom path for the file (default '.')
- 
1154    def tan(self):
-1155        return self._apply_func_to_corr(np.tan)
+            
1155    def tan(self):
+1156        return self._apply_func_to_corr(np.tan)
 

 
 
@@ -4594,8 +4597,8 @@ specifies a custom path for the file (default '.')
- 
1157    def sinh(self):
-1158        return self._apply_func_to_corr(np.sinh)
+            
1158    def sinh(self):
+1159        return self._apply_func_to_corr(np.sinh)
 

 
 
@@ -4613,8 +4616,8 @@ specifies a custom path for the file (default '.')
- 
1160    def cosh(self):
-1161        return self._apply_func_to_corr(np.cosh)
+            
1161    def cosh(self):
+1162        return self._apply_func_to_corr(np.cosh)
 

 
 
@@ -4632,8 +4635,8 @@ specifies a custom path for the file (default '.')
- 
1163    def tanh(self):
-1164        return self._apply_func_to_corr(np.tanh)
+            
1164    def tanh(self):
+1165        return self._apply_func_to_corr(np.tanh)
 

 
 
@@ -4651,8 +4654,8 @@ specifies a custom path for the file (default '.')
- 
1166    def arcsin(self):
-1167        return self._apply_func_to_corr(np.arcsin)
+            
1167    def arcsin(self):
+1168        return self._apply_func_to_corr(np.arcsin)
 

 
 
@@ -4670,8 +4673,8 @@ specifies a custom path for the file (default '.')
- 
1169    def arccos(self):
-1170        return self._apply_func_to_corr(np.arccos)
+            
1170    def arccos(self):
+1171        return self._apply_func_to_corr(np.arccos)
 

 
 
@@ -4689,8 +4692,8 @@ specifies a custom path for the file (default '.')
- 
1172    def arctan(self):
-1173        return self._apply_func_to_corr(np.arctan)
+            
1173    def arctan(self):
+1174        return self._apply_func_to_corr(np.arctan)
 

 
 
@@ -4708,8 +4711,8 @@ specifies a custom path for the file (default '.')
- 
1175    def arcsinh(self):
-1176        return self._apply_func_to_corr(np.arcsinh)
+            
1176    def arcsinh(self):
+1177        return self._apply_func_to_corr(np.arcsinh)
 

 
 
@@ -4727,8 +4730,8 @@ specifies a custom path for the file (default '.')
- 
1178    def arccosh(self):
-1179        return self._apply_func_to_corr(np.arccosh)
+            
1179    def arccosh(self):
+1180        return self._apply_func_to_corr(np.arccosh)
 

 
 
@@ -4746,8 +4749,8 @@ specifies a custom path for the file (default '.')
- 
1181    def arctanh(self):
-1182        return self._apply_func_to_corr(np.arctanh)
+            
1182    def arctanh(self):
+1183        return self._apply_func_to_corr(np.arctanh)
 

 
 
@@ -4765,62 +4768,62 @@ specifies a custom path for the file (default '.')
- 
1217    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1218        r''' Project large correlation matrix to lowest states
-1219
-1220        This method can be used to reduce the size of an (N x N) correlation matrix
-1221        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1222        is still small.
-1223
-1224        Parameters
-1225        ----------
-1226        Ntrunc: int
-1227            Rank of the target matrix.
-1228        tproj: int
-1229            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1230            The default value is 3.
-1231        t0proj: int
-1232            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1233            discouraged for O(a) improved theories, since the correctness of the procedure
-1234            cannot be granted in this case. The default value is 2.
-1235        basematrix : Corr
-1236            Correlation matrix that is used to determine the eigenvectors of the
-1237            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1238            is is not specified.
-1239
-1240        Notes
-1241        -----
-1242        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1243        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}$
-1244        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1245        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1246        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1247        correlation matrix and to remove some noise that is added by irrelevant operators.
-1248        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1249        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1250        '''
-1251
-1252        if self.N == 1:
-1253            raise Exception('Method cannot be applied to one-dimensional correlators.')
-1254        if basematrix is None:
-1255            basematrix = self
-1256        if Ntrunc >= basematrix.N:
-1257            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
-1258        if basematrix.N != self.N:
-1259            raise Exception('basematrix and targetmatrix have to be of the same size.')
-1260
-1261        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
-1262
-1263        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1264        rmat = []
-1265        for t in range(basematrix.T):
-1266            for i in range(Ntrunc):
-1267                for j in range(Ntrunc):
-1268                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1269            rmat.append(np.copy(tmpmat))
-1270
-1271        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
-1272        return Corr(newcontent)
+            
1218    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1219        r''' Project large correlation matrix to lowest states
+1220
+1221        This method can be used to reduce the size of an (N x N) correlation matrix
+1222        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1223        is still small.
+1224
+1225        Parameters
+1226        ----------
+1227        Ntrunc: int
+1228            Rank of the target matrix.
+1229        tproj: int
+1230            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1231            The default value is 3.
+1232        t0proj: int
+1233            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1234            discouraged for O(a) improved theories, since the correctness of the procedure
+1235            cannot be granted in this case. The default value is 2.
+1236        basematrix : Corr
+1237            Correlation matrix that is used to determine the eigenvectors of the
+1238            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1239            is is not specified.
+1240
+1241        Notes
+1242        -----
+1243        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1244        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}$
+1245        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1246        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1247        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1248        correlation matrix and to remove some noise that is added by irrelevant operators.
+1249        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1250        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1251        '''
+1252
+1253        if self.N == 1:
+1254            raise Exception('Method cannot be applied to one-dimensional correlators.')
+1255        if basematrix is None:
+1256            basematrix = self
+1257        if Ntrunc >= basematrix.N:
+1258            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
+1259        if basematrix.N != self.N:
+1260            raise Exception('basematrix and targetmatrix have to be of the same size.')
+1261
+1262        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
+1263
+1264        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1265        rmat = []
+1266        for t in range(basematrix.T):
+1267            for i in range(Ntrunc):
+1268                for j in range(Ntrunc):
+1269                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1270            rmat.append(np.copy(tmpmat))
+1271
+1272        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
+1273        return Corr(newcontent)