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

202 if self.T % 2 != 0: 203 raise Exception("Can not symmetrize odd T") 204 -205 if np.argmax(np.abs(self.content)) != 0: -206 warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning) -207 -208 newcontent = [self.content[0]] -209 for t in range(1, self.T): -210 if (self.content[t] is None) or (self.content[self.T - t] is None): -211 newcontent.append(None) -212 else: -213 newcontent.append(0.5 * (self.content[t] + self.content[self.T - t])) -214 if (all([x is None for x in newcontent])): -215 raise Exception("Corr could not be symmetrized: No redundant values") -216 return Corr(newcontent, prange=self.prange) +205 if self.content[0] is not None: +206 if np.argmax(np.abs([o[0].value if o is not None else 0 for o in self.content])) != 0: +207 warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning) +208 +209 newcontent = [self.content[0]] +210 for t in range(1, self.T): +211 if (self.content[t] is None) or (self.content[self.T - t] is None): +212 newcontent.append(None) +213 else: +214 newcontent.append(0.5 * (self.content[t] + self.content[self.T - t])) +215 if (all([x is None for x in newcontent])): +216 raise Exception("Corr could not be symmetrized: No redundant values") +217 return Corr(newcontent, prange=self.prange) @@ -3170,27 +3173,27 @@ timeslice and the error on each timeslice.

- 
218    def anti_symmetric(self):
-219        """Anti-symmetrize the correlator around x0=0."""
-220        if self.N != 1:
-221            raise Exception('anti_symmetric cannot be safely applied to multi-dimensional correlators.')
-222        if self.T % 2 != 0:
-223            raise Exception("Can not symmetrize odd T")
-224
-225        test = 1 * self
-226        test.gamma_method()
-227        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
-228            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning)
-229
-230        newcontent = [self.content[0]]
-231        for t in range(1, self.T):
-232            if (self.content[t] is None) or (self.content[self.T - t] is None):
-233                newcontent.append(None)
-234            else:
-235                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
-236        if (all([x is None for x in newcontent])):
-237            raise Exception("Corr could not be symmetrized: No redundant values")
-238        return Corr(newcontent, prange=self.prange)
+            
219    def anti_symmetric(self):
+220        """Anti-symmetrize the correlator around x0=0."""
+221        if self.N != 1:
+222            raise Exception('anti_symmetric cannot be safely applied to multi-dimensional correlators.')
+223        if self.T % 2 != 0:
+224            raise Exception("Can not symmetrize odd T")
+225
+226        test = 1 * self
+227        test.gamma_method()
+228        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
+229            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning)
+230
+231        newcontent = [self.content[0]]
+232        for t in range(1, self.T):
+233            if (self.content[t] is None) or (self.content[self.T - t] is None):
+234                newcontent.append(None)
+235            else:
+236                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
+237        if (all([x is None for x in newcontent])):
+238            raise Exception("Corr could not be symmetrized: No redundant values")
+239        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -3210,20 +3213,20 @@ timeslice and the error on each timeslice.
- 
240    def is_matrix_symmetric(self):
-241        """Checks whether a correlator matrices is symmetric on every timeslice."""
-242        if self.N == 1:
-243            raise Exception("Only works for correlator matrices.")
-244        for t in range(self.T):
-245            if self[t] is None:
-246                continue
-247            for i in range(self.N):
-248                for j in range(i + 1, self.N):
-249                    if self[t][i, j] is self[t][j, i]:
-250                        continue
-251                    if hash(self[t][i, j]) != hash(self[t][j, i]):
-252                        return False
-253        return True
+            
241    def is_matrix_symmetric(self):
+242        """Checks whether a correlator matrices is symmetric on every timeslice."""
+243        if self.N == 1:
+244            raise Exception("Only works for correlator matrices.")
+245        for t in range(self.T):
+246            if self[t] is None:
+247                continue
+248            for i in range(self.N):
+249                for j in range(i + 1, self.N):
+250                    if self[t][i, j] is self[t][j, i]:
+251                        continue
+252                    if hash(self[t][i, j]) != hash(self[t][j, i]):
+253                        return False
+254        return True
 

 
 
@@ -3243,15 +3246,15 @@ timeslice and the error on each timeslice.
- 
255    def matrix_symmetric(self):
-256        """Symmetrizes the correlator matrices on every timeslice."""
-257        if self.N == 1:
-258            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
-259        if self.is_matrix_symmetric():
-260            return 1.0 * self
-261        else:
-262            transposed = [None if _check_for_none(self, G) else G.T for G in self.content]
-263            return 0.5 * (Corr(transposed) + self)
+            
256    def matrix_symmetric(self):
+257        """Symmetrizes the correlator matrices on every timeslice."""
+258        if self.N == 1:
+259            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
+260        if self.is_matrix_symmetric():
+261            return 1.0 * self
+262        else:
+263            transposed = [None if _check_for_none(self, G) else G.T for G in self.content]
+264            return 0.5 * (Corr(transposed) + self)
 

 
 
