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

267 symmetric_corr = self.matrix_symmetric() 268 if sorted_list is None: 269 if (ts is None): -270 raise Exception("ts is required if sorted_list=None") -271 if (self.content[t0] is None) or (self.content[ts] is None): -272 raise Exception("Corr not defined at t0/ts") -273 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") -274 for i in range(self.N): -275 for j in range(self.N): -276 G0[i, j] = symmetric_corr[t0][i, j].value -277 Gt[i, j] = symmetric_corr[ts][i, j].value -278 -279 sp_vecs = _GEVP_solver(Gt, G0) -280 sp_vec = sp_vecs[state] -281 return sp_vec -282 elif sorted_list in ["Eigenvalue", "Eigenvector"]: -283 if sorted_list == "Eigenvalue" and ts is not None: -284 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) -285 all_vecs = [] -286 for t in range(self.T): -287 try: -288 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") -289 for i in range(self.N): -290 for j in range(self.N): -291 G0[i, j] = symmetric_corr[t0][i, j].value -292 Gt[i, j] = symmetric_corr[t][i, j].value -293 -294 sp_vecs = _GEVP_solver(Gt, G0) -295 if sorted_list == "Eigenvalue": -296 sp_vec = sp_vecs[state] -297 all_vecs.append(sp_vec) -298 else: -299 all_vecs.append(sp_vecs) -300 except Exception: -301 all_vecs.append(None) -302 if sorted_list == "Eigenvector": -303 if (ts is None): -304 raise Exception("ts is required for the Eigenvector sorting method.") -305 all_vecs = _sort_vectors(all_vecs, ts) -306 all_vecs = [a[state] for a in all_vecs] -307 else: -308 raise Exception("Unkown value for 'sorted_list'.") -309 -310 return all_vecs +270 raise Exception("ts is required if sorted_list=None.") +271 if (ts <= t0): +272 raise Exception("ts has to be larger than t0.") +273 if (self.content[t0] is None) or (self.content[ts] is None): +274 raise Exception("Corr not defined at t0/ts.") +275 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") +276 for i in range(self.N): +277 for j in range(self.N): +278 G0[i, j] = symmetric_corr[t0][i, j].value +279 Gt[i, j] = symmetric_corr[ts][i, j].value +280 +281 sp_vecs = _GEVP_solver(Gt, G0) +282 sp_vec = sp_vecs[state] +283 return sp_vec +284 elif sorted_list in ["Eigenvalue", "Eigenvector"]: +285 if sorted_list == "Eigenvalue" and ts is not None: +286 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) +287 all_vecs = [None] * (t0 + 1) +288 for t in range(t0 + 1, self.T): +289 try: +290 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") +291 for i in range(self.N): +292 for j in range(self.N): +293 G0[i, j] = symmetric_corr[t0][i, j].value +294 Gt[i, j] = symmetric_corr[t][i, j].value +295 +296 sp_vecs = _GEVP_solver(Gt, G0) +297 if sorted_list == "Eigenvalue": +298 sp_vec = sp_vecs[state] +299 all_vecs.append(sp_vec) +300 else: +301 all_vecs.append(sp_vecs) +302 except Exception: +303 all_vecs.append(None) +304 if sorted_list == "Eigenvector": +305 if (ts is None): +306 raise Exception("ts is required for the Eigenvector sorting method.") +307 all_vecs = _sort_vectors(all_vecs, ts) +308 all_vecs = [a[state] if a is not None else None for a in all_vecs] +309 else: +310 raise Exception("Unkown value for 'sorted_list'.") +311 +312 return all_vecs @@ -3140,26 +3146,26 @@ if this argument is set, a list of vectors (len=self.T) is returned. If it is le - 
312    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
-313        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
-314
-315        Parameters
-316        ----------
-317        t0 : int
-318            The time t0 for G(t)v= lambda G(t_0)v
-319        ts : int
-320            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
-321            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
-322        state : int
-323            The state one is interested in ordered by energy. The lowest state is zero.
-324        sorted_list : string
-325            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.
-326             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
-327             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
-328                              The reference state is identified by its eigenvalue at t=ts
-329        """
-330        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
-331        return self.projected(vec)
+            
314    def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
+315        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
+316
+317        Parameters
+318        ----------
+319        t0 : int
+320            The time t0 for G(t)v= lambda G(t_0)v
+321        ts : int
+322            fixed time G(t_s)v= lambda G(t_0)v  if return_list=False
+323            If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
+324        state : int
+325            The state one is interested in ordered by energy. The lowest state is zero.
+326        sorted_list : string
+327            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.
+328             "Eigenvalue"  -  The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
+329             "Eigenvector" -  Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
+330                              The reference state is identified by its eigenvalue at t=ts
+331        """
+332        vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
+333        return self.projected(vec)
 

 
 
