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

274 sp_vecs = _GEVP_solver(Gt, G0) 275 sp_vec = sp_vecs[state] 276 return sp_vec -277 else: -278 -279 all_vecs = [] -280 for t in range(self.T): -281 try: -282 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") -283 for i in range(self.N): -284 for j in range(self.N): -285 G0[i, j] = self.content[t0][i, j].value -286 Gt[i, j] = self.content[t][i, j].value -287 -288 sp_vecs = _GEVP_solver(Gt, G0) -289 if sorted_list == "Eigenvalue": -290 sp_vec = sp_vecs[state] -291 all_vecs.append(sp_vec) -292 else: -293 all_vecs.append(sp_vecs) -294 except Exception: -295 all_vecs.append(None) -296 if sorted_list == "Eigenvector": -297 if (ts is None): -298 raise Exception("ts is required for the Eigenvector sorting method.") -299 all_vecs = _sort_vectors(all_vecs, ts) -300 all_vecs = [a[state] for a in all_vecs] -301 -302 return all_vecs +277 elif sorted_list in ["Eigenvalue", "Eigenvector"]: +278 all_vecs = [] +279 for t in range(self.T): +280 try: +281 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") +282 for i in range(self.N): +283 for j in range(self.N): +284 G0[i, j] = self.content[t0][i, j].value +285 Gt[i, j] = self.content[t][i, j].value +286 +287 sp_vecs = _GEVP_solver(Gt, G0) +288 if sorted_list == "Eigenvalue": +289 sp_vec = sp_vecs[state] +290 all_vecs.append(sp_vec) +291 else: +292 all_vecs.append(sp_vecs) +293 except Exception: +294 all_vecs.append(None) +295 if sorted_list == "Eigenvector": +296 if (ts is None): +297 raise Exception("ts is required for the Eigenvector sorting method.") +298 all_vecs = _sort_vectors(all_vecs, ts) +299 all_vecs = [a[state] for a in all_vecs] +300 else: +301 raise Exception("Unkown value for 'sorted_list'.") +302 +303 return all_vecs @@ -3117,26 +3120,26 @@ if this argument is set, a list of vectors (len=self.T) is returned. If it is le
View Source - 
304    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
-305        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
-306
-307        Parameters
-308        ----------
-309        t0 : int
-310            The time t0 for G(t)v= lambda G(t_0)v
-311        ts : int
-312            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
-313            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
-314        state : int
-315            The state one is interested in ordered by energy. The lowest state is zero.
-316        sorted_list : string
-317            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
-318             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
-319             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
-320                              The reference state is identified by its eigenvalue at t=ts
-321        """
-322        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
-323        return self.projected(vec)
+            
305    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
+306        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
+307
+308        Parameters
+309        ----------
+310        t0 : int
+311            The time t0 for G(t)v= lambda G(t_0)v
+312        ts : int
+313            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
+314            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
+315        state : int
+316            The state one is interested in ordered by energy. The lowest state is zero.
+317        sorted_list : string
+318            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
+319             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
+320             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
+321                              The reference state is identified by its eigenvalue at t=ts
+322        """
+323        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
+324        return self.projected(vec)
@@ -3173,46 +3176,46 @@ if this argument is set, a list of vectors (len=self.T) is returned. If it is le
View Source - 
325    def Hankel(self, N, periodic=False):
-326        """Constructs an NxN Hankel matrix
-327
-328        C(t) c(t+1) ... c(t+n-1)
-329        C(t+1) c(t+2) ... c(t+n)
-330        .................
-331        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
-332
-333        Parameters
-334        ----------
-335        N : int
-336            Dimension of the Hankel matrix
-337        periodic : bool, optional
-338            determines whether the matrix is extended periodically
-339        """
-340
-341        if self.N != 1:
-342            raise Exception("Multi-operator Prony not implemented!")
-343
-344        array = np.empty([N, N], dtype="object")
-345        new_content = []
-346        for t in range(self.T):
-347            new_content.append(array.copy())
-348
-349        def wrap(i):
-350            while i >= self.T:
-351                i -= self.T
-352            return i
-353
-354        for t in range(self.T):
-355            for i in range(N):
-356                for j in range(N):
-357                    if periodic:
-358                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
-359                    elif (t + i + j) >= self.T:
-360                        new_content[t] = None
-361                    else:
-362                        new_content[t][i, j] = self.content[t + i + j][0]
-363
-364        return Corr(new_content)
+            
326    def Hankel(self, N, periodic=False):
+327        """Constructs an NxN Hankel matrix
+328
+329        C(t) c(t+1) ... c(t+n-1)
+330        C(t+1) c(t+2) ... c(t+n)
+331        .................
+332        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
+333
+334        Parameters
+335        ----------
+336        N : int
+337            Dimension of the Hankel matrix
+338        periodic : bool, optional
+339            determines whether the matrix is extended periodically
+340        """
+341
+342        if self.N != 1:
+343            raise Exception("Multi-operator Prony not implemented!")
+344
+345        array = np.empty([N, N], dtype="object")
+346        new_content = []
+347        for t in range(self.T):
+348            new_content.append(array.copy())
+349
+350        def wrap(i):
+351            while i >= self.T:
+352                i -= self.T
+353            return i
+354
+355        for t in range(self.T):
+356            for i in range(N):
+357                for j in range(N):
+358                    if periodic:
+359                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
+360                    elif (t + i + j) >= self.T:
+361                        new_content[t] = None
+362                    else:
+363                        new_content[t][i, j] = self.content[t + i + j][0]
+364
+365        return Corr(new_content)
@@ -3246,15 +3249,15 @@ determines whether the matrix is extended periodically
View Source - 
366    def roll(self, dt):
-367        """Periodically shift the correlator by dt timeslices
-368
-369        Parameters
-370        ----------
-371        dt : int
-372            number of timeslices
-373        """
-374        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
+            
367    def roll(self, dt):
+368        """Periodically shift the correlator by dt timeslices
+369
+370        Parameters
+371        ----------
+372        dt : int
+373            number of timeslices
+374        """
+375        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
@@ -3281,9 +3284,9 @@ number of timeslices
View Source - 
376    def reverse(self):
-377        """Reverse the time ordering of the Corr"""
-378        return Corr(self.content[:: -1])
+            
377    def reverse(self):
+378        """Reverse the time ordering of the Corr"""
+379        return Corr(self.content[:: -1])
@@ -3303,23 +3306,23 @@ number of timeslices
View Source - 
380    def thin(self, spacing=2, offset=0):
-381        """Thin out a correlator to suppress correlations
-382
-383        Parameters
-384        ----------
-385        spacing : int
-386            Keep only every 'spacing'th entry of the correlator
-387        offset : int
-388            Offset the equal spacing
-389        """
-390        new_content = []
-391        for t in range(self.T):
-392            if (offset + t) % spacing != 0:
-393                new_content.append(None)
-394            else:
-395                new_content.append(self.content[t])
-396        return Corr(new_content)
+            
381    def thin(self, spacing=2, offset=0):
+382        """Thin out a correlator to suppress correlations
+383
+384        Parameters
+385        ----------
+386        spacing : int
+387            Keep only every 'spacing'th entry of the correlator
+388        offset : int
+389            Offset the equal spacing
+390        """
+391        new_content = []
+392        for t in range(self.T):
+393            if (offset + t) % spacing != 0:
+394                new_content.append(None)
+395            else:
+396                new_content.append(self.content[t])
+397        return Corr(new_content)
@@ -3348,32 +3351,32 @@ Offset the equal spacing
View Source - 
398    def correlate(self, partner):
-399        """Correlate the correlator with another correlator or Obs
-400
-401        Parameters
-402        ----------
-403        partner : Obs or Corr
-404            partner to correlate the correlator with.
-405            Can either be an Obs which is correlated with all entries of the
-406            correlator or a Corr of same length.
-407        """
-408        new_content = []
-409        for x0, t_slice in enumerate(self.content):
-410            if t_slice is None:
-411                new_content.append(None)
-412            else:
-413                if isinstance(partner, Corr):
-414                    if partner.content[x0] is None:
-415                        new_content.append(None)
-416                    else:
-417                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
-418                elif isinstance(partner, Obs):  # Should this include CObs?
-419                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
-420                else:
-421                    raise Exception("Can only correlate with an Obs or a Corr.")
-422
-423        return Corr(new_content)
+            
399    def correlate(self, partner):
+400        """Correlate the correlator with another correlator or Obs
+401
+402        Parameters
+403        ----------
+404        partner : Obs or Corr
+405            partner to correlate the correlator with.
+406            Can either be an Obs which is correlated with all entries of the
+407            correlator or a Corr of same length.
+408        """
+409        new_content = []
+410        for x0, t_slice in enumerate(self.content):
+411            if t_slice is None:
+412                new_content.append(None)
+413            else:
+414                if isinstance(partner, Corr):
+415                    if partner.content[x0] is None:
+416                        new_content.append(None)
+417                    else:
+418                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
+419                elif isinstance(partner, Obs):  # Should this include CObs?
+420                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
+421                else:
+422                    raise Exception("Can only correlate with an Obs or a Corr.")
+423
+424        return Corr(new_content)
@@ -3402,26 +3405,26 @@ correlator or a Corr of same length.
View Source - 
425    def reweight(self, weight, **kwargs):
-426        """Reweight the correlator.
-427
-428        Parameters
-429        ----------
-430        weight : Obs
-431            Reweighting factor. An Observable that has to be defined on a superset of the
-432            configurations in obs[i].idl for all i.
-433        all_configs : bool
-434            if True, the reweighted observables are normalized by the average of
-435            the reweighting factor on all configurations in weight.idl and not
-436            on the configurations in obs[i].idl.
-437        """
-438        new_content = []
-439        for t_slice in self.content:
-440            if t_slice is None:
-441                new_content.append(None)
-442            else:
-443                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
-444        return Corr(new_content)
+            
426    def reweight(self, weight, **kwargs):
+427        """Reweight the correlator.
+428
+429        Parameters
+430        ----------
+431        weight : Obs
+432            Reweighting factor. An Observable that has to be defined on a superset of the
+433            configurations in obs[i].idl for all i.
+434        all_configs : bool
+435            if True, the reweighted observables are normalized by the average of
+436            the reweighting factor on all configurations in weight.idl and not
+437            on the configurations in obs[i].idl.
+438        """
+439        new_content = []
+440        for t_slice in self.content:
+441            if t_slice is None:
+442                new_content.append(None)
+443            else:
+444                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
+445        return Corr(new_content)
@@ -3453,33 +3456,33 @@ on the configurations in obs[i].idl.
View Source - 
446    def T_symmetry(self, partner, parity=+1):
-447        """Return the time symmetry average of the correlator and its partner
-448
-449        Parameters
-450        ----------
-451        partner : Corr
-452            Time symmetry partner of the Corr
-453        partity : int
-454            Parity quantum number of the correlator, can be +1 or -1
-455        """
-456        if not isinstance(partner, Corr):
-457            raise Exception("T partner has to be a Corr object.")
-458        if parity not in [+1, -1]:
-459            raise Exception("Parity has to be +1 or -1.")
-460        T_partner = parity * partner.reverse()
-461
-462        t_slices = []
-463        test = (self - T_partner)
-464        test.gamma_method()
-465        for x0, t_slice in enumerate(test.content):
-466            if t_slice is not None:
-467                if not t_slice[0].is_zero_within_error(5):
-468                    t_slices.append(x0)
-469        if t_slices:
-470            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
-471
-472        return (self + T_partner) / 2
+            
447    def T_symmetry(self, partner, parity=+1):
+448        """Return the time symmetry average of the correlator and its partner
+449
+450        Parameters
+451        ----------
+452        partner : Corr
+453            Time symmetry partner of the Corr
+454        partity : int
+455            Parity quantum number of the correlator, can be +1 or -1
+456        """
+457        if not isinstance(partner, Corr):
+458            raise Exception("T partner has to be a Corr object.")
+459        if parity not in [+1, -1]:
+460            raise Exception("Parity has to be +1 or -1.")
+461        T_partner = parity * partner.reverse()
+462
+463        t_slices = []
+464        test = (self - T_partner)
+465        test.gamma_method()
+466        for x0, t_slice in enumerate(test.content):
+467            if t_slice is not None:
+468                if not t_slice[0].is_zero_within_error(5):
+469                    t_slices.append(x0)
+470        if t_slices:
+471            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
+472
+473        return (self + T_partner) / 2
@@ -3508,57 +3511,57 @@ Parity quantum number of the correlator, can be +1 or -1
View Source - 
474    def deriv(self, variant="symmetric"):
-475        """Return the first derivative of the correlator with respect to x0.
-476
-477        Parameters
-478        ----------
-479        variant : str
-480            decides which definition of the finite differences derivative is used.
-481            Available choice: symmetric, forward, backward, improved, default: symmetric
-482        """
-483        if variant == "symmetric":
-484            newcontent = []
-485            for t in range(1, self.T - 1):
-486                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-487                    newcontent.append(None)
-488                else:
-489                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
-490            if(all([x is None for x in newcontent])):
-491                raise Exception('Derivative is undefined at all timeslices')
-492            return Corr(newcontent, padding=[1, 1])
-493        elif variant == "forward":
-494            newcontent = []
-495            for t in range(self.T - 1):
-496                if (self.content[t] is None) or (self.content[t + 1] is None):
-497                    newcontent.append(None)
-498                else:
-499                    newcontent.append(self.content[t + 1] - self.content[t])
-500            if(all([x is None for x in newcontent])):
-501                raise Exception("Derivative is undefined at all timeslices")
-502            return Corr(newcontent, padding=[0, 1])
-503        elif variant == "backward":
-504            newcontent = []
-505            for t in range(1, self.T):
-506                if (self.content[t - 1] is None) or (self.content[t] is None):
-507                    newcontent.append(None)
-508                else:
-509                    newcontent.append(self.content[t] - self.content[t - 1])
-510            if(all([x is None for x in newcontent])):
-511                raise Exception("Derivative is undefined at all timeslices")
-512            return Corr(newcontent, padding=[1, 0])
-513        elif variant == "improved":
-514            newcontent = []
-515            for t in range(2, self.T - 2):
-516                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):
-517                    newcontent.append(None)
-518                else:
-519                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
-520            if(all([x is None for x in newcontent])):
-521                raise Exception('Derivative is undefined at all timeslices')
-522            return Corr(newcontent, padding=[2, 2])
-523        else:
-524            raise Exception("Unknown variant.")
+            
475    def deriv(self, variant="symmetric"):
+476        """Return the first derivative of the correlator with respect to x0.
+477
+478        Parameters
+479        ----------
+480        variant : str
+481            decides which definition of the finite differences derivative is used.
+482            Available choice: symmetric, forward, backward, improved, default: symmetric
+483        """
+484        if variant == "symmetric":
+485            newcontent = []
+486            for t in range(1, self.T - 1):
+487                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+488                    newcontent.append(None)
+489                else:
+490                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
+491            if(all([x is None for x in newcontent])):
+492                raise Exception('Derivative is undefined at all timeslices')
+493            return Corr(newcontent, padding=[1, 1])
+494        elif variant == "forward":
+495            newcontent = []
+496            for t in range(self.T - 1):
+497                if (self.content[t] is None) or (self.content[t + 1] is None):
+498                    newcontent.append(None)
+499                else:
+500                    newcontent.append(self.content[t + 1] - self.content[t])
+501            if(all([x is None for x in newcontent])):
+502                raise Exception("Derivative is undefined at all timeslices")
+503            return Corr(newcontent, padding=[0, 1])
+504        elif variant == "backward":
+505            newcontent = []
+506            for t in range(1, self.T):
+507                if (self.content[t - 1] is None) or (self.content[t] is None):
+508                    newcontent.append(None)
+509                else:
+510                    newcontent.append(self.content[t] - self.content[t - 1])
+511            if(all([x is None for x in newcontent])):
+512                raise Exception("Derivative is undefined at all timeslices")
+513            return Corr(newcontent, padding=[1, 0])
+514        elif variant == "improved":
+515            newcontent = []
+516            for t in range(2, self.T - 2):
+517                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):
+518                    newcontent.append(None)
+519                else:
+520                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
+521            if(all([x is None for x in newcontent])):
+522                raise Exception('Derivative is undefined at all timeslices')
+523            return Corr(newcontent, padding=[2, 2])
+524        else:
+525            raise Exception("Unknown variant.")
@@ -3586,37 +3589,37 @@ Available choice: symmetric, forward, backward, improved, default: symmetric View Source - 
526    def second_deriv(self, variant="symmetric"):
-527        """Return the second derivative of the correlator with respect to x0.
-528
-529        Parameters
-530        ----------
-531        variant : str
-532            decides which definition of the finite differences derivative is used.
-533            Available choice: symmetric, improved, default: symmetric
-534        """
-535        if variant == "symmetric":
-536            newcontent = []
-537            for t in range(1, self.T - 1):
-538                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-539                    newcontent.append(None)
-540                else:
-541                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
-542            if(all([x is None for x in newcontent])):
-543                raise Exception("Derivative is undefined at all timeslices")
-544            return Corr(newcontent, padding=[1, 1])
-545        elif variant == "improved":
-546            newcontent = []
-547            for t in range(2, self.T - 2):
-548                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):
-549                    newcontent.append(None)
-550                else:
-551                    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]))
-552            if(all([x is None for x in newcontent])):
-553                raise Exception("Derivative is undefined at all timeslices")
-554            return Corr(newcontent, padding=[2, 2])
-555        else:
-556            raise Exception("Unknown variant.")
+            
527    def second_deriv(self, variant="symmetric"):
+528        """Return the second derivative of the correlator with respect to x0.
+529
+530        Parameters
+531        ----------
+532        variant : str
+533            decides which definition of the finite differences derivative is used.
+534            Available choice: symmetric, improved, default: symmetric
+535        """
+536        if variant == "symmetric":
+537            newcontent = []
+538            for t in range(1, self.T - 1):
+539                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+540                    newcontent.append(None)
+541                else:
+542                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
+543            if(all([x is None for x in newcontent])):
+544                raise Exception("Derivative is undefined at all timeslices")
+545            return Corr(newcontent, padding=[1, 1])
+546        elif variant == "improved":
+547            newcontent = []
+548            for t in range(2, self.T - 2):
+549                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):
+550                    newcontent.append(None)
+551                else:
+552                    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]))
+553            if(all([x is None for x in newcontent])):
+554                raise Exception("Derivative is undefined at all timeslices")
+555            return Corr(newcontent, padding=[2, 2])
+556        else:
+557            raise Exception("Unknown variant.")
 

 
         
