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