@@ -3196,46 +3202,46 @@ if this argument is set, a list of vectors (len=self.T) is returned. If it is le
- 
333    def Hankel(self, N, periodic=False):
-334        """Constructs an NxN Hankel matrix
-335
-336        C(t) c(t+1) ... c(t+n-1)
-337        C(t+1) c(t+2) ... c(t+n)
-338        .................
-339        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
-340
-341        Parameters
-342        ----------
-343        N : int
-344            Dimension of the Hankel matrix
-345        periodic : bool, optional
-346            determines whether the matrix is extended periodically
-347        """
-348
-349        if self.N != 1:
-350            raise Exception("Multi-operator Prony not implemented!")
-351
-352        array = np.empty([N, N], dtype="object")
-353        new_content = []
-354        for t in range(self.T):
-355            new_content.append(array.copy())
-356
-357        def wrap(i):
-358            while i >= self.T:
-359                i -= self.T
-360            return i
-361
-362        for t in range(self.T):
-363            for i in range(N):
-364                for j in range(N):
-365                    if periodic:
-366                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
-367                    elif (t + i + j) >= self.T:
-368                        new_content[t] = None
-369                    else:
-370                        new_content[t][i, j] = self.content[t + i + j][0]
-371
-372        return Corr(new_content)
+            
335    def Hankel(self, N, periodic=False):
+336        """Constructs an NxN Hankel matrix
+337
+338        C(t) c(t+1) ... c(t+n-1)
+339        C(t+1) c(t+2) ... c(t+n)
+340        .................
+341        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
+342
+343        Parameters
+344        ----------
+345        N : int
+346            Dimension of the Hankel matrix
+347        periodic : bool, optional
+348            determines whether the matrix is extended periodically
+349        """
+350
+351        if self.N != 1:
+352            raise Exception("Multi-operator Prony not implemented!")
+353
+354        array = np.empty([N, N], dtype="object")
+355        new_content = []
+356        for t in range(self.T):
+357            new_content.append(array.copy())
+358
+359        def wrap(i):
+360            while i >= self.T:
+361                i -= self.T
+362            return i
+363
+364        for t in range(self.T):
+365            for i in range(N):
+366                for j in range(N):
+367                    if periodic:
+368                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
+369                    elif (t + i + j) >= self.T:
+370                        new_content[t] = None
+371                    else:
+372                        new_content[t][i, j] = self.content[t + i + j][0]
+373
+374        return Corr(new_content)
 

 
 
@@ -3269,15 +3275,15 @@ determines whether the matrix is extended periodically
- 
374    def roll(self, dt):
-375        """Periodically shift the correlator by dt timeslices
-376
-377        Parameters
-378        ----------
-379        dt : int
-380            number of timeslices
-381        """
-382        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
+            
376    def roll(self, dt):
+377        """Periodically shift the correlator by dt timeslices
+378
+379        Parameters
+380        ----------
+381        dt : int
+382            number of timeslices
+383        """
+384        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
 

 
 
@@ -3304,9 +3310,9 @@ number of timeslices
- 
384    def reverse(self):
-385        """Reverse the time ordering of the Corr"""
-386        return Corr(self.content[:: -1])
+            
386    def reverse(self):
+387        """Reverse the time ordering of the Corr"""
+388        return Corr(self.content[:: -1])
 

 
 