@@ -3644,70 +3647,70 @@ Available choice: symmetric, improved, default: symmetric
 
             

             View Source
-            
558    def m_eff(self, variant='log', guess=1.0):
-559        """Returns the effective mass of the correlator as correlator object
-560
-561        Parameters
-562        ----------
-563        variant : str
-564            log : uses the standard effective mass log(C(t) / C(t+1))
-565            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.
-566            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.
-567            See, e.g., arXiv:1205.5380
-568            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
-569        guess : float
-570            guess for the root finder, only relevant for the root variant
-571        """
-572        if self.N != 1:
-573            raise Exception('Correlator must be projected before getting m_eff')
-574        if variant == 'log':
-575            newcontent = []
-576            for t in range(self.T - 1):
-577                if (self.content[t] is None) or (self.content[t + 1] is None):
-578                    newcontent.append(None)
-579                else:
-580                    newcontent.append(self.content[t] / self.content[t + 1])
-581            if(all([x is None for x in newcontent])):
-582                raise Exception('m_eff is undefined at all timeslices')
-583
-584            return np.log(Corr(newcontent, padding=[0, 1]))
-585
-586        elif variant in ['periodic', 'cosh', 'sinh']:
-587            if variant in ['periodic', 'cosh']:
-588                func = anp.cosh
-589            else:
-590                func = anp.sinh
-591
-592            def root_function(x, d):
-593                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
-594
-595            newcontent = []
-596            for t in range(self.T - 1):
-597                if (self.content[t] is None) or (self.content[t + 1] is None):
-598                    newcontent.append(None)
-599                # Fill the two timeslices in the middle of the lattice with their predecessors
-600                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
-601                    newcontent.append(newcontent[-1])
-602                else:
-603                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
-604            if(all([x is None for x in newcontent])):
-605                raise Exception('m_eff is undefined at all timeslices')
-606
-607            return Corr(newcontent, padding=[0, 1])
-608
-609        elif variant == 'arccosh':
-610            newcontent = []
-611            for t in range(1, self.T - 1):
-612                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
-613                    newcontent.append(None)
-614                else:
-615                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
-616            if(all([x is None for x in newcontent])):
-617                raise Exception("m_eff is undefined at all timeslices")
-618            return np.arccosh(Corr(newcontent, padding=[1, 1]))
-619
-620        else:
-621            raise Exception('Unknown variant.')
+            
559    def m_eff(self, variant='log', guess=1.0):
+560        """Returns the effective mass of the correlator as correlator object
+561
+562        Parameters
+563        ----------
+564        variant : str
+565            log : uses the standard effective mass log(C(t) / C(t+1))
+566            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.
+567            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.
+568            See, e.g., arXiv:1205.5380
+569            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
+570        guess : float
+571            guess for the root finder, only relevant for the root variant
+572        """
+573        if self.N != 1:
+574            raise Exception('Correlator must be projected before getting m_eff')
+575        if variant == 'log':
+576            newcontent = []
+577            for t in range(self.T - 1):
+578                if (self.content[t] is None) or (self.content[t + 1] is None):
+579                    newcontent.append(None)
+580                else:
+581                    newcontent.append(self.content[t] / self.content[t + 1])
+582            if(all([x is None for x in newcontent])):
+583                raise Exception('m_eff is undefined at all timeslices')
+584
+585            return np.log(Corr(newcontent, padding=[0, 1]))
+586
+587        elif variant in ['periodic', 'cosh', 'sinh']:
+588            if variant in ['periodic', 'cosh']:
+589                func = anp.cosh
+590            else:
+591                func = anp.sinh
+592
+593            def root_function(x, d):
+594                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
+595
+596            newcontent = []
+597            for t in range(self.T - 1):
+598                if (self.content[t] is None) or (self.content[t + 1] is None):
+599                    newcontent.append(None)
+600                # Fill the two timeslices in the middle of the lattice with their predecessors
+601                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+602                    newcontent.append(newcontent[-1])
+603                else:
+604                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+605            if(all([x is None for x in newcontent])):
+606                raise Exception('m_eff is undefined at all timeslices')
+607
+608            return Corr(newcontent, padding=[0, 1])
+609
+610        elif variant == 'arccosh':
+611            newcontent = []
+612            for t in range(1, self.T - 1):
+613                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
+614                    newcontent.append(None)
+615                else:
+616                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+617            if(all([x is None for x in newcontent])):
+618                raise Exception("m_eff is undefined at all timeslices")
+619            return np.arccosh(Corr(newcontent, padding=[1, 1]))
+620
+621        else:
+622            raise Exception('Unknown variant.')
 

 
         

