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

324 if sort == "Eigenvalue" and ts is not None: 325 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) 326 all_vecs = [None] * (t0 + 1) -327 for t in range(t0 + 1, self.T): -328 try: -329 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") -330 for i in range(self.N): -331 for j in range(self.N): -332 G0[i, j] = symmetric_corr[t0][i, j].value -333 Gt[i, j] = symmetric_corr[t][i, j].value -334 -335 all_vecs.append(_GEVP_solver(Gt, G0)) -336 except Exception: -337 all_vecs.append(None) -338 if sort == "Eigenvector": -339 if (ts is None): -340 raise Exception("ts is required for the Eigenvector sorting method.") -341 all_vecs = _sort_vectors(all_vecs, ts) +327 G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0]) +328 for t in range(t0 + 1, self.T): +329 try: +330 Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t]) +331 all_vecs.append(_GEVP_solver(Gt, G0)) +332 except Exception: +333 all_vecs.append(None) +334 if sort == "Eigenvector": +335 if (ts is None): +336 raise Exception("ts is required for the Eigenvector sorting method.") +337 all_vecs = _sort_vectors(all_vecs, ts) +338 +339 reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] +340 else: +341 raise Exception("Unkown value for 'sort'.") 342 -343 reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)] -344 else: -345 raise Exception("Unkown value for 'sort'.") -346 -347 if "state" in kwargs: -348 return reordered_vecs[kwargs.get("state")] -349 else: -350 return reordered_vecs +343 if "state" in kwargs: +344 return reordered_vecs[kwargs.get("state")] +345 else: +346 return reordered_vecs @@ -3387,18 +3375,18 @@ Returns only the vector(s) for a specified state. The lowest state is zero. - 
352    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
-353        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
-354
-355        Parameters
-356        ----------
-357        state : int
-358            The state one is interested in ordered by energy. The lowest state is zero.
-359
-360        All other parameters are identical to the ones of Corr.GEVP.
-361        """
-362        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
-363        return self.projected(vec)
+            
348    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
+349        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
+350
+351        Parameters
+352        ----------
+353        state : int
+354            The state one is interested in ordered by energy. The lowest state is zero.
+355
+356        All other parameters are identical to the ones of Corr.GEVP.
+357        """
+358        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
+359        return self.projected(vec)
 

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

 
 
@@ -3499,15 +3487,15 @@ determines whether the matrix is extended periodically
- 
406    def roll(self, dt):
-407        """Periodically shift the correlator by dt timeslices
-408
-409        Parameters
-410        ----------
-411        dt : int
-412            number of timeslices
-413        """
-414        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
+            
402    def roll(self, dt):
+403        """Periodically shift the correlator by dt timeslices
+404
+405        Parameters
+406        ----------
+407        dt : int
+408            number of timeslices
+409        """
+410        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
 

 
 
@@ -3534,9 +3522,9 @@ number of timeslices
- 
416    def reverse(self):
-417        """Reverse the time ordering of the Corr"""
-418        return Corr(self.content[:: -1])
+            
412    def reverse(self):
+413        """Reverse the time ordering of the Corr"""
+414        return Corr(self.content[:: -1])
 

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

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

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

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

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

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

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

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

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

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

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

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

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

 
 
@@ -4483,8 +4471,8 @@ specifies a custom path for the file (default '.')
- 
969    def print(self, print_range=None):
-970        print(self.__repr__(print_range))
+            
965    def print(self, print_range=None):
+966        print(self.__repr__(print_range))
 

 
 
@@ -4502,8 +4490,8 @@ specifies a custom path for the file (default '.')
- 
1134    def sqrt(self):
-1135        return self ** 0.5
+            
1130    def sqrt(self):
+1131        return self ** 0.5
 

 
 