@@ -3326,23 +3332,23 @@ number of timeslices
- 
388    def thin(self, spacing=2, offset=0):
-389        """Thin out a correlator to suppress correlations
-390
-391        Parameters
-392        ----------
-393        spacing : int
-394            Keep only every 'spacing'th entry of the correlator
-395        offset : int
-396            Offset the equal spacing
-397        """
-398        new_content = []
-399        for t in range(self.T):
-400            if (offset + t) % spacing != 0:
-401                new_content.append(None)
-402            else:
-403                new_content.append(self.content[t])
-404        return Corr(new_content)
+            
390    def thin(self, spacing=2, offset=0):
+391        """Thin out a correlator to suppress correlations
+392
+393        Parameters
+394        ----------
+395        spacing : int
+396            Keep only every 'spacing'th entry of the correlator
+397        offset : int
+398            Offset the equal spacing
+399        """
+400        new_content = []
+401        for t in range(self.T):
+402            if (offset + t) % spacing != 0:
+403                new_content.append(None)
+404            else:
+405                new_content.append(self.content[t])
+406        return Corr(new_content)
 

 
 
@@ -3371,32 +3377,32 @@ Offset the equal spacing
- 
406    def correlate(self, partner):
-407        """Correlate the correlator with another correlator or Obs
-408
-409        Parameters
-410        ----------
-411        partner : Obs or Corr
-412            partner to correlate the correlator with.
-413            Can either be an Obs which is correlated with all entries of the
-414            correlator or a Corr of same length.
-415        """
-416        new_content = []
-417        for x0, t_slice in enumerate(self.content):
-418            if t_slice is None:
-419                new_content.append(None)
-420            else:
-421                if isinstance(partner, Corr):
-422                    if partner.content[x0] is None:
-423                        new_content.append(None)
-424                    else:
-425                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
-426                elif isinstance(partner, Obs):  # Should this include CObs?
-427                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
-428                else:
-429                    raise Exception("Can only correlate with an Obs or a Corr.")
-430
-431        return Corr(new_content)
+            
408    def correlate(self, partner):
+409        """Correlate the correlator with another correlator or Obs
+410
+411        Parameters
+412        ----------
+413        partner : Obs or Corr
+414            partner to correlate the correlator with.
+415            Can either be an Obs which is correlated with all entries of the
+416            correlator or a Corr of same length.
+417        """
+418        new_content = []
+419        for x0, t_slice in enumerate(self.content):
+420            if t_slice is None:
+421                new_content.append(None)
+422            else:
+423                if isinstance(partner, Corr):
+424                    if partner.content[x0] is None:
+425                        new_content.append(None)
+426                    else:
+427                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
+428                elif isinstance(partner, Obs):  # Should this include CObs?
+429                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
+430                else:
+431                    raise Exception("Can only correlate with an Obs or a Corr.")
+432
+433        return Corr(new_content)
 

 
 
@@ -3425,26 +3431,26 @@ correlator or a Corr of same length.
- 
433    def reweight(self, weight, **kwargs):
-434        """Reweight the correlator.
-435
-436        Parameters
-437        ----------
-438        weight : Obs
-439            Reweighting factor. An Observable that has to be defined on a superset of the
-440            configurations in obs[i].idl for all i.
-441        all_configs : bool
-442            if True, the reweighted observables are normalized by the average of
-443            the reweighting factor on all configurations in weight.idl and not
-444            on the configurations in obs[i].idl.
-445        """
-446        new_content = []
-447        for t_slice in self.content:
-448            if t_slice is None:
-449                new_content.append(None)
-450            else:
-451                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
-452        return Corr(new_content)
+            
435    def reweight(self, weight, **kwargs):
+436        """Reweight the correlator.
+437
+438        Parameters
+439        ----------
+440        weight : Obs
+441            Reweighting factor. An Observable that has to be defined on a superset of the
+442            configurations in obs[i].idl for all i.
+443        all_configs : bool
+444            if True, the reweighted observables are normalized by the average of
+445            the reweighting factor on all configurations in weight.idl and not
+446            on the configurations in obs[i].idl.
+447        """
+448        new_content = []
+449        for t_slice in self.content:
+450            if t_slice is None:
+451                new_content.append(None)
+452            else:
+453                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
+454        return Corr(new_content)
 

 
 