@@ -3740,39 +3743,39 @@ guess for the root finder, only relevant for the root variant
 
             

             View Source
-            
623    def fit(self, function, fitrange=None, silent=False, **kwargs):
-624        r'''Fits function to the data
-625
-626        Parameters
-627        ----------
-628        function : obj
-629            function to fit to the data. See fits.least_squares for details.
-630        fitrange : list
-631            Two element list containing the timeslices on which the fit is supposed to start and stop.
-632            Caution: This range is inclusive as opposed to standard python indexing.
-633            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
-634            If not specified, self.prange or all timeslices are used.
-635        silent : bool
-636            Decides whether output is printed to the standard output.
-637        '''
-638        if self.N != 1:
-639            raise Exception("Correlator must be projected before fitting")
-640
-641        if fitrange is None:
-642            if self.prange:
-643                fitrange = self.prange
-644            else:
-645                fitrange = [0, self.T - 1]
-646        else:
-647            if not isinstance(fitrange, list):
-648                raise Exception("fitrange has to be a list with two elements")
-649            if len(fitrange) != 2:
-650                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
-651
-652        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-653        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-654        result = least_squares(xs, ys, function, silent=silent, **kwargs)
-655        return result
+            
624    def fit(self, function, fitrange=None, silent=False, **kwargs):
+625        r'''Fits function to the data
+626
+627        Parameters
+628        ----------
+629        function : obj
+630            function to fit to the data. See fits.least_squares for details.
+631        fitrange : list
+632            Two element list containing the timeslices on which the fit is supposed to start and stop.
+633            Caution: This range is inclusive as opposed to standard python indexing.
+634            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+635            If not specified, self.prange or all timeslices are used.
+636        silent : bool
+637            Decides whether output is printed to the standard output.
+638        '''
+639        if self.N != 1:
+640            raise Exception("Correlator must be projected before fitting")
+641
+642        if fitrange is None:
+643            if self.prange:
+644                fitrange = self.prange
+645            else:
+646                fitrange = [0, self.T - 1]
+647        else:
+648            if not isinstance(fitrange, list):
+649                raise Exception("fitrange has to be a list with two elements")
+650            if len(fitrange) != 2:
+651                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
+652
+653        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+654        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+655        result = least_squares(xs, ys, function, silent=silent, **kwargs)
+656        return result
 

 
         

