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

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

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

             View Source
-            
1040    def sqrt(self):
-1041        return self**0.5
+            
1042    def sqrt(self):
+1043        return self**0.5
 

 
         

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

             View Source
-            
1043    def log(self):
-1044        newcontent = [None if (item is None) else np.log(item) for item in self.content]
-1045        return Corr(newcontent, prange=self.prange)
+            
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)
 

 
         

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

             View Source
-            
1047    def exp(self):
-1048        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
-1049        return Corr(newcontent, prange=self.prange)
+            
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)
 

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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

 
         

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

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