@@ -3476,33 +3482,33 @@ on the configurations in obs[i].idl.
- 
454    def T_symmetry(self, partner, parity=+1):
-455        """Return the time symmetry average of the correlator and its partner
-456
-457        Parameters
-458        ----------
-459        partner : Corr
-460            Time symmetry partner of the Corr
-461        partity : int
-462            Parity quantum number of the correlator, can be +1 or -1
-463        """
-464        if not isinstance(partner, Corr):
-465            raise Exception("T partner has to be a Corr object.")
-466        if parity not in [+1, -1]:
-467            raise Exception("Parity has to be +1 or -1.")
-468        T_partner = parity * partner.reverse()
-469
-470        t_slices = []
-471        test = (self - T_partner)
-472        test.gamma_method()
-473        for x0, t_slice in enumerate(test.content):
-474            if t_slice is not None:
-475                if not t_slice[0].is_zero_within_error(5):
-476                    t_slices.append(x0)
-477        if t_slices:
-478            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
-479
-480        return (self + T_partner) / 2
+            
456    def T_symmetry(self, partner, parity=+1):
+457        """Return the time symmetry average of the correlator and its partner
+458
+459        Parameters
+460        ----------
+461        partner : Corr
+462            Time symmetry partner of the Corr
+463        partity : int
+464            Parity quantum number of the correlator, can be +1 or -1
+465        """
+466        if not isinstance(partner, Corr):
+467            raise Exception("T partner has to be a Corr object.")
+468        if parity not in [+1, -1]:
+469            raise Exception("Parity has to be +1 or -1.")
+470        T_partner = parity * partner.reverse()
+471
+472        t_slices = []
+473        test = (self - T_partner)
+474        test.gamma_method()
+475        for x0, t_slice in enumerate(test.content):
+476            if t_slice is not None:
+477                if not t_slice[0].is_zero_within_error(5):
+478                    t_slices.append(x0)
+479        if t_slices:
+480            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
+481
+482        return (self + T_partner) / 2
 

 
 