@@ -3806,42 +3809,42 @@ Decides whether output is printed to the standard output.
 
             

             View Source
-            
657    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
-658        """ Extract a plateau value from a Corr object
-659
-660        Parameters
-661        ----------
-662        plateau_range : list
-663            list with two entries, indicating the first and the last timeslice
-664            of the plateau region.
-665        method : str
-666            method to extract the plateau.
-667                'fit' fits a constant to the plateau region
-668                'avg', 'average' or 'mean' just average over the given timeslices.
-669        auto_gamma : bool
-670            apply gamma_method with default parameters to the Corr. Defaults to None
-671        """
-672        if not plateau_range:
-673            if self.prange:
-674                plateau_range = self.prange
-675            else:
-676                raise Exception("no plateau range provided")
-677        if self.N != 1:
-678            raise Exception("Correlator must be projected before getting a plateau.")
-679        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
-680            raise Exception("plateau is undefined at all timeslices in plateaurange.")
-681        if auto_gamma:
-682            self.gamma_method()
-683        if method == "fit":
-684            def const_func(a, t):
-685                return a[0]
-686            return self.fit(const_func, plateau_range)[0]
-687        elif method in ["avg", "average", "mean"]:
-688            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
-689            return returnvalue
-690
-691        else:
-692            raise Exception("Unsupported plateau method: " + method)
+            
658    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+659        """ Extract a plateau value from a Corr object
+660
+661        Parameters
+662        ----------
+663        plateau_range : list
+664            list with two entries, indicating the first and the last timeslice
+665            of the plateau region.
+666        method : str
+667            method to extract the plateau.
+668                'fit' fits a constant to the plateau region
+669                'avg', 'average' or 'mean' just average over the given timeslices.
+670        auto_gamma : bool
+671            apply gamma_method with default parameters to the Corr. Defaults to None
+672        """
+673        if not plateau_range:
+674            if self.prange:
+675                plateau_range = self.prange
+676            else:
+677                raise Exception("no plateau range provided")
+678        if self.N != 1:
+679            raise Exception("Correlator must be projected before getting a plateau.")
+680        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+681            raise Exception("plateau is undefined at all timeslices in plateaurange.")
+682        if auto_gamma:
+683            self.gamma_method()
+684        if method == "fit":
+685            def const_func(a, t):
+686                return a[0]
+687            return self.fit(const_func, plateau_range)[0]
+688        elif method in ["avg", "average", "mean"]:
+689            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+690            return returnvalue
+691
+692        else:
+693            raise Exception("Unsupported plateau method: " + method)
 

 
         

