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