@@ -3531,57 +3537,57 @@ Parity quantum number of the correlator, can be +1 or -1
- 
482    def deriv(self, variant="symmetric"):
-483        """Return the first derivative of the correlator with respect to x0.
-484
-485        Parameters
-486        ----------
-487        variant : str
-488            decides which definition of the finite differences derivative is used.
-489            Available choice: symmetric, forward, backward, improved, default: symmetric
-490        """
-491        if variant == "symmetric":
-492            newcontent = []
-493            for t in range(1, self.T - 1):
-494                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-495                    newcontent.append(None)
-496                else:
-497                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
-498            if(all([x is None for x in newcontent])):
-499                raise Exception('Derivative is undefined at all timeslices')
-500            return Corr(newcontent, padding=[1, 1])
-501        elif variant == "forward":
-502            newcontent = []
-503            for t in range(self.T - 1):
-504                if (self.content[t] is None) or (self.content[t + 1] is None):
-505                    newcontent.append(None)
-506                else:
-507                    newcontent.append(self.content[t + 1] - self.content[t])
-508            if(all([x is None for x in newcontent])):
-509                raise Exception("Derivative is undefined at all timeslices")
-510            return Corr(newcontent, padding=[0, 1])
-511        elif variant == "backward":
-512            newcontent = []
-513            for t in range(1, self.T):
-514                if (self.content[t - 1] is None) or (self.content[t] is None):
-515                    newcontent.append(None)
-516                else:
-517                    newcontent.append(self.content[t] - self.content[t - 1])
-518            if(all([x is None for x in newcontent])):
-519                raise Exception("Derivative is undefined at all timeslices")
-520            return Corr(newcontent, padding=[1, 0])
-521        elif variant == "improved":
-522            newcontent = []
-523            for t in range(2, self.T - 2):
-524                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):
-525                    newcontent.append(None)
-526                else:
-527                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
-528            if(all([x is None for x in newcontent])):
-529                raise Exception('Derivative is undefined at all timeslices')
-530            return Corr(newcontent, padding=[2, 2])
-531        else:
-532            raise Exception("Unknown variant.")
+            
484    def deriv(self, variant="symmetric"):
+485        """Return the first derivative of the correlator with respect to x0.
+486
+487        Parameters
+488        ----------
+489        variant : str
+490            decides which definition of the finite differences derivative is used.
+491            Available choice: symmetric, forward, backward, improved, default: symmetric
+492        """
+493        if variant == "symmetric":
+494            newcontent = []
+495            for t in range(1, self.T - 1):
+496                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+497                    newcontent.append(None)
+498                else:
+499                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
+500            if(all([x is None for x in newcontent])):
+501                raise Exception('Derivative is undefined at all timeslices')
+502            return Corr(newcontent, padding=[1, 1])
+503        elif variant == "forward":
+504            newcontent = []
+505            for t in range(self.T - 1):
+506                if (self.content[t] is None) or (self.content[t + 1] is None):
+507                    newcontent.append(None)
+508                else:
+509                    newcontent.append(self.content[t + 1] - self.content[t])
+510            if(all([x is None for x in newcontent])):
+511                raise Exception("Derivative is undefined at all timeslices")
+512            return Corr(newcontent, padding=[0, 1])
+513        elif variant == "backward":
+514            newcontent = []
+515            for t in range(1, self.T):
+516                if (self.content[t - 1] is None) or (self.content[t] is None):
+517                    newcontent.append(None)
+518                else:
+519                    newcontent.append(self.content[t] - self.content[t - 1])
+520            if(all([x is None for x in newcontent])):
+521                raise Exception("Derivative is undefined at all timeslices")
+522            return Corr(newcontent, padding=[1, 0])
+523        elif variant == "improved":
+524            newcontent = []
+525            for t in range(2, self.T - 2):
+526                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):
+527                    newcontent.append(None)
+528                else:
+529                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
+530            if(all([x is None for x in newcontent])):
+531                raise Exception('Derivative is undefined at all timeslices')
+532            return Corr(newcontent, padding=[2, 2])
+533        else:
+534            raise Exception("Unknown variant.")
 

 
 
@@ -3609,37 +3615,37 @@ Available choice: symmetric, forward, backward, improved, default: symmetric
     
-            
534    def second_deriv(self, variant="symmetric"):
-535        """Return the second derivative of the correlator with respect to x0.
-536
-537        Parameters
-538        ----------
-539        variant : str
-540            decides which definition of the finite differences derivative is used.
-541            Available choice: symmetric, improved, default: symmetric
-542        """
-543        if variant == "symmetric":
-544            newcontent = []
-545            for t in range(1, self.T - 1):
-546                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-547                    newcontent.append(None)
-548                else:
-549                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
-550            if(all([x is None for x in newcontent])):
-551                raise Exception("Derivative is undefined at all timeslices")
-552            return Corr(newcontent, padding=[1, 1])
-553        elif variant == "improved":
-554            newcontent = []
-555            for t in range(2, self.T - 2):
-556                if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
-557                    newcontent.append(None)
-558                else:
-559                    newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2]))
-560            if(all([x is None for x in newcontent])):
-561                raise Exception("Derivative is undefined at all timeslices")
-562            return Corr(newcontent, padding=[2, 2])
-563        else:
-564            raise Exception("Unknown variant.")
+            
536    def second_deriv(self, variant="symmetric"):
+537        """Return the second derivative of the correlator with respect to x0.
+538
+539        Parameters
+540        ----------
+541        variant : str
+542            decides which definition of the finite differences derivative is used.
+543            Available choice: symmetric, improved, default: symmetric
+544        """
+545        if variant == "symmetric":
+546            newcontent = []
+547            for t in range(1, self.T - 1):
+548                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+549                    newcontent.append(None)
+550                else:
+551                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
+552            if(all([x is None for x in newcontent])):
+553                raise Exception("Derivative is undefined at all timeslices")
+554            return Corr(newcontent, padding=[1, 1])
+555        elif variant == "improved":
+556            newcontent = []
+557            for t in range(2, self.T - 2):
+558                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):
+559                    newcontent.append(None)
+560                else:
+561                    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]))
+562            if(all([x is None for x in newcontent])):
+563                raise Exception("Derivative is undefined at all timeslices")
+564            return Corr(newcontent, padding=[2, 2])
+565        else:
+566            raise Exception("Unknown variant.")
 

 
 