@@ -3875,17 +3878,17 @@ apply gamma_method with default parameters to the Corr. Defaults to None
 
             

             View Source
-            
694    def set_prange(self, prange):
-695        """Sets the attribute prange of the Corr object."""
-696        if not len(prange) == 2:
-697            raise Exception("prange must be a list or array with two values")
-698        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
-699            raise Exception("Start and end point must be integers")
-700        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
-701            raise Exception("Start and end point must define a range in the interval 0,T")
-702
-703        self.prange = prange
-704        return
+            
695    def set_prange(self, prange):
+696        """Sets the attribute prange of the Corr object."""
+697        if not len(prange) == 2:
+698            raise Exception("prange must be a list or array with two values")
+699        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+700            raise Exception("Start and end point must be integers")
+701        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
+702            raise Exception("Start and end point must define a range in the interval 0,T")
+703
+704        self.prange = prange
+705        return
 

 
         

@@ -3918,118 +3921,118 @@ apply gamma_method with default parameters to the Corr. Defaults to None
 
             

             View Source
-            
706    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):
-707        """Plots the correlator using the tag of the correlator as label if available.
-708
-709        Parameters
-710        ----------
-711        x_range : list
-712            list of two values, determining the range of the x-axis e.g. [4, 8]
-713        comp : Corr or list of Corr
-714            Correlator or list of correlators which are plotted for comparison.
-715            The tags of these correlators are used as labels if available.
-716        logscale : bool
-717            Sets y-axis to logscale
-718        plateau : Obs
-719            Plateau value to be visualized in the figure
-720        fit_res : Fit_result
-721            Fit_result object to be visualized
-722        ylabel : str
-723            Label for the y-axis
-724        save : str
-725            path to file in which the figure should be saved
-726        auto_gamma : bool
-727            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
-728        hide_sigma : float
-729            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
-730        references : list
-731            List of floating point values that are displayed as horizontal lines for reference.
-732        """
-733        if self.N != 1:
-734            raise Exception("Correlator must be projected before plotting")
-735
-736        if auto_gamma:
-737            self.gamma_method()
-738
-739        if x_range is None:
-740            x_range = [0, self.T - 1]
-741
-742        fig = plt.figure()
-743        ax1 = fig.add_subplot(111)
-744
-745        x, y, y_err = self.plottable()
-746        if hide_sigma:
-747            hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
-748        else:
-749            hide_from = None
-750        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
-751        if logscale:
-752            ax1.set_yscale('log')
-753        else:
-754            if y_range is None:
-755                try:
-756                    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)])
-757                    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)])
-758                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
-759                except Exception:
-760                    pass
-761            else:
-762                ax1.set_ylim(y_range)
-763        if comp:
-764            if isinstance(comp, (Corr, list)):
-765                for corr in comp if isinstance(comp, list) else [comp]:
-766                    if auto_gamma:
-767                        corr.gamma_method()
-768                    x, y, y_err = corr.plottable()
-769                    if hide_sigma:
-770                        hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
-771                    else:
-772                        hide_from = None
-773                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
-774            else:
-775                raise Exception("'comp' must be a correlator or a list of correlators.")
-776
-777        if plateau:
-778            if isinstance(plateau, Obs):
-779                if auto_gamma:
-780                    plateau.gamma_method()
-781                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
-782                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
-783            else:
-784                raise Exception("'plateau' must be an Obs")
-785
-786        if references:
-787            if isinstance(references, list):
-788                for ref in references:
-789                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
-790            else:
-791                raise Exception("'references' must be a list of floating pint values.")
-792
-793        if self.prange:
-794            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
-795            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
-796
-797        if fit_res:
-798            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
-799            ax1.plot(x_samples,
-800                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
-801                     ls='-', marker=',', lw=2)
-802
-803        ax1.set_xlabel(r'$x_0 / a$')
-804        if ylabel:
-805            ax1.set_ylabel(ylabel)
-806        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
-807
-808        handles, labels = ax1.get_legend_handles_labels()
-809        if labels:
-810            ax1.legend()
-811        plt.draw()
-812
-813        if save:
-814            if isinstance(save, str):
-815                fig.savefig(save)
-816            else:
-817                raise Exception("'save' has to be a string.")
+            
707    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):
+708        """Plots the correlator using the tag of the correlator as label if available.
+709
+710        Parameters
+711        ----------
+712        x_range : list
+713            list of two values, determining the range of the x-axis e.g. [4, 8]
+714        comp : Corr or list of Corr
+715            Correlator or list of correlators which are plotted for comparison.
+716            The tags of these correlators are used as labels if available.
+717        logscale : bool
+718            Sets y-axis to logscale
+719        plateau : Obs
+720            Plateau value to be visualized in the figure
+721        fit_res : Fit_result
+722            Fit_result object to be visualized
+723        ylabel : str
+724            Label for the y-axis
+725        save : str
+726            path to file in which the figure should be saved
+727        auto_gamma : bool
+728            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+729        hide_sigma : float
+730            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+731        references : list
+732            List of floating point values that are displayed as horizontal lines for reference.
+733        """
+734        if self.N != 1:
+735            raise Exception("Correlator must be projected before plotting")
+736
+737        if auto_gamma:
+738            self.gamma_method()
+739
+740        if x_range is None:
+741            x_range = [0, self.T - 1]
+742
+743        fig = plt.figure()
+744        ax1 = fig.add_subplot(111)
+745
+746        x, y, y_err = self.plottable()
+747        if hide_sigma:
+748            hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
+749        else:
+750            hide_from = None
+751        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+752        if logscale:
+753            ax1.set_yscale('log')
+754        else:
+755            if y_range is None:
+756                try:
+757                    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)])
+758                    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)])
+759                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+760                except Exception:
+761                    pass
+762            else:
+763                ax1.set_ylim(y_range)
+764        if comp:
+765            if isinstance(comp, (Corr, list)):
+766                for corr in comp if isinstance(comp, list) else [comp]:
+767                    if auto_gamma:
+768                        corr.gamma_method()
+769                    x, y, y_err = corr.plottable()
+770                    if hide_sigma:
+771                        hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
+772                    else:
+773                        hide_from = None
+774                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+775            else:
+776                raise Exception("'comp' must be a correlator or a list of correlators.")
+777
+778        if plateau:
+779            if isinstance(plateau, Obs):
+780                if auto_gamma:
+781                    plateau.gamma_method()
+782                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+783                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+784            else:
+785                raise Exception("'plateau' must be an Obs")
+786
+787        if references:
+788            if isinstance(references, list):
+789                for ref in references:
+790                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+791            else:
+792                raise Exception("'references' must be a list of floating pint values.")
+793
+794        if self.prange:
+795            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
+796            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+797
+798        if fit_res:
+799            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+800            ax1.plot(x_samples,
+801                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
+802                     ls='-', marker=',', lw=2)
+803
+804        ax1.set_xlabel(r'$x_0 / a$')
+805        if ylabel:
+806            ax1.set_ylabel(ylabel)
+807        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
+808
+809        handles, labels = ax1.get_legend_handles_labels()
+810        if labels:
+811            ax1.legend()
+812        plt.draw()
+813
+814        if save:
+815            if isinstance(save, str):
+816                fig.savefig(save)
+817            else:
+818                raise Exception("'save' has to be a string.")
 

 
         

