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)
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
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)
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
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)
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)
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)))
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])
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)
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)
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)
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
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.")
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.")
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.')
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
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)
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
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.")
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()
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))
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))
1128 def sqrt(self):
-1129 return self ** 0.5
+ 1129 def sqrt(self):
+1130 return self ** 0.5
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)