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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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