@@ -4075,34 +4078,34 @@ List of floating point values that are displayed as horizontal lines for referen
 
             

             View Source
-            
819    def spaghetti_plot(self, logscale=True):
-820        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-821
-822        Parameters
-823        ----------
-824        logscale : bool
-825            Determines whether the scale of the y-axis is logarithmic or standard.
-826        """
-827        if self.N != 1:
-828            raise Exception("Correlator needs to be projected first.")
-829
-830        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]))
-831        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
-832
-833        for name in mc_names:
-834            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
-835
-836            fig = plt.figure()
-837            ax = fig.add_subplot(111)
-838            for dat in data:
-839                ax.plot(x0_vals, dat, ls='-', marker='')
-840
-841            if logscale is True:
-842                ax.set_yscale('log')
-843
-844            ax.set_xlabel(r'$x_0 / a$')
-845            plt.title(name)
-846            plt.draw()
+            
820    def spaghetti_plot(self, logscale=True):
+821        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+822
+823        Parameters
+824        ----------
+825        logscale : bool
+826            Determines whether the scale of the y-axis is logarithmic or standard.
+827        """
+828        if self.N != 1:
+829            raise Exception("Correlator needs to be projected first.")
+830
+831        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]))
+832        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
+833
+834        for name in mc_names:
+835            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+836
+837            fig = plt.figure()
+838            ax = fig.add_subplot(111)
+839            for dat in data:
+840                ax.plot(x0_vals, dat, ls='-', marker='')
+841
+842            if logscale is True:
+843                ax.set_yscale('log')
+844
+845            ax.set_xlabel(r'$x_0 / a$')
+846            plt.title(name)
+847            plt.draw()
 

 
         

