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

259 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. 260 The reference state is identified by its eigenvalue at t=ts 261 """ -262 symmetric_corr = self.matrix_symmetric() -263 if sorted_list is None: -264 if (ts is None): -265 raise Exception("ts is required if sorted_list=None") -266 if (self.content[t0] is None) or (self.content[ts] is None): -267 raise Exception("Corr not defined at t0/ts") -268 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") -269 for i in range(self.N): -270 for j in range(self.N): -271 G0[i, j] = symmetric_corr[t0][i, j].value -272 Gt[i, j] = symmetric_corr[ts][i, j].value -273 -274 sp_vecs = _GEVP_solver(Gt, G0) -275 sp_vec = sp_vecs[state] -276 return sp_vec -277 elif sorted_list in ["Eigenvalue", "Eigenvector"]: -278 all_vecs = [] -279 for t in range(self.T): -280 try: -281 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") -282 for i in range(self.N): -283 for j in range(self.N): -284 G0[i, j] = symmetric_corr[t0][i, j].value -285 Gt[i, j] = symmetric_corr[t][i, j].value -286 -287 sp_vecs = _GEVP_solver(Gt, G0) -288 if sorted_list == "Eigenvalue": -289 sp_vec = sp_vecs[state] -290 all_vecs.append(sp_vec) -291 else: -292 all_vecs.append(sp_vecs) -293 except Exception: -294 all_vecs.append(None) -295 if sorted_list == "Eigenvector": -296 if (ts is None): -297 raise Exception("ts is required for the Eigenvector sorting method.") -298 all_vecs = _sort_vectors(all_vecs, ts) -299 all_vecs = [a[state] for a in all_vecs] -300 else: -301 raise Exception("Unkown value for 'sorted_list'.") -302 -303 return all_vecs +262 +263 if self.N == 1: +264 raise Exception("GEVP methods only works on correlator matrices and not single correlators.") +265 +266 symmetric_corr = self.matrix_symmetric() +267 if sorted_list is None: +268 if (ts is None): +269 raise Exception("ts is required if sorted_list=None") +270 if (self.content[t0] is None) or (self.content[ts] is None): +271 raise Exception("Corr not defined at t0/ts") +272 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") +273 for i in range(self.N): +274 for j in range(self.N): +275 G0[i, j] = symmetric_corr[t0][i, j].value +276 Gt[i, j] = symmetric_corr[ts][i, j].value +277 +278 sp_vecs = _GEVP_solver(Gt, G0) +279 sp_vec = sp_vecs[state] +280 return sp_vec +281 elif sorted_list in ["Eigenvalue", "Eigenvector"]: +282 all_vecs = [] +283 for t in range(self.T): +284 try: +285 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") +286 for i in range(self.N): +287 for j in range(self.N): +288 G0[i, j] = symmetric_corr[t0][i, j].value +289 Gt[i, j] = symmetric_corr[t][i, j].value +290 +291 sp_vecs = _GEVP_solver(Gt, G0) +292 if sorted_list == "Eigenvalue": +293 sp_vec = sp_vecs[state] +294 all_vecs.append(sp_vec) +295 else: +296 all_vecs.append(sp_vecs) +297 except Exception: +298 all_vecs.append(None) +299 if sorted_list == "Eigenvector": +300 if (ts is None): +301 raise Exception("ts is required for the Eigenvector sorting method.") +302 all_vecs = _sort_vectors(all_vecs, ts) +303 all_vecs = [a[state] for a in all_vecs] +304 else: +305 raise Exception("Unkown value for 'sorted_list'.") +306 +307 return all_vecs @@ -3120,26 +3132,26 @@ if this argument is set, a list of vectors (len=self.T) is returned. If it is le
View Source - 
305    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
-306        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
-307
-308        Parameters
-309        ----------
-310        t0 : int
-311            The time t0 for G(t)v= lambda G(t_0)v
-312        ts : int
-313            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
-314            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
-315        state : int
-316            The state one is interested in ordered by energy. The lowest state is zero.
-317        sorted_list : string
-318            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
-319             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
-320             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
-321                              The reference state is identified by its eigenvalue at t=ts
-322        """
-323        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
-324        return self.projected(vec)
+            
309    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
+310        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
+311
+312        Parameters
+313        ----------
+314        t0 : int
+315            The time t0 for G(t)v= lambda G(t_0)v
+316        ts : int
+317            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
+318            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
+319        state : int
+320            The state one is interested in ordered by energy. The lowest state is zero.
+321        sorted_list : string
+322            if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
+323             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
+324             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
+325                              The reference state is identified by its eigenvalue at t=ts
+326        """
+327        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
+328        return self.projected(vec)
@@ -3176,46 +3188,46 @@ if this argument is set, a list of vectors (len=self.T) is returned. If it is le
View Source - 
326    def Hankel(self, N, periodic=False):
-327        """Constructs an NxN Hankel matrix
-328
-329        C(t) c(t+1) ... c(t+n-1)
-330        C(t+1) c(t+2) ... c(t+n)
-331        .................
-332        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
-333
-334        Parameters
-335        ----------
-336        N : int
-337            Dimension of the Hankel matrix
-338        periodic : bool, optional
-339            determines whether the matrix is extended periodically
-340        """
-341
-342        if self.N != 1:
-343            raise Exception("Multi-operator Prony not implemented!")
-344
-345        array = np.empty([N, N], dtype="object")
-346        new_content = []
-347        for t in range(self.T):
-348            new_content.append(array.copy())
-349
-350        def wrap(i):
-351            while i >= self.T:
-352                i -= self.T
-353            return i
-354
-355        for t in range(self.T):
-356            for i in range(N):
-357                for j in range(N):
-358                    if periodic:
-359                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
-360                    elif (t + i + j) >= self.T:
-361                        new_content[t] = None
-362                    else:
-363                        new_content[t][i, j] = self.content[t + i + j][0]
-364
-365        return Corr(new_content)
+            
330    def Hankel(self, N, periodic=False):
+331        """Constructs an NxN Hankel matrix
+332
+333        C(t) c(t+1) ... c(t+n-1)
+334        C(t+1) c(t+2) ... c(t+n)
+335        .................
+336        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
+337
+338        Parameters
+339        ----------
+340        N : int
+341            Dimension of the Hankel matrix
+342        periodic : bool, optional
+343            determines whether the matrix is extended periodically
+344        """
+345
+346        if self.N != 1:
+347            raise Exception("Multi-operator Prony not implemented!")
+348
+349        array = np.empty([N, N], dtype="object")
+350        new_content = []
+351        for t in range(self.T):
+352            new_content.append(array.copy())
+353
+354        def wrap(i):
+355            while i >= self.T:
+356                i -= self.T
+357            return i
+358
+359        for t in range(self.T):
+360            for i in range(N):
+361                for j in range(N):
+362                    if periodic:
+363                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
+364                    elif (t + i + j) >= self.T:
+365                        new_content[t] = None
+366                    else:
+367                        new_content[t][i, j] = self.content[t + i + j][0]
+368
+369        return Corr(new_content)
@@ -3249,15 +3261,15 @@ determines whether the matrix is extended periodically
View Source - 
367    def roll(self, dt):
-368        """Periodically shift the correlator by dt timeslices
-369
-370        Parameters
-371        ----------
-372        dt : int
-373            number of timeslices
-374        """
-375        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
+            
371    def roll(self, dt):
+372        """Periodically shift the correlator by dt timeslices
+373
+374        Parameters
+375        ----------
+376        dt : int
+377            number of timeslices
+378        """
+379        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
@@ -3284,9 +3296,9 @@ number of timeslices
View Source - 
377    def reverse(self):
-378        """Reverse the time ordering of the Corr"""
-379        return Corr(self.content[:: -1])
+            
381    def reverse(self):
+382        """Reverse the time ordering of the Corr"""
+383        return Corr(self.content[:: -1])
@@ -3306,23 +3318,23 @@ number of timeslices
View Source - 
381    def thin(self, spacing=2, offset=0):
-382        """Thin out a correlator to suppress correlations
-383
-384        Parameters
-385        ----------
-386        spacing : int
-387            Keep only every 'spacing'th entry of the correlator
-388        offset : int
-389            Offset the equal spacing
-390        """
-391        new_content = []
-392        for t in range(self.T):
-393            if (offset + t) % spacing != 0:
-394                new_content.append(None)
-395            else:
-396                new_content.append(self.content[t])
-397        return Corr(new_content)
+            
385    def thin(self, spacing=2, offset=0):
+386        """Thin out a correlator to suppress correlations
+387
+388        Parameters
+389        ----------
+390        spacing : int
+391            Keep only every 'spacing'th entry of the correlator
+392        offset : int
+393            Offset the equal spacing
+394        """
+395        new_content = []
+396        for t in range(self.T):
+397            if (offset + t) % spacing != 0:
+398                new_content.append(None)
+399            else:
+400                new_content.append(self.content[t])
+401        return Corr(new_content)
@@ -3351,32 +3363,32 @@ Offset the equal spacing
View Source - 
399    def correlate(self, partner):
-400        """Correlate the correlator with another correlator or Obs
-401
-402        Parameters
-403        ----------
-404        partner : Obs or Corr
-405            partner to correlate the correlator with.
-406            Can either be an Obs which is correlated with all entries of the
-407            correlator or a Corr of same length.
-408        """
-409        new_content = []
-410        for x0, t_slice in enumerate(self.content):
-411            if t_slice is None:
-412                new_content.append(None)
-413            else:
-414                if isinstance(partner, Corr):
-415                    if partner.content[x0] is None:
-416                        new_content.append(None)
-417                    else:
-418                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
-419                elif isinstance(partner, Obs):  # Should this include CObs?
-420                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
-421                else:
-422                    raise Exception("Can only correlate with an Obs or a Corr.")
-423
-424        return Corr(new_content)
+            
403    def correlate(self, partner):
+404        """Correlate the correlator with another correlator or Obs
+405
+406        Parameters
+407        ----------
+408        partner : Obs or Corr
+409            partner to correlate the correlator with.
+410            Can either be an Obs which is correlated with all entries of the
+411            correlator or a Corr of same length.
+412        """
+413        new_content = []
+414        for x0, t_slice in enumerate(self.content):
+415            if t_slice is None:
+416                new_content.append(None)
+417            else:
+418                if isinstance(partner, Corr):
+419                    if partner.content[x0] is None:
+420                        new_content.append(None)
+421                    else:
+422                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
+423                elif isinstance(partner, Obs):  # Should this include CObs?
+424                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
+425                else:
+426                    raise Exception("Can only correlate with an Obs or a Corr.")
+427
+428        return Corr(new_content)
@@ -3405,26 +3417,26 @@ correlator or a Corr of same length.
View Source - 
426    def reweight(self, weight, **kwargs):
-427        """Reweight the correlator.
-428
-429        Parameters
-430        ----------
-431        weight : Obs
-432            Reweighting factor. An Observable that has to be defined on a superset of the
-433            configurations in obs[i].idl for all i.
-434        all_configs : bool
-435            if True, the reweighted observables are normalized by the average of
-436            the reweighting factor on all configurations in weight.idl and not
-437            on the configurations in obs[i].idl.
-438        """
-439        new_content = []
-440        for t_slice in self.content:
-441            if t_slice is None:
-442                new_content.append(None)
-443            else:
-444                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
-445        return Corr(new_content)
+            
430    def reweight(self, weight, **kwargs):
+431        """Reweight the correlator.
+432
+433        Parameters
+434        ----------
+435        weight : Obs
+436            Reweighting factor. An Observable that has to be defined on a superset of the
+437            configurations in obs[i].idl for all i.
+438        all_configs : bool
+439            if True, the reweighted observables are normalized by the average of
+440            the reweighting factor on all configurations in weight.idl and not
+441            on the configurations in obs[i].idl.
+442        """
+443        new_content = []
+444        for t_slice in self.content:
+445            if t_slice is None:
+446                new_content.append(None)
+447            else:
+448                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
+449        return Corr(new_content)
@@ -3456,33 +3468,33 @@ on the configurations in obs[i].idl.
View Source - 
447    def T_symmetry(self, partner, parity=+1):
-448        """Return the time symmetry average of the correlator and its partner
-449
-450        Parameters
-451        ----------
-452        partner : Corr
-453            Time symmetry partner of the Corr
-454        partity : int
-455            Parity quantum number of the correlator, can be +1 or -1
-456        """
-457        if not isinstance(partner, Corr):
-458            raise Exception("T partner has to be a Corr object.")
-459        if parity not in [+1, -1]:
-460            raise Exception("Parity has to be +1 or -1.")
-461        T_partner = parity * partner.reverse()
-462
-463        t_slices = []
-464        test = (self - T_partner)
-465        test.gamma_method()
-466        for x0, t_slice in enumerate(test.content):
-467            if t_slice is not None:
-468                if not t_slice[0].is_zero_within_error(5):
-469                    t_slices.append(x0)
-470        if t_slices:
-471            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
-472
-473        return (self + T_partner) / 2
+            
451    def T_symmetry(self, partner, parity=+1):
+452        """Return the time symmetry average of the correlator and its partner
+453
+454        Parameters
+455        ----------
+456        partner : Corr
+457            Time symmetry partner of the Corr
+458        partity : int
+459            Parity quantum number of the correlator, can be +1 or -1
+460        """
+461        if not isinstance(partner, Corr):
+462            raise Exception("T partner has to be a Corr object.")
+463        if parity not in [+1, -1]:
+464            raise Exception("Parity has to be +1 or -1.")
+465        T_partner = parity * partner.reverse()
+466
+467        t_slices = []
+468        test = (self - T_partner)
+469        test.gamma_method()
+470        for x0, t_slice in enumerate(test.content):
+471            if t_slice is not None:
+472                if not t_slice[0].is_zero_within_error(5):
+473                    t_slices.append(x0)
+474        if t_slices:
+475            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
+476
+477        return (self + T_partner) / 2
@@ -3511,57 +3523,57 @@ Parity quantum number of the correlator, can be +1 or -1
View Source - 
475    def deriv(self, variant="symmetric"):
-476        """Return the first derivative of the correlator with respect to x0.
-477
-478        Parameters
-479        ----------
-480        variant : str
-481            decides which definition of the finite differences derivative is used.
-482            Available choice: symmetric, forward, backward, improved, default: symmetric
-483        """
-484        if variant == "symmetric":
-485            newcontent = []
-486            for t in range(1, self.T - 1):
-487                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-488                    newcontent.append(None)
-489                else:
-490                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
-491            if(all([x is None for x in newcontent])):
-492                raise Exception('Derivative is undefined at all timeslices')
-493            return Corr(newcontent, padding=[1, 1])
-494        elif variant == "forward":
-495            newcontent = []
-496            for t in range(self.T - 1):
-497                if (self.content[t] is None) or (self.content[t + 1] is None):
-498                    newcontent.append(None)
-499                else:
-500                    newcontent.append(self.content[t + 1] - self.content[t])
-501            if(all([x is None for x in newcontent])):
-502                raise Exception("Derivative is undefined at all timeslices")
-503            return Corr(newcontent, padding=[0, 1])
-504        elif variant == "backward":
-505            newcontent = []
-506            for t in range(1, self.T):
-507                if (self.content[t - 1] is None) or (self.content[t] is None):
-508                    newcontent.append(None)
-509                else:
-510                    newcontent.append(self.content[t] - self.content[t - 1])
-511            if(all([x is None for x in newcontent])):
-512                raise Exception("Derivative is undefined at all timeslices")
-513            return Corr(newcontent, padding=[1, 0])
-514        elif variant == "improved":
-515            newcontent = []
-516            for t in range(2, self.T - 2):
-517                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):
-518                    newcontent.append(None)
-519                else:
-520                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
-521            if(all([x is None for x in newcontent])):
-522                raise Exception('Derivative is undefined at all timeslices')
-523            return Corr(newcontent, padding=[2, 2])
-524        else:
-525            raise Exception("Unknown variant.")
+            
479    def deriv(self, variant="symmetric"):
+480        """Return the first derivative of the correlator with respect to x0.
+481
+482        Parameters
+483        ----------
+484        variant : str
+485            decides which definition of the finite differences derivative is used.
+486            Available choice: symmetric, forward, backward, improved, default: symmetric
+487        """
+488        if variant == "symmetric":
+489            newcontent = []
+490            for t in range(1, self.T - 1):
+491                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+492                    newcontent.append(None)
+493                else:
+494                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
+495            if(all([x is None for x in newcontent])):
+496                raise Exception('Derivative is undefined at all timeslices')
+497            return Corr(newcontent, padding=[1, 1])
+498        elif variant == "forward":
+499            newcontent = []
+500            for t in range(self.T - 1):
+501                if (self.content[t] is None) or (self.content[t + 1] is None):
+502                    newcontent.append(None)
+503                else:
+504                    newcontent.append(self.content[t + 1] - self.content[t])
+505            if(all([x is None for x in newcontent])):
+506                raise Exception("Derivative is undefined at all timeslices")
+507            return Corr(newcontent, padding=[0, 1])
+508        elif variant == "backward":
+509            newcontent = []
+510            for t in range(1, self.T):
+511                if (self.content[t - 1] is None) or (self.content[t] is None):
+512                    newcontent.append(None)
+513                else:
+514                    newcontent.append(self.content[t] - self.content[t - 1])
+515            if(all([x is None for x in newcontent])):
+516                raise Exception("Derivative is undefined at all timeslices")
+517            return Corr(newcontent, padding=[1, 0])
+518        elif variant == "improved":
+519            newcontent = []
+520            for t in range(2, self.T - 2):
+521                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):
+522                    newcontent.append(None)
+523                else:
+524                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
+525            if(all([x is None for x in newcontent])):
+526                raise Exception('Derivative is undefined at all timeslices')
+527            return Corr(newcontent, padding=[2, 2])
+528        else:
+529            raise Exception("Unknown variant.")
@@ -3589,37 +3601,37 @@ Available choice: symmetric, forward, backward, improved, default: symmetric View Source - 
527    def second_deriv(self, variant="symmetric"):
-528        """Return the second derivative of the correlator with respect to x0.
-529
-530        Parameters
-531        ----------
-532        variant : str
-533            decides which definition of the finite differences derivative is used.
-534            Available choice: symmetric, improved, default: symmetric
-535        """
-536        if variant == "symmetric":
-537            newcontent = []
-538            for t in range(1, self.T - 1):
-539                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-540                    newcontent.append(None)
-541                else:
-542                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
-543            if(all([x is None for x in newcontent])):
-544                raise Exception("Derivative is undefined at all timeslices")
-545            return Corr(newcontent, padding=[1, 1])
-546        elif variant == "improved":
-547            newcontent = []
-548            for t in range(2, self.T - 2):
-549                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):
-550                    newcontent.append(None)
-551                else:
-552                    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]))
-553            if(all([x is None for x in newcontent])):
-554                raise Exception("Derivative is undefined at all timeslices")
-555            return Corr(newcontent, padding=[2, 2])
-556        else:
-557            raise Exception("Unknown variant.")
+            
531    def second_deriv(self, variant="symmetric"):
+532        """Return the second derivative of the correlator with respect to x0.
+533
+534        Parameters
+535        ----------
+536        variant : str
+537            decides which definition of the finite differences derivative is used.
+538            Available choice: symmetric, improved, default: symmetric
+539        """
+540        if variant == "symmetric":
+541            newcontent = []
+542            for t in range(1, self.T - 1):
+543                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+544                    newcontent.append(None)
+545                else:
+546                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
+547            if(all([x is None for x in newcontent])):
+548                raise Exception("Derivative is undefined at all timeslices")
+549            return Corr(newcontent, padding=[1, 1])
+550        elif variant == "improved":
+551            newcontent = []
+552            for t in range(2, self.T - 2):
+553                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):
+554                    newcontent.append(None)
+555                else:
+556                    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]))
+557            if(all([x is None for x in newcontent])):
+558                raise Exception("Derivative is undefined at all timeslices")
+559            return Corr(newcontent, padding=[2, 2])
+560        else:
+561            raise Exception("Unknown variant.")
 

 
         
@@ -3647,70 +3659,70 @@ Available choice: symmetric, improved, default: symmetric
 
             

             View Source
-            
559    def m_eff(self, variant='log', guess=1.0):
-560        """Returns the effective mass of the correlator as correlator object
-561
-562        Parameters
-563        ----------
-564        variant : str
-565            log : uses the standard effective mass log(C(t) / C(t+1))
-566            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.
-567            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.
-568            See, e.g., arXiv:1205.5380
-569            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
-570        guess : float
-571            guess for the root finder, only relevant for the root variant
-572        """
-573        if self.N != 1:
-574            raise Exception('Correlator must be projected before getting m_eff')
-575        if variant == 'log':
-576            newcontent = []
-577            for t in range(self.T - 1):
-578                if (self.content[t] is None) or (self.content[t + 1] is None):
-579                    newcontent.append(None)
-580                else:
-581                    newcontent.append(self.content[t] / self.content[t + 1])
-582            if(all([x is None for x in newcontent])):
-583                raise Exception('m_eff is undefined at all timeslices')
-584
-585            return np.log(Corr(newcontent, padding=[0, 1]))
-586
-587        elif variant in ['periodic', 'cosh', 'sinh']:
-588            if variant in ['periodic', 'cosh']:
-589                func = anp.cosh
-590            else:
-591                func = anp.sinh
-592
-593            def root_function(x, d):
-594                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
-595
-596            newcontent = []
-597            for t in range(self.T - 1):
-598                if (self.content[t] is None) or (self.content[t + 1] is None):
-599                    newcontent.append(None)
-600                # Fill the two timeslices in the middle of the lattice with their predecessors
-601                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
-602                    newcontent.append(newcontent[-1])
-603                else:
-604                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
-605            if(all([x is None for x in newcontent])):
-606                raise Exception('m_eff is undefined at all timeslices')
-607
-608            return Corr(newcontent, padding=[0, 1])
-609
-610        elif variant == 'arccosh':
-611            newcontent = []
-612            for t in range(1, self.T - 1):
-613                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
-614                    newcontent.append(None)
-615                else:
-616                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
-617            if(all([x is None for x in newcontent])):
-618                raise Exception("m_eff is undefined at all timeslices")
-619            return np.arccosh(Corr(newcontent, padding=[1, 1]))
-620
-621        else:
-622            raise Exception('Unknown variant.')
+            
563    def m_eff(self, variant='log', guess=1.0):
+564        """Returns the effective mass of the correlator as correlator object
+565
+566        Parameters
+567        ----------
+568        variant : str
+569            log : uses the standard effective mass log(C(t) / C(t+1))
+570            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.
+571            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.
+572            See, e.g., arXiv:1205.5380
+573            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
+574        guess : float
+575            guess for the root finder, only relevant for the root variant
+576        """
+577        if self.N != 1:
+578            raise Exception('Correlator must be projected before getting m_eff')
+579        if variant == 'log':
+580            newcontent = []
+581            for t in range(self.T - 1):
+582                if (self.content[t] is None) or (self.content[t + 1] is None):
+583                    newcontent.append(None)
+584                else:
+585                    newcontent.append(self.content[t] / self.content[t + 1])
+586            if(all([x is None for x in newcontent])):
+587                raise Exception('m_eff is undefined at all timeslices')
+588
+589            return np.log(Corr(newcontent, padding=[0, 1]))
+590
+591        elif variant in ['periodic', 'cosh', 'sinh']:
+592            if variant in ['periodic', 'cosh']:
+593                func = anp.cosh
+594            else:
+595                func = anp.sinh
+596
+597            def root_function(x, d):
+598                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
+599
+600            newcontent = []
+601            for t in range(self.T - 1):
+602                if (self.content[t] is None) or (self.content[t + 1] is None):
+603                    newcontent.append(None)
+604                # Fill the two timeslices in the middle of the lattice with their predecessors
+605                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+606                    newcontent.append(newcontent[-1])
+607                else:
+608                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+609            if(all([x is None for x in newcontent])):
+610                raise Exception('m_eff is undefined at all timeslices')
+611
+612            return Corr(newcontent, padding=[0, 1])
+613
+614        elif variant == 'arccosh':
+615            newcontent = []
+616            for t in range(1, self.T - 1):
+617                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
+618                    newcontent.append(None)
+619                else:
+620                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+621            if(all([x is None for x in newcontent])):
+622                raise Exception("m_eff is undefined at all timeslices")
+623            return np.arccosh(Corr(newcontent, padding=[1, 1]))
+624
+625        else:
+626            raise Exception('Unknown variant.')
 

 
         

@@ -3743,39 +3755,39 @@ guess for the root finder, only relevant for the root variant
 
             

             View Source
-            
624    def fit(self, function, fitrange=None, silent=False, **kwargs):
-625        r'''Fits function to the data
-626
-627        Parameters
-628        ----------
-629        function : obj
-630            function to fit to the data. See fits.least_squares for details.
-631        fitrange : list
-632            Two element list containing the timeslices on which the fit is supposed to start and stop.
-633            Caution: This range is inclusive as opposed to standard python indexing.
-634            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
-635            If not specified, self.prange or all timeslices are used.
-636        silent : bool
-637            Decides whether output is printed to the standard output.
-638        '''
-639        if self.N != 1:
-640            raise Exception("Correlator must be projected before fitting")
-641
-642        if fitrange is None:
-643            if self.prange:
-644                fitrange = self.prange
-645            else:
-646                fitrange = [0, self.T - 1]
-647        else:
-648            if not isinstance(fitrange, list):
-649                raise Exception("fitrange has to be a list with two elements")
-650            if len(fitrange) != 2:
-651                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
-652
-653        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-654        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-655        result = least_squares(xs, ys, function, silent=silent, **kwargs)
-656        return result
+            
628    def fit(self, function, fitrange=None, silent=False, **kwargs):
+629        r'''Fits function to the data
+630
+631        Parameters
+632        ----------
+633        function : obj
+634            function to fit to the data. See fits.least_squares for details.
+635        fitrange : list
+636            Two element list containing the timeslices on which the fit is supposed to start and stop.
+637            Caution: This range is inclusive as opposed to standard python indexing.
+638            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+639            If not specified, self.prange or all timeslices are used.
+640        silent : bool
+641            Decides whether output is printed to the standard output.
+642        '''
+643        if self.N != 1:
+644            raise Exception("Correlator must be projected before fitting")
+645
+646        if fitrange is None:
+647            if self.prange:
+648                fitrange = self.prange
+649            else:
+650                fitrange = [0, self.T - 1]
+651        else:
+652            if not isinstance(fitrange, list):
+653                raise Exception("fitrange has to be a list with two elements")
+654            if len(fitrange) != 2:
+655                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
+656
+657        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+658        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+659        result = least_squares(xs, ys, function, silent=silent, **kwargs)
+660        return result
 

 
         

@@ -3809,42 +3821,42 @@ Decides whether output is printed to the standard output.
 
             

             View Source
-            
658    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
-659        """ Extract a plateau value from a Corr object
-660
-661        Parameters
-662        ----------
-663        plateau_range : list
-664            list with two entries, indicating the first and the last timeslice
-665            of the plateau region.
-666        method : str
-667            method to extract the plateau.
-668                'fit' fits a constant to the plateau region
-669                'avg', 'average' or 'mean' just average over the given timeslices.
-670        auto_gamma : bool
-671            apply gamma_method with default parameters to the Corr. Defaults to None
-672        """
-673        if not plateau_range:
-674            if self.prange:
-675                plateau_range = self.prange
-676            else:
-677                raise Exception("no plateau range provided")
-678        if self.N != 1:
-679            raise Exception("Correlator must be projected before getting a plateau.")
-680        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
-681            raise Exception("plateau is undefined at all timeslices in plateaurange.")
-682        if auto_gamma:
-683            self.gamma_method()
-684        if method == "fit":
-685            def const_func(a, t):
-686                return a[0]
-687            return self.fit(const_func, plateau_range)[0]
-688        elif method in ["avg", "average", "mean"]:
-689            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
-690            return returnvalue
-691
-692        else:
-693            raise Exception("Unsupported plateau method: " + method)
+            
662    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+663        """ Extract a plateau value from a Corr object
+664
+665        Parameters
+666        ----------
+667        plateau_range : list
+668            list with two entries, indicating the first and the last timeslice
+669            of the plateau region.
+670        method : str
+671            method to extract the plateau.
+672                'fit' fits a constant to the plateau region
+673                'avg', 'average' or 'mean' just average over the given timeslices.
+674        auto_gamma : bool
+675            apply gamma_method with default parameters to the Corr. Defaults to None
+676        """
+677        if not plateau_range:
+678            if self.prange:
+679                plateau_range = self.prange
+680            else:
+681                raise Exception("no plateau range provided")
+682        if self.N != 1:
+683            raise Exception("Correlator must be projected before getting a plateau.")
+684        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+685            raise Exception("plateau is undefined at all timeslices in plateaurange.")
+686        if auto_gamma:
+687            self.gamma_method()
+688        if method == "fit":
+689            def const_func(a, t):
+690                return a[0]
+691            return self.fit(const_func, plateau_range)[0]
+692        elif method in ["avg", "average", "mean"]:
+693            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+694            return returnvalue
+695
+696        else:
+697            raise Exception("Unsupported plateau method: " + method)
 

 
         

@@ -3878,17 +3890,17 @@ apply gamma_method with default parameters to the Corr. Defaults to None
 
             

             View Source
-            
695    def set_prange(self, prange):
-696        """Sets the attribute prange of the Corr object."""
-697        if not len(prange) == 2:
-698            raise Exception("prange must be a list or array with two values")
-699        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
-700            raise Exception("Start and end point must be integers")
-701        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
-702            raise Exception("Start and end point must define a range in the interval 0,T")
-703
-704        self.prange = prange
-705        return
+            
699    def set_prange(self, prange):
+700        """Sets the attribute prange of the Corr object."""
+701        if not len(prange) == 2:
+702            raise Exception("prange must be a list or array with two values")
+703        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+704            raise Exception("Start and end point must be integers")
+705        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
+706            raise Exception("Start and end point must define a range in the interval 0,T")
+707
+708        self.prange = prange
+709        return
 

 
         

@@ -3921,118 +3933,118 @@ apply gamma_method with default parameters to the Corr. Defaults to None
 
             

             View Source
-            
707    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):
-708        """Plots the correlator using the tag of the correlator as label if available.
-709
-710        Parameters
-711        ----------
-712        x_range : list
-713            list of two values, determining the range of the x-axis e.g. [4, 8]
-714        comp : Corr or list of Corr
-715            Correlator or list of correlators which are plotted for comparison.
-716            The tags of these correlators are used as labels if available.
-717        logscale : bool
-718            Sets y-axis to logscale
-719        plateau : Obs
-720            Plateau value to be visualized in the figure
-721        fit_res : Fit_result
-722            Fit_result object to be visualized
-723        ylabel : str
-724            Label for the y-axis
-725        save : str
-726            path to file in which the figure should be saved
-727        auto_gamma : bool
-728            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
-729        hide_sigma : float
-730            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
-731        references : list
-732            List of floating point values that are displayed as horizontal lines for reference.
-733        """
-734        if self.N != 1:
-735            raise Exception("Correlator must be projected before plotting")
-736
-737        if auto_gamma:
-738            self.gamma_method()
-739
-740        if x_range is None:
-741            x_range = [0, self.T - 1]
-742
-743        fig = plt.figure()
-744        ax1 = fig.add_subplot(111)
-745
-746        x, y, y_err = self.plottable()
-747        if hide_sigma:
-748            hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
-749        else:
-750            hide_from = None
-751        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
-752        if logscale:
-753            ax1.set_yscale('log')
-754        else:
-755            if y_range is None:
-756                try:
-757                    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)])
-758                    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)])
-759                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
-760                except Exception:
-761                    pass
-762            else:
-763                ax1.set_ylim(y_range)
-764        if comp:
-765            if isinstance(comp, (Corr, list)):
-766                for corr in comp if isinstance(comp, list) else [comp]:
-767                    if auto_gamma:
-768                        corr.gamma_method()
-769                    x, y, y_err = corr.plottable()
-770                    if hide_sigma:
-771                        hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
-772                    else:
-773                        hide_from = None
-774                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
-775            else:
-776                raise Exception("'comp' must be a correlator or a list of correlators.")
-777
-778        if plateau:
-779            if isinstance(plateau, Obs):
-780                if auto_gamma:
-781                    plateau.gamma_method()
-782                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
-783                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
-784            else:
-785                raise Exception("'plateau' must be an Obs")
-786
-787        if references:
-788            if isinstance(references, list):
-789                for ref in references:
-790                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
-791            else:
-792                raise Exception("'references' must be a list of floating pint values.")
-793
-794        if self.prange:
-795            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
-796            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+            
711    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):
+712        """Plots the correlator using the tag of the correlator as label if available.
+713
+714        Parameters
+715        ----------
+716        x_range : list
+717            list of two values, determining the range of the x-axis e.g. [4, 8]
+718        comp : Corr or list of Corr
+719            Correlator or list of correlators which are plotted for comparison.
+720            The tags of these correlators are used as labels if available.
+721        logscale : bool
+722            Sets y-axis to logscale
+723        plateau : Obs
+724            Plateau value to be visualized in the figure
+725        fit_res : Fit_result
+726            Fit_result object to be visualized
+727        ylabel : str
+728            Label for the y-axis
+729        save : str
+730            path to file in which the figure should be saved
+731        auto_gamma : bool
+732            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+733        hide_sigma : float
+734            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+735        references : list
+736            List of floating point values that are displayed as horizontal lines for reference.
+737        """
+738        if self.N != 1:
+739            raise Exception("Correlator must be projected before plotting")
+740
+741        if auto_gamma:
+742            self.gamma_method()
+743
+744        if x_range is None:
+745            x_range = [0, self.T - 1]
+746
+747        fig = plt.figure()
+748        ax1 = fig.add_subplot(111)
+749
+750        x, y, y_err = self.plottable()
+751        if hide_sigma:
+752            hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
+753        else:
+754            hide_from = None
+755        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+756        if logscale:
+757            ax1.set_yscale('log')
+758        else:
+759            if y_range is None:
+760                try:
+761                    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)])
+762                    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)])
+763                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+764                except Exception:
+765                    pass
+766            else:
+767                ax1.set_ylim(y_range)
+768        if comp:
+769            if isinstance(comp, (Corr, list)):
+770                for corr in comp if isinstance(comp, list) else [comp]:
+771                    if auto_gamma:
+772                        corr.gamma_method()
+773                    x, y, y_err = corr.plottable()
+774                    if hide_sigma:
+775                        hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
+776                    else:
+777                        hide_from = None
+778                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+779            else:
+780                raise Exception("'comp' must be a correlator or a list of correlators.")
+781
+782        if plateau:
+783            if isinstance(plateau, Obs):
+784                if auto_gamma:
+785                    plateau.gamma_method()
+786                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+787                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+788            else:
+789                raise Exception("'plateau' must be an Obs")
+790
+791        if references:
+792            if isinstance(references, list):
+793                for ref in references:
+794                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+795            else:
+796                raise Exception("'references' must be a list of floating pint values.")
 797
-798        if fit_res:
-799            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
-800            ax1.plot(x_samples,
-801                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
-802                     ls='-', marker=',', lw=2)
-803
-804        ax1.set_xlabel(r'$x_0 / a$')
-805        if ylabel:
-806            ax1.set_ylabel(ylabel)
-807        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
-808
-809        handles, labels = ax1.get_legend_handles_labels()
-810        if labels:
-811            ax1.legend()
-812        plt.draw()
-813
-814        if save:
-815            if isinstance(save, str):
-816                fig.savefig(save)
-817            else:
-818                raise Exception("'save' has to be a string.")
+798        if self.prange:
+799            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
+800            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+801
+802        if fit_res:
+803            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+804            ax1.plot(x_samples,
+805                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
+806                     ls='-', marker=',', lw=2)
+807
+808        ax1.set_xlabel(r'$x_0 / a$')
+809        if ylabel:
+810            ax1.set_ylabel(ylabel)
+811        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
+812
+813        handles, labels = ax1.get_legend_handles_labels()
+814        if labels:
+815            ax1.legend()
+816        plt.draw()
+817
+818        if save:
+819            if isinstance(save, str):
+820                fig.savefig(save)
+821            else:
+822                raise Exception("'save' has to be a string.")
 

 
         

@@ -4078,34 +4090,34 @@ List of floating point values that are displayed as horizontal lines for referen
 
             

             View Source
-            
820    def spaghetti_plot(self, logscale=True):
-821        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-822
-823        Parameters
-824        ----------
-825        logscale : bool
-826            Determines whether the scale of the y-axis is logarithmic or standard.
-827        """
-828        if self.N != 1:
-829            raise Exception("Correlator needs to be projected first.")
-830
-831        mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
-832        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
-833
-834        for name in mc_names:
-835            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
-836
-837            fig = plt.figure()
-838            ax = fig.add_subplot(111)
-839            for dat in data:
-840                ax.plot(x0_vals, dat, ls='-', marker='')
-841
-842            if logscale is True:
-843                ax.set_yscale('log')
-844
-845            ax.set_xlabel(r'$x_0 / a$')
-846            plt.title(name)
-847            plt.draw()
+            
824    def spaghetti_plot(self, logscale=True):
+825        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+826
+827        Parameters
+828        ----------
+829        logscale : bool
+830            Determines whether the scale of the y-axis is logarithmic or standard.
+831        """
+832        if self.N != 1:
+833            raise Exception("Correlator needs to be projected first.")
+834
+835        mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
+836        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
+837
+838        for name in mc_names:
+839            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+840
+841            fig = plt.figure()
+842            ax = fig.add_subplot(111)
+843            for dat in data:
+844                ax.plot(x0_vals, dat, ls='-', marker='')
+845
+846            if logscale is True:
+847                ax.set_yscale('log')
+848
+849            ax.set_xlabel(r'$x_0 / a$')
+850            plt.title(name)
+851            plt.draw()
 

 
         

@@ -4132,29 +4144,29 @@ Determines whether the scale of the y-axis is logarithmic or standard.
 
             

             View Source
-            
849    def dump(self, filename, datatype="json.gz", **kwargs):
-850        """Dumps the Corr into a file of chosen type
-851        Parameters
-852        ----------
-853        filename : str
-854            Name of the file to be saved.
-855        datatype : str
-856            Format of the exported file. Supported formats include
-857            "json.gz" and "pickle"
-858        path : str
-859            specifies a custom path for the file (default '.')
-860        """
-861        if datatype == "json.gz":
-862            from .input.json import dump_to_json
-863            if 'path' in kwargs:
-864                file_name = kwargs.get('path') + '/' + filename
-865            else:
-866                file_name = filename
-867            dump_to_json(self, file_name)
-868        elif datatype == "pickle":
-869            dump_object(self, filename, **kwargs)
-870        else:
-871            raise Exception("Unknown datatype " + str(datatype))
+            
853    def dump(self, filename, datatype="json.gz", **kwargs):
+854        """Dumps the Corr into a file of chosen type
+855        Parameters
+856        ----------
+857        filename : str
+858            Name of the file to be saved.
+859        datatype : str
+860            Format of the exported file. Supported formats include
+861            "json.gz" and "pickle"
+862        path : str
+863            specifies a custom path for the file (default '.')
+864        """
+865        if datatype == "json.gz":
+866            from .input.json import dump_to_json
+867            if 'path' in kwargs:
+868                file_name = kwargs.get('path') + '/' + filename
+869            else:
+870                file_name = filename
+871            dump_to_json(self, file_name)
+872        elif datatype == "pickle":
+873            dump_object(self, filename, **kwargs)
+874        else:
+875            raise Exception("Unknown datatype " + str(datatype))
 

 
         

@@ -4186,8 +4198,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
873    def print(self, range=[0, None]):
-874        print(self.__repr__(range))
+            
877    def print(self, range=[0, None]):
+878        print(self.__repr__(range))
 

 
         

@@ -4205,8 +4217,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1036    def sqrt(self):
-1037        return self**0.5
+            
1040    def sqrt(self):
+1041        return self**0.5
 

 
         

@@ -4224,9 +4236,9 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1039    def log(self):
-1040        newcontent = [None if (item is None) else np.log(item) for item in self.content]
-1041        return Corr(newcontent, prange=self.prange)
+            
1043    def log(self):
+1044        newcontent = [None if (item is None) else np.log(item) for item in self.content]
+1045        return Corr(newcontent, prange=self.prange)
 

 
         

@@ -4244,9 +4256,9 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1043    def exp(self):
-1044        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
-1045        return Corr(newcontent, prange=self.prange)
+            
1047    def exp(self):
+1048        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
+1049        return Corr(newcontent, prange=self.prange)
 

 
         

@@ -4264,8 +4276,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1058    def sin(self):
-1059        return self._apply_func_to_corr(np.sin)
+            
1062    def sin(self):
+1063        return self._apply_func_to_corr(np.sin)
 

 
         

@@ -4283,8 +4295,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1061    def cos(self):
-1062        return self._apply_func_to_corr(np.cos)
+            
1065    def cos(self):
+1066        return self._apply_func_to_corr(np.cos)
 

 
         

@@ -4302,8 +4314,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1064    def tan(self):
-1065        return self._apply_func_to_corr(np.tan)
+            
1068    def tan(self):
+1069        return self._apply_func_to_corr(np.tan)
 

 
         

@@ -4321,8 +4333,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1067    def sinh(self):
-1068        return self._apply_func_to_corr(np.sinh)
+            
1071    def sinh(self):
+1072        return self._apply_func_to_corr(np.sinh)
 

 
         

@@ -4340,8 +4352,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1070    def cosh(self):
-1071        return self._apply_func_to_corr(np.cosh)
+            
1074    def cosh(self):
+1075        return self._apply_func_to_corr(np.cosh)
 

 
         

@@ -4359,8 +4371,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1073    def tanh(self):
-1074        return self._apply_func_to_corr(np.tanh)
+            
1077    def tanh(self):
+1078        return self._apply_func_to_corr(np.tanh)
 

 
         

@@ -4378,8 +4390,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1076    def arcsin(self):
-1077        return self._apply_func_to_corr(np.arcsin)
+            
1080    def arcsin(self):
+1081        return self._apply_func_to_corr(np.arcsin)
 

 
         

@@ -4397,8 +4409,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1079    def arccos(self):
-1080        return self._apply_func_to_corr(np.arccos)
+            
1083    def arccos(self):
+1084        return self._apply_func_to_corr(np.arccos)
 

 
         

@@ -4416,8 +4428,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1082    def arctan(self):
-1083        return self._apply_func_to_corr(np.arctan)
+            
1086    def arctan(self):
+1087        return self._apply_func_to_corr(np.arctan)
 

 
         

@@ -4435,8 +4447,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1085    def arcsinh(self):
-1086        return self._apply_func_to_corr(np.arcsinh)
+            
1089    def arcsinh(self):
+1090        return self._apply_func_to_corr(np.arcsinh)
 

 
         

@@ -4454,8 +4466,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1088    def arccosh(self):
-1089        return self._apply_func_to_corr(np.arccosh)
+            
1092    def arccosh(self):
+1093        return self._apply_func_to_corr(np.arccosh)
 

 
         

@@ -4473,8 +4485,8 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1091    def arctanh(self):
-1092        return self._apply_func_to_corr(np.arctanh)
+            
1095    def arctanh(self):
+1096        return self._apply_func_to_corr(np.arctanh)
 

 
         

@@ -4512,64 +4524,64 @@ specifies a custom path for the file (default '.')
 
             

             View Source
-            
1127    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1128        r''' Project large correlation matrix to lowest states
-1129
-1130        This method can be used to reduce the size of an (N x N) correlation matrix
-1131        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1132        is still small.
+            
1131    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1132        r''' Project large correlation matrix to lowest states
 1133
-1134        Parameters
-1135        ----------
-1136        Ntrunc: int
-1137            Rank of the target matrix.
-1138        tproj: int
-1139            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1140            The default value is 3.
-1141        t0proj: int
-1142            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1143            discouraged for O(a) improved theories, since the correctness of the procedure
-1144            cannot be granted in this case. The default value is 2.
-1145        basematrix : Corr
-1146            Correlation matrix that is used to determine the eigenvectors of the
-1147            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1148            is is not specified.
-1149
-1150        Notes
-1151        -----
-1152        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1153        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}$
-1154        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1155        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1156        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1157        correlation matrix and to remove some noise that is added by irrelevant operators.
-1158        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1159        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1160        '''
-1161
-1162        if self.N == 1:
-1163            raise Exception('Method cannot be applied to one-dimensional correlators.')
-1164        if basematrix is None:
-1165            basematrix = self
-1166        if Ntrunc >= basematrix.N:
-1167            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
-1168        if basematrix.N != self.N:
-1169            raise Exception('basematrix and targetmatrix have to be of the same size.')
-1170
-1171        evecs = []
-1172        for i in range(Ntrunc):
-1173            evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
+1134        This method can be used to reduce the size of an (N x N) correlation matrix
+1135        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1136        is still small.
+1137
+1138        Parameters
+1139        ----------
+1140        Ntrunc: int
+1141            Rank of the target matrix.
+1142        tproj: int
+1143            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1144            The default value is 3.
+1145        t0proj: int
+1146            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1147            discouraged for O(a) improved theories, since the correctness of the procedure
+1148            cannot be granted in this case. The default value is 2.
+1149        basematrix : Corr
+1150            Correlation matrix that is used to determine the eigenvectors of the
+1151            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1152            is is not specified.
+1153
+1154        Notes
+1155        -----
+1156        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1157        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}$
+1158        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1159        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1160        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1161        correlation matrix and to remove some noise that is added by irrelevant operators.
+1162        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1163        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1164        '''
+1165
+1166        if self.N == 1:
+1167            raise Exception('Method cannot be applied to one-dimensional correlators.')
+1168        if basematrix is None:
+1169            basematrix = self
+1170        if Ntrunc >= basematrix.N:
+1171            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
+1172        if basematrix.N != self.N:
+1173            raise Exception('basematrix and targetmatrix have to be of the same size.')
 1174
-1175        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1176        rmat = []
-1177        for t in range(basematrix.T):
-1178            for i in range(Ntrunc):
-1179                for j in range(Ntrunc):
-1180                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1181            rmat.append(np.copy(tmpmat))
-1182
-1183        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
-1184        return Corr(newcontent)
+1175        evecs = []
+1176        for i in range(Ntrunc):
+1177            evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
+1178
+1179        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1180        rmat = []
+1181        for t in range(basematrix.T):
+1182            for i in range(Ntrunc):
+1183                for j in range(Ntrunc):
+1184                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1185            rmat.append(np.copy(tmpmat))
+1186
+1187        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
+1188        return Corr(newcontent)