@@ -3271,84 +3274,84 @@ timeslice and the error on each timeslice.
- 
265    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
-266        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
-267
-268        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
-269        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
-270        ```python
-271        C.GEVP(t0=2)[0]  # Ground state vector(s)
-272        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
-273        ```
-274
-275        Parameters
-276        ----------
-277        t0 : int
-278            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
-279        ts : int
-280            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
-281            If sort="Eigenvector" it gives a reference point for the sorting method.
-282        sort : string
-283            If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned.
-284            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
-285            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
-286              The reference state is identified by its eigenvalue at $t=t_s$.
-287
-288        Other Parameters
-289        ----------------
-290        state : int
-291           Returns only the vector(s) for a specified state. The lowest state is zero.
-292        '''
-293
-294        if self.N == 1:
-295            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
-296        if ts is not None:
-297            if (ts <= t0):
-298                raise Exception("ts has to be larger than t0.")
-299
-300        if "sorted_list" in kwargs:
-301            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
-302            sort = kwargs.get("sorted_list")
-303
-304        if self.is_matrix_symmetric():
-305            symmetric_corr = self
-306        else:
-307            symmetric_corr = self.matrix_symmetric()
-308
-309        G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0])
-310        np.linalg.cholesky(G0)  # Check if matrix G0 is positive-semidefinite.
-311
-312        if sort is None:
-313            if (ts is None):
-314                raise Exception("ts is required if sort=None.")
-315            if (self.content[t0] is None) or (self.content[ts] is None):
-316                raise Exception("Corr not defined at t0/ts.")
-317            Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts])
-318            reordered_vecs = _GEVP_solver(Gt, G0)
-319
-320        elif sort in ["Eigenvalue", "Eigenvector"]:
-321            if sort == "Eigenvalue" and ts is not None:
-322                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
-323            all_vecs = [None] * (t0 + 1)
-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        if "state" in kwargs:
-340            return reordered_vecs[kwargs.get("state")]
-341        else:
-342            return reordered_vecs
+            
266    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
+267        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
+268
+269        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
+270        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
+271        ```python
+272        C.GEVP(t0=2)[0]  # Ground state vector(s)
+273        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
+274        ```
+275
+276        Parameters
+277        ----------
+278        t0 : int
+279            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
+280        ts : int
+281            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
+282            If sort="Eigenvector" it gives a reference point for the sorting method.
+283        sort : string
+284            If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned.
+285            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
+286            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
+287              The reference state is identified by its eigenvalue at $t=t_s$.
+288
+289        Other Parameters
+290        ----------------
+291        state : int
+292           Returns only the vector(s) for a specified state. The lowest state is zero.
+293        '''
+294
+295        if self.N == 1:
+296            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
+297        if ts is not None:
+298            if (ts <= t0):
+299                raise Exception("ts has to be larger than t0.")
+300
+301        if "sorted_list" in kwargs:
+302            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
+303            sort = kwargs.get("sorted_list")
+304
+305        if self.is_matrix_symmetric():
+306            symmetric_corr = self
+307        else:
+308            symmetric_corr = self.matrix_symmetric()
+309
+310        G0 = np.vectorize(lambda x: x.value)(symmetric_corr[t0])
+311        np.linalg.cholesky(G0)  # Check if matrix G0 is positive-semidefinite.
+312
+313        if sort is None:
+314            if (ts is None):
+315                raise Exception("ts is required if sort=None.")
+316            if (self.content[t0] is None) or (self.content[ts] is None):
+317                raise Exception("Corr not defined at t0/ts.")
+318            Gt = np.vectorize(lambda x: x.value)(symmetric_corr[ts])
+319            reordered_vecs = _GEVP_solver(Gt, G0)
+320
+321        elif sort in ["Eigenvalue", "Eigenvector"]:
+322            if sort == "Eigenvalue" and ts is not None:
+323                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
+324            all_vecs = [None] * (t0 + 1)
+325            for t in range(t0 + 1, self.T):
+326                try:
+327                    Gt = np.vectorize(lambda x: x.value)(symmetric_corr[t])
+328                    all_vecs.append(_GEVP_solver(Gt, G0))
+329                except Exception:
+330                    all_vecs.append(None)
+331            if sort == "Eigenvector":
+332                if ts is None:
+333                    raise Exception("ts is required for the Eigenvector sorting method.")
+334                all_vecs = _sort_vectors(all_vecs, ts)
+335
+336            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
+337        else:
+338            raise Exception("Unkown value for 'sort'.")
+339
+340        if "state" in kwargs:
+341            return reordered_vecs[kwargs.get("state")]
+342        else:
+343            return reordered_vecs
 

 
 