@@ -4129,29 +4132,29 @@ Determines whether the scale of the y-axis is logarithmic or standard.
 
             

             View Source
-            
848    def dump(self, filename, datatype="json.gz", **kwargs):
-849        """Dumps the Corr into a file of chosen type
-850        Parameters
-851        ----------
-852        filename : str
-853            Name of the file to be saved.
-854        datatype : str
-855            Format of the exported file. Supported formats include
-856            "json.gz" and "pickle"
-857        path : str
-858            specifies a custom path for the file (default '.')
-859        """
-860        if datatype == "json.gz":
-861            from .input.json import dump_to_json
-862            if 'path' in kwargs:
-863                file_name = kwargs.get('path') + '/' + filename
-864            else:
-865                file_name = filename
-866            dump_to_json(self, file_name)
-867        elif datatype == "pickle":
-868            dump_object(self, filename, **kwargs)
-869        else:
-870            raise Exception("Unknown datatype " + str(datatype))
+            
849    def dump(self, filename, datatype="json.gz", **kwargs):
+850        """Dumps the Corr into a file of chosen type
+851        Parameters
+852        ----------
+853        filename : str
+854            Name of the file to be saved.
+855        datatype : str
+856            Format of the exported file. Supported formats include
+857            "json.gz" and "pickle"
+858        path : str
+859            specifies a custom path for the file (default '.')
+860        """
+861        if datatype == "json.gz":
+862            from .input.json import dump_to_json
+863            if 'path' in kwargs:
+864                file_name = kwargs.get('path') + '/' + filename
+865            else:
+866                file_name = filename
+867            dump_to_json(self, file_name)
+868        elif datatype == "pickle":
+869            dump_object(self, filename, **kwargs)
+870        else:
+871            raise Exception("Unknown datatype " + str(datatype))
 

 
         

@@ -4183,8 +4186,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
872    def print(self, range=[0, None]):
-873        print(self.__repr__(range))
+            
873    def print(self, range=[0, None]):
+874        print(self.__repr__(range))
 

 
         

@@ -4202,8 +4205,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1035    def sqrt(self):
-1036        return self**0.5
+            
1036    def sqrt(self):
+1037        return self**0.5
 

 
         

@@ -4221,9 +4224,9 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1038    def log(self):
-1039        newcontent = [None if (item is None) else np.log(item) for item in self.content]
-1040        return Corr(newcontent, prange=self.prange)
+            
1039    def log(self):
+1040        newcontent = [None if (item is None) else np.log(item) for item in self.content]
+1041        return Corr(newcontent, prange=self.prange)
 

 
         

@@ -4241,9 +4244,9 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1042    def exp(self):
-1043        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
-1044        return Corr(newcontent, prange=self.prange)
+            
1043    def exp(self):
+1044        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
+1045        return Corr(newcontent, prange=self.prange)
 

 
         

@@ -4261,8 +4264,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1057    def sin(self):
-1058        return self._apply_func_to_corr(np.sin)
+            
1058    def sin(self):
+1059        return self._apply_func_to_corr(np.sin)
 

 
         

@@ -4280,8 +4283,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1060    def cos(self):
-1061        return self._apply_func_to_corr(np.cos)
+            
1061    def cos(self):
+1062        return self._apply_func_to_corr(np.cos)
 

 
         