@@ -4521,9 +4509,9 @@ specifies a custom path for the file (default '.')
- 
1137    def log(self):
-1138        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
-1139        return Corr(newcontent, prange=self.prange)
+            
1133    def log(self):
+1134        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
+1135        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4541,9 +4529,9 @@ specifies a custom path for the file (default '.')
- 
1141    def exp(self):
-1142        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
-1143        return Corr(newcontent, prange=self.prange)
+            
1137    def exp(self):
+1138        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
+1139        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4561,8 +4549,8 @@ specifies a custom path for the file (default '.')
- 
1156    def sin(self):
-1157        return self._apply_func_to_corr(np.sin)
+            
1152    def sin(self):
+1153        return self._apply_func_to_corr(np.sin)
 

 
 
@@ -4580,8 +4568,8 @@ specifies a custom path for the file (default '.')
- 
1159    def cos(self):
-1160        return self._apply_func_to_corr(np.cos)
+            
1155    def cos(self):
+1156        return self._apply_func_to_corr(np.cos)
 

 
 
@@ -4599,8 +4587,8 @@ specifies a custom path for the file (default '.')
- 
1162    def tan(self):
-1163        return self._apply_func_to_corr(np.tan)
+            
1158    def tan(self):
+1159        return self._apply_func_to_corr(np.tan)
 

 
 
@@ -4618,8 +4606,8 @@ specifies a custom path for the file (default '.')
- 
1165    def sinh(self):
-1166        return self._apply_func_to_corr(np.sinh)
+            
1161    def sinh(self):
+1162        return self._apply_func_to_corr(np.sinh)
 

 
 
@@ -4637,8 +4625,8 @@ specifies a custom path for the file (default '.')
- 
1168    def cosh(self):
-1169        return self._apply_func_to_corr(np.cosh)
+            
1164    def cosh(self):
+1165        return self._apply_func_to_corr(np.cosh)
 

 
 
@@ -4656,8 +4644,8 @@ specifies a custom path for the file (default '.')
- 
1171    def tanh(self):
-1172        return self._apply_func_to_corr(np.tanh)
+            
1167    def tanh(self):
+1168        return self._apply_func_to_corr(np.tanh)
 

 
 
@@ -4675,8 +4663,8 @@ specifies a custom path for the file (default '.')
- 
1174    def arcsin(self):
-1175        return self._apply_func_to_corr(np.arcsin)
+            
1170    def arcsin(self):
+1171        return self._apply_func_to_corr(np.arcsin)
 

 
 
@@ -4694,8 +4682,8 @@ specifies a custom path for the file (default '.')
- 
1177    def arccos(self):
-1178        return self._apply_func_to_corr(np.arccos)
+            
1173    def arccos(self):
+1174        return self._apply_func_to_corr(np.arccos)
 

 
 
@@ -4713,8 +4701,8 @@ specifies a custom path for the file (default '.')
- 
1180    def arctan(self):
-1181        return self._apply_func_to_corr(np.arctan)
+            
1176    def arctan(self):
+1177        return self._apply_func_to_corr(np.arctan)
 

 
 
@@ -4732,8 +4720,8 @@ specifies a custom path for the file (default '.')
- 
1183    def arcsinh(self):
-1184        return self._apply_func_to_corr(np.arcsinh)
+            
1179    def arcsinh(self):
+1180        return self._apply_func_to_corr(np.arcsinh)
 

 
 
@@ -4751,8 +4739,8 @@ specifies a custom path for the file (default '.')
- 
1186    def arccosh(self):
-1187        return self._apply_func_to_corr(np.arccosh)
+            
1182    def arccosh(self):
+1183        return self._apply_func_to_corr(np.arccosh)
 

 
 
@@ -4770,8 +4758,8 @@ specifies a custom path for the file (default '.')
- 
1189    def arctanh(self):
-1190        return self._apply_func_to_corr(np.arctanh)
+            
1185    def arctanh(self):
+1186        return self._apply_func_to_corr(np.arctanh)
 

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