@@ -3667,70 +3673,70 @@ Available choice: symmetric, improved, default: symmetric
 
     

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

 
 
@@ -3763,39 +3769,39 @@ guess for the root finder, only relevant for the root variant
 
     

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

 
 
@@ -3829,42 +3835,42 @@ Decides whether output is printed to the standard output.
 
     

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

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

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

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

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

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

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

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

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

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

     
-            
880    def print(self, range=[0, None]):
-881        print(self.__repr__(range))
+            
882    def print(self, range=[0, None]):
+883        print(self.__repr__(range))
 

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

     
-            
1043    def sqrt(self):
-1044        return self**0.5
+            
1045    def sqrt(self):
+1046        return self**0.5
 

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

     
-            
1046    def log(self):
-1047        newcontent = [None if (item is None) else np.log(item) for item in self.content]
-1048        return Corr(newcontent, prange=self.prange)
+            
1048    def log(self):
+1049        newcontent = [None if (item is None) else np.log(item) for item in self.content]
+1050        return Corr(newcontent, prange=self.prange)
 

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

     
-            
1050    def exp(self):
-1051        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
-1052        return Corr(newcontent, prange=self.prange)
+            
1052    def exp(self):
+1053        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
+1054        return Corr(newcontent, prange=self.prange)
 

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

     
-            
1065    def sin(self):
-1066        return self._apply_func_to_corr(np.sin)
+            
1067    def sin(self):
+1068        return self._apply_func_to_corr(np.sin)
 

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

     
-            
1068    def cos(self):
-1069        return self._apply_func_to_corr(np.cos)
+            
1070    def cos(self):
+1071        return self._apply_func_to_corr(np.cos)
 

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

     
-            
1071    def tan(self):
-1072        return self._apply_func_to_corr(np.tan)
+            
1073    def tan(self):
+1074        return self._apply_func_to_corr(np.tan)
 

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

     
-            
1074    def sinh(self):
-1075        return self._apply_func_to_corr(np.sinh)
+            
1076    def sinh(self):
+1077        return self._apply_func_to_corr(np.sinh)
 

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

     
-            
1077    def cosh(self):
-1078        return self._apply_func_to_corr(np.cosh)
+            
1079    def cosh(self):
+1080        return self._apply_func_to_corr(np.cosh)
 

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

     
-            
1080    def tanh(self):
-1081        return self._apply_func_to_corr(np.tanh)
+            
1082    def tanh(self):
+1083        return self._apply_func_to_corr(np.tanh)
 

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

     
-            
1083    def arcsin(self):
-1084        return self._apply_func_to_corr(np.arcsin)
+            
1085    def arcsin(self):
+1086        return self._apply_func_to_corr(np.arcsin)
 

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

     
-            
1086    def arccos(self):
-1087        return self._apply_func_to_corr(np.arccos)
+            
1088    def arccos(self):
+1089        return self._apply_func_to_corr(np.arccos)
 

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

     
-            
1089    def arctan(self):
-1090        return self._apply_func_to_corr(np.arctan)
+            
1091    def arctan(self):
+1092        return self._apply_func_to_corr(np.arctan)
 

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

     
-            
1092    def arcsinh(self):
-1093        return self._apply_func_to_corr(np.arcsinh)
+            
1094    def arcsinh(self):
+1095        return self._apply_func_to_corr(np.arcsinh)
 

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

     
-            
1095    def arccosh(self):
-1096        return self._apply_func_to_corr(np.arccosh)
+            
1097    def arccosh(self):
+1098        return self._apply_func_to_corr(np.arccosh)
 

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

     
-            
1098    def arctanh(self):
-1099        return self._apply_func_to_corr(np.arctanh)
+            
1100    def arctanh(self):
+1101        return self._apply_func_to_corr(np.arctanh)
 

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

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