@@ -4299,8 +4302,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1063    def tan(self):
-1064        return self._apply_func_to_corr(np.tan)
+            
1064    def tan(self):
+1065        return self._apply_func_to_corr(np.tan)
 

 
         

@@ -4318,8 +4321,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1066    def sinh(self):
-1067        return self._apply_func_to_corr(np.sinh)
+            
1067    def sinh(self):
+1068        return self._apply_func_to_corr(np.sinh)
 

 
         

@@ -4337,8 +4340,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1069    def cosh(self):
-1070        return self._apply_func_to_corr(np.cosh)
+            
1070    def cosh(self):
+1071        return self._apply_func_to_corr(np.cosh)
 

 
         

@@ -4356,8 +4359,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1072    def tanh(self):
-1073        return self._apply_func_to_corr(np.tanh)
+            
1073    def tanh(self):
+1074        return self._apply_func_to_corr(np.tanh)
 

 
         

@@ -4375,8 +4378,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1075    def arcsin(self):
-1076        return self._apply_func_to_corr(np.arcsin)
+            
1076    def arcsin(self):
+1077        return self._apply_func_to_corr(np.arcsin)
 

 
         

@@ -4394,8 +4397,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1078    def arccos(self):
-1079        return self._apply_func_to_corr(np.arccos)
+            
1079    def arccos(self):
+1080        return self._apply_func_to_corr(np.arccos)
 

 
         

@@ -4413,8 +4416,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1081    def arctan(self):
-1082        return self._apply_func_to_corr(np.arctan)
+            
1082    def arctan(self):
+1083        return self._apply_func_to_corr(np.arctan)
 

 
         

@@ -4432,8 +4435,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1084    def arcsinh(self):
-1085        return self._apply_func_to_corr(np.arcsinh)
+            
1085    def arcsinh(self):
+1086        return self._apply_func_to_corr(np.arcsinh)
 

 
         

@@ -4451,8 +4454,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1087    def arccosh(self):
-1088        return self._apply_func_to_corr(np.arccosh)
+            
1088    def arccosh(self):
+1089        return self._apply_func_to_corr(np.arccosh)
 

 
         

@@ -4470,8 +4473,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1090    def arctanh(self):
-1091        return self._apply_func_to_corr(np.arctanh)
+            
1091    def arctanh(self):
+1092        return self._apply_func_to_corr(np.arctanh)
 

 
         

@@ -4509,64 +4512,64 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1126    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1127        r''' Project large correlation matrix to lowest states
-1128
-1129        This method can be used to reduce the size of an (N x N) correlation matrix
-1130        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1131        is still small.
-1132
-1133        Parameters
-1134        ----------
-1135        Ntrunc: int
-1136            Rank of the target matrix.
-1137        tproj: int
-1138            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1139            The default value is 3.
-1140        t0proj: int
-1141            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1142            discouraged for O(a) improved theories, since the correctness of the procedure
-1143            cannot be granted in this case. The default value is 2.
-1144        basematrix : Corr
-1145            Correlation matrix that is used to determine the eigenvectors of the
-1146            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1147            is is not specified.
-1148
-1149        Notes
-1150        -----
-1151        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1152        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}$
-1153        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1154        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1155        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1156        correlation matrix and to remove some noise that is added by irrelevant operators.
-1157        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1158        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1159        '''
-1160
-1161        if self.N == 1:
-1162            raise Exception('Method cannot be applied to one-dimensional correlators.')
-1163        if basematrix is None:
-1164            basematrix = self
-1165        if Ntrunc >= basematrix.N:
-1166            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
-1167        if basematrix.N != self.N:
-1168            raise Exception('basematrix and targetmatrix have to be of the same size.')
-1169
-1170        evecs = []
-1171        for i in range(Ntrunc):
-1172            evecs.append(basematrix.GEVP(t0proj, tproj, state=i))
-1173
-1174        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1175        rmat = []
-1176        for t in range(basematrix.T):
-1177            for i in range(Ntrunc):
-1178                for j in range(Ntrunc):
-1179                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1180            rmat.append(np.copy(tmpmat))
-1181
-1182        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
-1183        return Corr(newcontent)
+            
1127    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1128        r''' Project large correlation matrix to lowest states
+1129
+1130        This method can be used to reduce the size of an (N x N) correlation matrix
+1131        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1132        is still small.
+1133
+1134        Parameters
+1135        ----------
+1136        Ntrunc: int
+1137            Rank of the target matrix.
+1138        tproj: int
+1139            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1140            The default value is 3.
+1141        t0proj: int
+1142            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1143            discouraged for O(a) improved theories, since the correctness of the procedure
+1144            cannot be granted in this case. The default value is 2.
+1145        basematrix : Corr
+1146            Correlation matrix that is used to determine the eigenvectors of the
+1147            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1148            is is not specified.
+1149
+1150        Notes
+1151        -----
+1152        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1153        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}$
+1154        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1155        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1156        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1157        correlation matrix and to remove some noise that is added by irrelevant operators.
+1158        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1159        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1160        '''
+1161
+1162        if self.N == 1:
+1163            raise Exception('Method cannot be applied to one-dimensional correlators.')
+1164        if basematrix is None:
+1165            basematrix = self
+1166        if Ntrunc >= basematrix.N:
+1167            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
+1168        if basematrix.N != self.N:
+1169            raise Exception('basematrix and targetmatrix have to be of the same size.')
+1170
+1171        evecs = []
+1172        for i in range(Ntrunc):
+1173            evecs.append(basematrix.GEVP(t0proj, tproj, state=i))
+1174
+1175        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1176        rmat = []
+1177        for t in range(basematrix.T):
+1178            for i in range(Ntrunc):
+1179                for j in range(Ntrunc):
+1180                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1181            rmat.append(np.copy(tmpmat))
+1182
+1183        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
+1184        return Corr(newcontent)