@@ -3401,18 +3404,18 @@ Returns only the vector(s) for a specified state. The lowest state is zero.
- 
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)
+            
345    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
+346        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
+347
+348        Parameters
+349        ----------
+350        state : int
+351            The state one is interested in ordered by energy. The lowest state is zero.
+352
+353        All other parameters are identical to the ones of Corr.GEVP.
+354        """
+355        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
+356        return self.projected(vec)
 

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

 
 
@@ -3513,15 +3516,15 @@ determines whether the matrix is extended periodically
- 
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)))
+            
399    def roll(self, dt):
+400        """Periodically shift the correlator by dt timeslices
+401
+402        Parameters
+403        ----------
+404        dt : int
+405            number of timeslices
+406        """
+407        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
 

 
 
@@ -3548,9 +3551,9 @@ number of timeslices
- 
408    def reverse(self):
-409        """Reverse the time ordering of the Corr"""
-410        return Corr(self.content[:: -1])
+            
409    def reverse(self):
+410        """Reverse the time ordering of the Corr"""
+411        return Corr(self.content[:: -1])
 

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

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

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

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

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

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

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

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

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

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

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

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

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

 
 
@@ -4497,8 +4500,8 @@ specifies a custom path for the file (default '.')
- 
961    def print(self, print_range=None):
-962        print(self.__repr__(print_range))
+            
962    def print(self, print_range=None):
+963        print(self.__repr__(print_range))
 

 
 
@@ -4516,8 +4519,8 @@ specifies a custom path for the file (default '.')
- 
1128    def sqrt(self):
-1129        return self ** 0.5
+            
1129    def sqrt(self):
+1130        return self ** 0.5
 

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

 
 
@@ -4555,9 +4558,9 @@ specifies a custom path for the file (default '.')
- 
1135    def exp(self):
-1136        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
-1137        return Corr(newcontent, prange=self.prange)
+            
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)
 

 
 
@@ -4575,8 +4578,8 @@ specifies a custom path for the file (default '.')
- 
1152    def sin(self):
-1153        return self._apply_func_to_corr(np.sin)
+            
1153    def sin(self):
+1154        return self._apply_func_to_corr(np.sin)
 

 
 
@@ -4594,8 +4597,8 @@ specifies a custom path for the file (default '.')
- 
1155    def cos(self):
-1156        return self._apply_func_to_corr(np.cos)
+            
1156    def cos(self):
+1157        return self._apply_func_to_corr(np.cos)
 

 
 
@@ -4613,8 +4616,8 @@ specifies a custom path for the file (default '.')
- 
1158    def tan(self):
-1159        return self._apply_func_to_corr(np.tan)
+            
1159    def tan(self):
+1160        return self._apply_func_to_corr(np.tan)
 

 
 
@@ -4632,8 +4635,8 @@ specifies a custom path for the file (default '.')
- 
1161    def sinh(self):
-1162        return self._apply_func_to_corr(np.sinh)
+            
1162    def sinh(self):
+1163        return self._apply_func_to_corr(np.sinh)
 

 
 
@@ -4651,8 +4654,8 @@ specifies a custom path for the file (default '.')
- 
1164    def cosh(self):
-1165        return self._apply_func_to_corr(np.cosh)
+            
1165    def cosh(self):
+1166        return self._apply_func_to_corr(np.cosh)
 

 
 
@@ -4670,8 +4673,8 @@ specifies a custom path for the file (default '.')
- 
1167    def tanh(self):
-1168        return self._apply_func_to_corr(np.tanh)
+            
1168    def tanh(self):
+1169        return self._apply_func_to_corr(np.tanh)
 

 
 
@@ -4689,8 +4692,8 @@ specifies a custom path for the file (default '.')
- 
1170    def arcsin(self):
-1171        return self._apply_func_to_corr(np.arcsin)
+            
1171    def arcsin(self):
+1172        return self._apply_func_to_corr(np.arcsin)
 

 
 
@@ -4708,8 +4711,8 @@ specifies a custom path for the file (default '.')
- 
1173    def arccos(self):
-1174        return self._apply_func_to_corr(np.arccos)
+            
1174    def arccos(self):
+1175        return self._apply_func_to_corr(np.arccos)
 

 
 
@@ -4727,8 +4730,8 @@ specifies a custom path for the file (default '.')
- 
1176    def arctan(self):
-1177        return self._apply_func_to_corr(np.arctan)
+            
1177    def arctan(self):
+1178        return self._apply_func_to_corr(np.arctan)
 

 
 
@@ -4746,8 +4749,8 @@ specifies a custom path for the file (default '.')
- 
1179    def arcsinh(self):
-1180        return self._apply_func_to_corr(np.arcsinh)
+            
1180    def arcsinh(self):
+1181        return self._apply_func_to_corr(np.arcsinh)
 

 
 
@@ -4765,8 +4768,8 @@ specifies a custom path for the file (default '.')
- 
1182    def arccosh(self):
-1183        return self._apply_func_to_corr(np.arccosh)
+            
1183    def arccosh(self):
+1184        return self._apply_func_to_corr(np.arccosh)
 

 
 
@@ -4784,8 +4787,8 @@ specifies a custom path for the file (default '.')
- 
1185    def arctanh(self):
-1186        return self._apply_func_to_corr(np.arctanh)
+            
1186    def arctanh(self):
+1187        return self._apply_func_to_corr(np.arctanh)
 

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