305 def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
-306 """Determines the eigenvalue of the GEVP by solving and projecting the correlator
-307
-308 Parameters
-309 ----------
-310 t0 : int
-311 The time t0 for G(t)v= lambda G(t_0)v
-312 ts : int
-313 fixed time G(t_s)v= lambda G(t_0)v if return_list=False
-314 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
-315 state : int
-316 The state one is interested in ordered by energy. The lowest state is zero.
-317 sorted_list : string
-318 if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
-319 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
-320 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
-321 The reference state is identified by its eigenvalue at t=ts
-322 """
-323 vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
-324 return self.projected(vec)
+ 309 def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None):
+310 """Determines the eigenvalue of the GEVP by solving and projecting the correlator
+311
+312 Parameters
+313 ----------
+314 t0 : int
+315 The time t0 for G(t)v= lambda G(t_0)v
+316 ts : int
+317 fixed time G(t_s)v= lambda G(t_0)v if return_list=False
+318 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
+319 state : int
+320 The state one is interested in ordered by energy. The lowest state is zero.
+321 sorted_list : string
+322 if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned.
+323 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
+324 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state.
+325 The reference state is identified by its eigenvalue at t=ts
+326 """
+327 vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list)
+328 return self.projected(vec)
View Source
- 326 def Hankel(self, N, periodic=False):
-327 """Constructs an NxN Hankel matrix
-328
-329 C(t) c(t+1) ... c(t+n-1)
-330 C(t+1) c(t+2) ... c(t+n)
-331 .................
-332 C(t+(n-1)) c(t+n) ... c(t+2(n-1))
-333
-334 Parameters
-335 ----------
-336 N : int
-337 Dimension of the Hankel matrix
-338 periodic : bool, optional
-339 determines whether the matrix is extended periodically
-340 """
-341
-342 if self.N != 1:
-343 raise Exception("Multi-operator Prony not implemented!")
-344
-345 array = np.empty([N, N], dtype="object")
-346 new_content = []
-347 for t in range(self.T):
-348 new_content.append(array.copy())
-349
-350 def wrap(i):
-351 while i >= self.T:
-352 i -= self.T
-353 return i
-354
-355 for t in range(self.T):
-356 for i in range(N):
-357 for j in range(N):
-358 if periodic:
-359 new_content[t][i, j] = self.content[wrap(t + i + j)][0]
-360 elif (t + i + j) >= self.T:
-361 new_content[t] = None
-362 else:
-363 new_content[t][i, j] = self.content[t + i + j][0]
-364
-365 return Corr(new_content)
+ 330 def Hankel(self, N, periodic=False):
+331 """Constructs an NxN Hankel matrix
+332
+333 C(t) c(t+1) ... c(t+n-1)
+334 C(t+1) c(t+2) ... c(t+n)
+335 .................
+336 C(t+(n-1)) c(t+n) ... c(t+2(n-1))
+337
+338 Parameters
+339 ----------
+340 N : int
+341 Dimension of the Hankel matrix
+342 periodic : bool, optional
+343 determines whether the matrix is extended periodically
+344 """
+345
+346 if self.N != 1:
+347 raise Exception("Multi-operator Prony not implemented!")
+348
+349 array = np.empty([N, N], dtype="object")
+350 new_content = []
+351 for t in range(self.T):
+352 new_content.append(array.copy())
+353
+354 def wrap(i):
+355 while i >= self.T:
+356 i -= self.T
+357 return i
+358
+359 for t in range(self.T):
+360 for i in range(N):
+361 for j in range(N):
+362 if periodic:
+363 new_content[t][i, j] = self.content[wrap(t + i + j)][0]
+364 elif (t + i + j) >= self.T:
+365 new_content[t] = None
+366 else:
+367 new_content[t][i, j] = self.content[t + i + j][0]
+368
+369 return Corr(new_content)
View Source
- 367 def roll(self, dt):
-368 """Periodically shift the correlator by dt timeslices
-369
-370 Parameters
-371 ----------
-372 dt : int
-373 number of timeslices
-374 """
-375 return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
+ 371 def roll(self, dt):
+372 """Periodically shift the correlator by dt timeslices
+373
+374 Parameters
+375 ----------
+376 dt : int
+377 number of timeslices
+378 """
+379 return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
View Source
- 377 def reverse(self):
-378 """Reverse the time ordering of the Corr"""
-379 return Corr(self.content[:: -1])
+ 381 def reverse(self):
+382 """Reverse the time ordering of the Corr"""
+383 return Corr(self.content[:: -1])
View Source
- 381 def thin(self, spacing=2, offset=0):
-382 """Thin out a correlator to suppress correlations
-383
-384 Parameters
-385 ----------
-386 spacing : int
-387 Keep only every 'spacing'th entry of the correlator
-388 offset : int
-389 Offset the equal spacing
-390 """
-391 new_content = []
-392 for t in range(self.T):
-393 if (offset + t) % spacing != 0:
-394 new_content.append(None)
-395 else:
-396 new_content.append(self.content[t])
-397 return Corr(new_content)
+ 385 def thin(self, spacing=2, offset=0):
+386 """Thin out a correlator to suppress correlations
+387
+388 Parameters
+389 ----------
+390 spacing : int
+391 Keep only every 'spacing'th entry of the correlator
+392 offset : int
+393 Offset the equal spacing
+394 """
+395 new_content = []
+396 for t in range(self.T):
+397 if (offset + t) % spacing != 0:
+398 new_content.append(None)
+399 else:
+400 new_content.append(self.content[t])
+401 return Corr(new_content)
View Source
- 399 def correlate(self, partner):
-400 """Correlate the correlator with another correlator or Obs
-401
-402 Parameters
-403 ----------
-404 partner : Obs or Corr
-405 partner to correlate the correlator with.
-406 Can either be an Obs which is correlated with all entries of the
-407 correlator or a Corr of same length.
-408 """
-409 new_content = []
-410 for x0, t_slice in enumerate(self.content):
-411 if t_slice is None:
-412 new_content.append(None)
-413 else:
-414 if isinstance(partner, Corr):
-415 if partner.content[x0] is None:
-416 new_content.append(None)
-417 else:
-418 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
-419 elif isinstance(partner, Obs): # Should this include CObs?
-420 new_content.append(np.array([correlate(o, partner) for o in t_slice]))
-421 else:
-422 raise Exception("Can only correlate with an Obs or a Corr.")
-423
-424 return Corr(new_content)
+ 403 def correlate(self, partner):
+404 """Correlate the correlator with another correlator or Obs
+405
+406 Parameters
+407 ----------
+408 partner : Obs or Corr
+409 partner to correlate the correlator with.
+410 Can either be an Obs which is correlated with all entries of the
+411 correlator or a Corr of same length.
+412 """
+413 new_content = []
+414 for x0, t_slice in enumerate(self.content):
+415 if t_slice is None:
+416 new_content.append(None)
+417 else:
+418 if isinstance(partner, Corr):
+419 if partner.content[x0] is None:
+420 new_content.append(None)
+421 else:
+422 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
+423 elif isinstance(partner, Obs): # Should this include CObs?
+424 new_content.append(np.array([correlate(o, partner) for o in t_slice]))
+425 else:
+426 raise Exception("Can only correlate with an Obs or a Corr.")
+427
+428 return Corr(new_content)
View Source
- 426 def reweight(self, weight, **kwargs):
-427 """Reweight the correlator.
-428
-429 Parameters
-430 ----------
-431 weight : Obs
-432 Reweighting factor. An Observable that has to be defined on a superset of the
-433 configurations in obs[i].idl for all i.
-434 all_configs : bool
-435 if True, the reweighted observables are normalized by the average of
-436 the reweighting factor on all configurations in weight.idl and not
-437 on the configurations in obs[i].idl.
-438 """
-439 new_content = []
-440 for t_slice in self.content:
-441 if t_slice is None:
-442 new_content.append(None)
-443 else:
-444 new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
-445 return Corr(new_content)
+ 430 def reweight(self, weight, **kwargs):
+431 """Reweight the correlator.
+432
+433 Parameters
+434 ----------
+435 weight : Obs
+436 Reweighting factor. An Observable that has to be defined on a superset of the
+437 configurations in obs[i].idl for all i.
+438 all_configs : bool
+439 if True, the reweighted observables are normalized by the average of
+440 the reweighting factor on all configurations in weight.idl and not
+441 on the configurations in obs[i].idl.
+442 """
+443 new_content = []
+444 for t_slice in self.content:
+445 if t_slice is None:
+446 new_content.append(None)
+447 else:
+448 new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
+449 return Corr(new_content)
View Source
- 447 def T_symmetry(self, partner, parity=+1):
-448 """Return the time symmetry average of the correlator and its partner
-449
-450 Parameters
-451 ----------
-452 partner : Corr
-453 Time symmetry partner of the Corr
-454 partity : int
-455 Parity quantum number of the correlator, can be +1 or -1
-456 """
-457 if not isinstance(partner, Corr):
-458 raise Exception("T partner has to be a Corr object.")
-459 if parity not in [+1, -1]:
-460 raise Exception("Parity has to be +1 or -1.")
-461 T_partner = parity * partner.reverse()
-462
-463 t_slices = []
-464 test = (self - T_partner)
-465 test.gamma_method()
-466 for x0, t_slice in enumerate(test.content):
-467 if t_slice is not None:
-468 if not t_slice[0].is_zero_within_error(5):
-469 t_slices.append(x0)
-470 if t_slices:
-471 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
-472
-473 return (self + T_partner) / 2
+ 451 def T_symmetry(self, partner, parity=+1):
+452 """Return the time symmetry average of the correlator and its partner
+453
+454 Parameters
+455 ----------
+456 partner : Corr
+457 Time symmetry partner of the Corr
+458 partity : int
+459 Parity quantum number of the correlator, can be +1 or -1
+460 """
+461 if not isinstance(partner, Corr):
+462 raise Exception("T partner has to be a Corr object.")
+463 if parity not in [+1, -1]:
+464 raise Exception("Parity has to be +1 or -1.")
+465 T_partner = parity * partner.reverse()
+466
+467 t_slices = []
+468 test = (self - T_partner)
+469 test.gamma_method()
+470 for x0, t_slice in enumerate(test.content):
+471 if t_slice is not None:
+472 if not t_slice[0].is_zero_within_error(5):
+473 t_slices.append(x0)
+474 if t_slices:
+475 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
+476
+477 return (self + T_partner) / 2
View Source
- 475 def deriv(self, variant="symmetric"):
-476 """Return the first derivative of the correlator with respect to x0.
-477
-478 Parameters
-479 ----------
-480 variant : str
-481 decides which definition of the finite differences derivative is used.
-482 Available choice: symmetric, forward, backward, improved, default: symmetric
-483 """
-484 if variant == "symmetric":
-485 newcontent = []
-486 for t in range(1, self.T - 1):
-487 if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-488 newcontent.append(None)
-489 else:
-490 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
-491 if(all([x is None for x in newcontent])):
-492 raise Exception('Derivative is undefined at all timeslices')
-493 return Corr(newcontent, padding=[1, 1])
-494 elif variant == "forward":
-495 newcontent = []
-496 for t in range(self.T - 1):
-497 if (self.content[t] is None) or (self.content[t + 1] is None):
-498 newcontent.append(None)
-499 else:
-500 newcontent.append(self.content[t + 1] - self.content[t])
-501 if(all([x is None for x in newcontent])):
-502 raise Exception("Derivative is undefined at all timeslices")
-503 return Corr(newcontent, padding=[0, 1])
-504 elif variant == "backward":
-505 newcontent = []
-506 for t in range(1, self.T):
-507 if (self.content[t - 1] is None) or (self.content[t] is None):
-508 newcontent.append(None)
-509 else:
-510 newcontent.append(self.content[t] - self.content[t - 1])
-511 if(all([x is None for x in newcontent])):
-512 raise Exception("Derivative is undefined at all timeslices")
-513 return Corr(newcontent, padding=[1, 0])
-514 elif variant == "improved":
-515 newcontent = []
-516 for t in range(2, self.T - 2):
-517 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
-518 newcontent.append(None)
-519 else:
-520 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
-521 if(all([x is None for x in newcontent])):
-522 raise Exception('Derivative is undefined at all timeslices')
-523 return Corr(newcontent, padding=[2, 2])
-524 else:
-525 raise Exception("Unknown variant.")
+ 479 def deriv(self, variant="symmetric"):
+480 """Return the first derivative of the correlator with respect to x0.
+481
+482 Parameters
+483 ----------
+484 variant : str
+485 decides which definition of the finite differences derivative is used.
+486 Available choice: symmetric, forward, backward, improved, default: symmetric
+487 """
+488 if variant == "symmetric":
+489 newcontent = []
+490 for t in range(1, self.T - 1):
+491 if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+492 newcontent.append(None)
+493 else:
+494 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
+495 if(all([x is None for x in newcontent])):
+496 raise Exception('Derivative is undefined at all timeslices')
+497 return Corr(newcontent, padding=[1, 1])
+498 elif variant == "forward":
+499 newcontent = []
+500 for t in range(self.T - 1):
+501 if (self.content[t] is None) or (self.content[t + 1] is None):
+502 newcontent.append(None)
+503 else:
+504 newcontent.append(self.content[t + 1] - self.content[t])
+505 if(all([x is None for x in newcontent])):
+506 raise Exception("Derivative is undefined at all timeslices")
+507 return Corr(newcontent, padding=[0, 1])
+508 elif variant == "backward":
+509 newcontent = []
+510 for t in range(1, self.T):
+511 if (self.content[t - 1] is None) or (self.content[t] is None):
+512 newcontent.append(None)
+513 else:
+514 newcontent.append(self.content[t] - self.content[t - 1])
+515 if(all([x is None for x in newcontent])):
+516 raise Exception("Derivative is undefined at all timeslices")
+517 return Corr(newcontent, padding=[1, 0])
+518 elif variant == "improved":
+519 newcontent = []
+520 for t in range(2, self.T - 2):
+521 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
+522 newcontent.append(None)
+523 else:
+524 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
+525 if(all([x is None for x in newcontent])):
+526 raise Exception('Derivative is undefined at all timeslices')
+527 return Corr(newcontent, padding=[2, 2])
+528 else:
+529 raise Exception("Unknown variant.")
View Source
- 527 def second_deriv(self, variant="symmetric"):
-528 """Return the second derivative of the correlator with respect to x0.
-529
-530 Parameters
-531 ----------
-532 variant : str
-533 decides which definition of the finite differences derivative is used.
-534 Available choice: symmetric, improved, default: symmetric
-535 """
-536 if variant == "symmetric":
-537 newcontent = []
-538 for t in range(1, self.T - 1):
-539 if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-540 newcontent.append(None)
-541 else:
-542 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
-543 if(all([x is None for x in newcontent])):
-544 raise Exception("Derivative is undefined at all timeslices")
-545 return Corr(newcontent, padding=[1, 1])
-546 elif variant == "improved":
-547 newcontent = []
-548 for t in range(2, self.T - 2):
-549 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
-550 newcontent.append(None)
-551 else:
-552 newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2]))
-553 if(all([x is None for x in newcontent])):
-554 raise Exception("Derivative is undefined at all timeslices")
-555 return Corr(newcontent, padding=[2, 2])
-556 else:
-557 raise Exception("Unknown variant.")
+ 531 def second_deriv(self, variant="symmetric"):
+532 """Return the second derivative of the correlator with respect to x0.
+533
+534 Parameters
+535 ----------
+536 variant : str
+537 decides which definition of the finite differences derivative is used.
+538 Available choice: symmetric, improved, default: symmetric
+539 """
+540 if variant == "symmetric":
+541 newcontent = []
+542 for t in range(1, self.T - 1):
+543 if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+544 newcontent.append(None)
+545 else:
+546 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
+547 if(all([x is None for x in newcontent])):
+548 raise Exception("Derivative is undefined at all timeslices")
+549 return Corr(newcontent, padding=[1, 1])
+550 elif variant == "improved":
+551 newcontent = []
+552 for t in range(2, self.T - 2):
+553 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
+554 newcontent.append(None)
+555 else:
+556 newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2]))
+557 if(all([x is None for x in newcontent])):
+558 raise Exception("Derivative is undefined at all timeslices")
+559 return Corr(newcontent, padding=[2, 2])
+560 else:
+561 raise Exception("Unknown variant.")
View Source
- 559 def m_eff(self, variant='log', guess=1.0):
-560 """Returns the effective mass of the correlator as correlator object
-561
-562 Parameters
-563 ----------
-564 variant : str
-565 log : uses the standard effective mass log(C(t) / C(t+1))
-566 cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
-567 sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
-568 See, e.g., arXiv:1205.5380
-569 arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
-570 guess : float
-571 guess for the root finder, only relevant for the root variant
-572 """
-573 if self.N != 1:
-574 raise Exception('Correlator must be projected before getting m_eff')
-575 if variant == 'log':
-576 newcontent = []
-577 for t in range(self.T - 1):
-578 if (self.content[t] is None) or (self.content[t + 1] is None):
-579 newcontent.append(None)
-580 else:
-581 newcontent.append(self.content[t] / self.content[t + 1])
-582 if(all([x is None for x in newcontent])):
-583 raise Exception('m_eff is undefined at all timeslices')
-584
-585 return np.log(Corr(newcontent, padding=[0, 1]))
-586
-587 elif variant in ['periodic', 'cosh', 'sinh']:
-588 if variant in ['periodic', 'cosh']:
-589 func = anp.cosh
-590 else:
-591 func = anp.sinh
-592
-593 def root_function(x, d):
-594 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
-595
-596 newcontent = []
-597 for t in range(self.T - 1):
-598 if (self.content[t] is None) or (self.content[t + 1] is None):
-599 newcontent.append(None)
-600 # Fill the two timeslices in the middle of the lattice with their predecessors
-601 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
-602 newcontent.append(newcontent[-1])
-603 else:
-604 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
-605 if(all([x is None for x in newcontent])):
-606 raise Exception('m_eff is undefined at all timeslices')
-607
-608 return Corr(newcontent, padding=[0, 1])
-609
-610 elif variant == 'arccosh':
-611 newcontent = []
-612 for t in range(1, self.T - 1):
-613 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
-614 newcontent.append(None)
-615 else:
-616 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
-617 if(all([x is None for x in newcontent])):
-618 raise Exception("m_eff is undefined at all timeslices")
-619 return np.arccosh(Corr(newcontent, padding=[1, 1]))
-620
-621 else:
-622 raise Exception('Unknown variant.')
+ 563 def m_eff(self, variant='log', guess=1.0):
+564 """Returns the effective mass of the correlator as correlator object
+565
+566 Parameters
+567 ----------
+568 variant : str
+569 log : uses the standard effective mass log(C(t) / C(t+1))
+570 cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
+571 sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
+572 See, e.g., arXiv:1205.5380
+573 arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
+574 guess : float
+575 guess for the root finder, only relevant for the root variant
+576 """
+577 if self.N != 1:
+578 raise Exception('Correlator must be projected before getting m_eff')
+579 if variant == 'log':
+580 newcontent = []
+581 for t in range(self.T - 1):
+582 if (self.content[t] is None) or (self.content[t + 1] is None):
+583 newcontent.append(None)
+584 else:
+585 newcontent.append(self.content[t] / self.content[t + 1])
+586 if(all([x is None for x in newcontent])):
+587 raise Exception('m_eff is undefined at all timeslices')
+588
+589 return np.log(Corr(newcontent, padding=[0, 1]))
+590
+591 elif variant in ['periodic', 'cosh', 'sinh']:
+592 if variant in ['periodic', 'cosh']:
+593 func = anp.cosh
+594 else:
+595 func = anp.sinh
+596
+597 def root_function(x, d):
+598 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
+599
+600 newcontent = []
+601 for t in range(self.T - 1):
+602 if (self.content[t] is None) or (self.content[t + 1] is None):
+603 newcontent.append(None)
+604 # Fill the two timeslices in the middle of the lattice with their predecessors
+605 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+606 newcontent.append(newcontent[-1])
+607 else:
+608 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+609 if(all([x is None for x in newcontent])):
+610 raise Exception('m_eff is undefined at all timeslices')
+611
+612 return Corr(newcontent, padding=[0, 1])
+613
+614 elif variant == 'arccosh':
+615 newcontent = []
+616 for t in range(1, self.T - 1):
+617 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
+618 newcontent.append(None)
+619 else:
+620 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+621 if(all([x is None for x in newcontent])):
+622 raise Exception("m_eff is undefined at all timeslices")
+623 return np.arccosh(Corr(newcontent, padding=[1, 1]))
+624
+625 else:
+626 raise Exception('Unknown variant.')
View Source
- 624 def fit(self, function, fitrange=None, silent=False, **kwargs):
-625 r'''Fits function to the data
-626
-627 Parameters
-628 ----------
-629 function : obj
-630 function to fit to the data. See fits.least_squares for details.
-631 fitrange : list
-632 Two element list containing the timeslices on which the fit is supposed to start and stop.
-633 Caution: This range is inclusive as opposed to standard python indexing.
-634 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
-635 If not specified, self.prange or all timeslices are used.
-636 silent : bool
-637 Decides whether output is printed to the standard output.
-638 '''
-639 if self.N != 1:
-640 raise Exception("Correlator must be projected before fitting")
-641
-642 if fitrange is None:
-643 if self.prange:
-644 fitrange = self.prange
-645 else:
-646 fitrange = [0, self.T - 1]
-647 else:
-648 if not isinstance(fitrange, list):
-649 raise Exception("fitrange has to be a list with two elements")
-650 if len(fitrange) != 2:
-651 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
-652
-653 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-654 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-655 result = least_squares(xs, ys, function, silent=silent, **kwargs)
-656 return result
+ 628 def fit(self, function, fitrange=None, silent=False, **kwargs):
+629 r'''Fits function to the data
+630
+631 Parameters
+632 ----------
+633 function : obj
+634 function to fit to the data. See fits.least_squares for details.
+635 fitrange : list
+636 Two element list containing the timeslices on which the fit is supposed to start and stop.
+637 Caution: This range is inclusive as opposed to standard python indexing.
+638 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+639 If not specified, self.prange or all timeslices are used.
+640 silent : bool
+641 Decides whether output is printed to the standard output.
+642 '''
+643 if self.N != 1:
+644 raise Exception("Correlator must be projected before fitting")
+645
+646 if fitrange is None:
+647 if self.prange:
+648 fitrange = self.prange
+649 else:
+650 fitrange = [0, self.T - 1]
+651 else:
+652 if not isinstance(fitrange, list):
+653 raise Exception("fitrange has to be a list with two elements")
+654 if len(fitrange) != 2:
+655 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
+656
+657 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+658 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+659 result = least_squares(xs, ys, function, silent=silent, **kwargs)
+660 return result
View Source
- 658 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
-659 """ Extract a plateau value from a Corr object
-660
-661 Parameters
-662 ----------
-663 plateau_range : list
-664 list with two entries, indicating the first and the last timeslice
-665 of the plateau region.
-666 method : str
-667 method to extract the plateau.
-668 'fit' fits a constant to the plateau region
-669 'avg', 'average' or 'mean' just average over the given timeslices.
-670 auto_gamma : bool
-671 apply gamma_method with default parameters to the Corr. Defaults to None
-672 """
-673 if not plateau_range:
-674 if self.prange:
-675 plateau_range = self.prange
-676 else:
-677 raise Exception("no plateau range provided")
-678 if self.N != 1:
-679 raise Exception("Correlator must be projected before getting a plateau.")
-680 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
-681 raise Exception("plateau is undefined at all timeslices in plateaurange.")
-682 if auto_gamma:
-683 self.gamma_method()
-684 if method == "fit":
-685 def const_func(a, t):
-686 return a[0]
-687 return self.fit(const_func, plateau_range)[0]
-688 elif method in ["avg", "average", "mean"]:
-689 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
-690 return returnvalue
-691
-692 else:
-693 raise Exception("Unsupported plateau method: " + method)
+ 662 def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+663 """ Extract a plateau value from a Corr object
+664
+665 Parameters
+666 ----------
+667 plateau_range : list
+668 list with two entries, indicating the first and the last timeslice
+669 of the plateau region.
+670 method : str
+671 method to extract the plateau.
+672 'fit' fits a constant to the plateau region
+673 'avg', 'average' or 'mean' just average over the given timeslices.
+674 auto_gamma : bool
+675 apply gamma_method with default parameters to the Corr. Defaults to None
+676 """
+677 if not plateau_range:
+678 if self.prange:
+679 plateau_range = self.prange
+680 else:
+681 raise Exception("no plateau range provided")
+682 if self.N != 1:
+683 raise Exception("Correlator must be projected before getting a plateau.")
+684 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+685 raise Exception("plateau is undefined at all timeslices in plateaurange.")
+686 if auto_gamma:
+687 self.gamma_method()
+688 if method == "fit":
+689 def const_func(a, t):
+690 return a[0]
+691 return self.fit(const_func, plateau_range)[0]
+692 elif method in ["avg", "average", "mean"]:
+693 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+694 return returnvalue
+695
+696 else:
+697 raise Exception("Unsupported plateau method: " + method)
View Source
- 695 def set_prange(self, prange):
-696 """Sets the attribute prange of the Corr object."""
-697 if not len(prange) == 2:
-698 raise Exception("prange must be a list or array with two values")
-699 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
-700 raise Exception("Start and end point must be integers")
-701 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
-702 raise Exception("Start and end point must define a range in the interval 0,T")
-703
-704 self.prange = prange
-705 return
+ 699 def set_prange(self, prange):
+700 """Sets the attribute prange of the Corr object."""
+701 if not len(prange) == 2:
+702 raise Exception("prange must be a list or array with two values")
+703 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+704 raise Exception("Start and end point must be integers")
+705 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
+706 raise Exception("Start and end point must define a range in the interval 0,T")
+707
+708 self.prange = prange
+709 return
View Source
- 707 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None):
-708 """Plots the correlator using the tag of the correlator as label if available.
-709
-710 Parameters
-711 ----------
-712 x_range : list
-713 list of two values, determining the range of the x-axis e.g. [4, 8]
-714 comp : Corr or list of Corr
-715 Correlator or list of correlators which are plotted for comparison.
-716 The tags of these correlators are used as labels if available.
-717 logscale : bool
-718 Sets y-axis to logscale
-719 plateau : Obs
-720 Plateau value to be visualized in the figure
-721 fit_res : Fit_result
-722 Fit_result object to be visualized
-723 ylabel : str
-724 Label for the y-axis
-725 save : str
-726 path to file in which the figure should be saved
-727 auto_gamma : bool
-728 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
-729 hide_sigma : float
-730 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
-731 references : list
-732 List of floating point values that are displayed as horizontal lines for reference.
-733 """
-734 if self.N != 1:
-735 raise Exception("Correlator must be projected before plotting")
-736
-737 if auto_gamma:
-738 self.gamma_method()
-739
-740 if x_range is None:
-741 x_range = [0, self.T - 1]
-742
-743 fig = plt.figure()
-744 ax1 = fig.add_subplot(111)
-745
-746 x, y, y_err = self.plottable()
-747 if hide_sigma:
-748 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
-749 else:
-750 hide_from = None
-751 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
-752 if logscale:
-753 ax1.set_yscale('log')
-754 else:
-755 if y_range is None:
-756 try:
-757 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
-758 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
-759 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
-760 except Exception:
-761 pass
-762 else:
-763 ax1.set_ylim(y_range)
-764 if comp:
-765 if isinstance(comp, (Corr, list)):
-766 for corr in comp if isinstance(comp, list) else [comp]:
-767 if auto_gamma:
-768 corr.gamma_method()
-769 x, y, y_err = corr.plottable()
-770 if hide_sigma:
-771 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
-772 else:
-773 hide_from = None
-774 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
-775 else:
-776 raise Exception("'comp' must be a correlator or a list of correlators.")
-777
-778 if plateau:
-779 if isinstance(plateau, Obs):
-780 if auto_gamma:
-781 plateau.gamma_method()
-782 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
-783 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
-784 else:
-785 raise Exception("'plateau' must be an Obs")
-786
-787 if references:
-788 if isinstance(references, list):
-789 for ref in references:
-790 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
-791 else:
-792 raise Exception("'references' must be a list of floating pint values.")
-793
-794 if self.prange:
-795 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
-796 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+ 711 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None):
+712 """Plots the correlator using the tag of the correlator as label if available.
+713
+714 Parameters
+715 ----------
+716 x_range : list
+717 list of two values, determining the range of the x-axis e.g. [4, 8]
+718 comp : Corr or list of Corr
+719 Correlator or list of correlators which are plotted for comparison.
+720 The tags of these correlators are used as labels if available.
+721 logscale : bool
+722 Sets y-axis to logscale
+723 plateau : Obs
+724 Plateau value to be visualized in the figure
+725 fit_res : Fit_result
+726 Fit_result object to be visualized
+727 ylabel : str
+728 Label for the y-axis
+729 save : str
+730 path to file in which the figure should be saved
+731 auto_gamma : bool
+732 Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+733 hide_sigma : float
+734 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+735 references : list
+736 List of floating point values that are displayed as horizontal lines for reference.
+737 """
+738 if self.N != 1:
+739 raise Exception("Correlator must be projected before plotting")
+740
+741 if auto_gamma:
+742 self.gamma_method()
+743
+744 if x_range is None:
+745 x_range = [0, self.T - 1]
+746
+747 fig = plt.figure()
+748 ax1 = fig.add_subplot(111)
+749
+750 x, y, y_err = self.plottable()
+751 if hide_sigma:
+752 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
+753 else:
+754 hide_from = None
+755 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+756 if logscale:
+757 ax1.set_yscale('log')
+758 else:
+759 if y_range is None:
+760 try:
+761 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+762 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+763 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+764 except Exception:
+765 pass
+766 else:
+767 ax1.set_ylim(y_range)
+768 if comp:
+769 if isinstance(comp, (Corr, list)):
+770 for corr in comp if isinstance(comp, list) else [comp]:
+771 if auto_gamma:
+772 corr.gamma_method()
+773 x, y, y_err = corr.plottable()
+774 if hide_sigma:
+775 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1
+776 else:
+777 hide_from = None
+778 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+779 else:
+780 raise Exception("'comp' must be a correlator or a list of correlators.")
+781
+782 if plateau:
+783 if isinstance(plateau, Obs):
+784 if auto_gamma:
+785 plateau.gamma_method()
+786 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+787 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+788 else:
+789 raise Exception("'plateau' must be an Obs")
+790
+791 if references:
+792 if isinstance(references, list):
+793 for ref in references:
+794 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+795 else:
+796 raise Exception("'references' must be a list of floating pint values.")
797
-798 if fit_res:
-799 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
-800 ax1.plot(x_samples,
-801 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
-802 ls='-', marker=',', lw=2)
-803
-804 ax1.set_xlabel(r'$x_0 / a$')
-805 if ylabel:
-806 ax1.set_ylabel(ylabel)
-807 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
-808
-809 handles, labels = ax1.get_legend_handles_labels()
-810 if labels:
-811 ax1.legend()
-812 plt.draw()
-813
-814 if save:
-815 if isinstance(save, str):
-816 fig.savefig(save)
-817 else:
-818 raise Exception("'save' has to be a string.")
+798 if self.prange:
+799 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
+800 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+801
+802 if fit_res:
+803 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+804 ax1.plot(x_samples,
+805 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
+806 ls='-', marker=',', lw=2)
+807
+808 ax1.set_xlabel(r'$x_0 / a$')
+809 if ylabel:
+810 ax1.set_ylabel(ylabel)
+811 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
+812
+813 handles, labels = ax1.get_legend_handles_labels()
+814 if labels:
+815 ax1.legend()
+816 plt.draw()
+817
+818 if save:
+819 if isinstance(save, str):
+820 fig.savefig(save)
+821 else:
+822 raise Exception("'save' has to be a string.")
View Source
- 820 def spaghetti_plot(self, logscale=True):
-821 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-822
-823 Parameters
-824 ----------
-825 logscale : bool
-826 Determines whether the scale of the y-axis is logarithmic or standard.
-827 """
-828 if self.N != 1:
-829 raise Exception("Correlator needs to be projected first.")
-830
-831 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
-832 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
-833
-834 for name in mc_names:
-835 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
-836
-837 fig = plt.figure()
-838 ax = fig.add_subplot(111)
-839 for dat in data:
-840 ax.plot(x0_vals, dat, ls='-', marker='')
-841
-842 if logscale is True:
-843 ax.set_yscale('log')
-844
-845 ax.set_xlabel(r'$x_0 / a$')
-846 plt.title(name)
-847 plt.draw()
+ 824 def spaghetti_plot(self, logscale=True):
+825 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
+826
+827 Parameters
+828 ----------
+829 logscale : bool
+830 Determines whether the scale of the y-axis is logarithmic or standard.
+831 """
+832 if self.N != 1:
+833 raise Exception("Correlator needs to be projected first.")
+834
+835 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
+836 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
+837
+838 for name in mc_names:
+839 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+840
+841 fig = plt.figure()
+842 ax = fig.add_subplot(111)
+843 for dat in data:
+844 ax.plot(x0_vals, dat, ls='-', marker='')
+845
+846 if logscale is True:
+847 ax.set_yscale('log')
+848
+849 ax.set_xlabel(r'$x_0 / a$')
+850 plt.title(name)
+851 plt.draw()
View Source
- 849 def dump(self, filename, datatype="json.gz", **kwargs):
-850 """Dumps the Corr into a file of chosen type
-851 Parameters
-852 ----------
-853 filename : str
-854 Name of the file to be saved.
-855 datatype : str
-856 Format of the exported file. Supported formats include
-857 "json.gz" and "pickle"
-858 path : str
-859 specifies a custom path for the file (default '.')
-860 """
-861 if datatype == "json.gz":
-862 from .input.json import dump_to_json
-863 if 'path' in kwargs:
-864 file_name = kwargs.get('path') + '/' + filename
-865 else:
-866 file_name = filename
-867 dump_to_json(self, file_name)
-868 elif datatype == "pickle":
-869 dump_object(self, filename, **kwargs)
-870 else:
-871 raise Exception("Unknown datatype " + str(datatype))
+ 853 def dump(self, filename, datatype="json.gz", **kwargs):
+854 """Dumps the Corr into a file of chosen type
+855 Parameters
+856 ----------
+857 filename : str
+858 Name of the file to be saved.
+859 datatype : str
+860 Format of the exported file. Supported formats include
+861 "json.gz" and "pickle"
+862 path : str
+863 specifies a custom path for the file (default '.')
+864 """
+865 if datatype == "json.gz":
+866 from .input.json import dump_to_json
+867 if 'path' in kwargs:
+868 file_name = kwargs.get('path') + '/' + filename
+869 else:
+870 file_name = filename
+871 dump_to_json(self, file_name)
+872 elif datatype == "pickle":
+873 dump_object(self, filename, **kwargs)
+874 else:
+875 raise Exception("Unknown datatype " + str(datatype))
View Source
- 873 def print(self, range=[0, None]):
-874 print(self.__repr__(range))
+ 877 def print(self, range=[0, None]):
+878 print(self.__repr__(range))
View Source
- 1036 def sqrt(self):
-1037 return self**0.5
+ 1040 def sqrt(self):
+1041 return self**0.5
View Source
- 1039 def log(self):
-1040 newcontent = [None if (item is None) else np.log(item) for item in self.content]
-1041 return Corr(newcontent, prange=self.prange)
+ 1043 def log(self):
+1044 newcontent = [None if (item is None) else np.log(item) for item in self.content]
+1045 return Corr(newcontent, prange=self.prange)
View Source
- 1043 def exp(self):
-1044 newcontent = [None if (item is None) else np.exp(item) for item in self.content]
-1045 return Corr(newcontent, prange=self.prange)
+ 1047 def exp(self):
+1048 newcontent = [None if (item is None) else np.exp(item) for item in self.content]
+1049 return Corr(newcontent, prange=self.prange)
View Source
- 1058 def sin(self):
-1059 return self._apply_func_to_corr(np.sin)
+ 1062 def sin(self):
+1063 return self._apply_func_to_corr(np.sin)
View Source
- 1061 def cos(self):
-1062 return self._apply_func_to_corr(np.cos)
+ 1065 def cos(self):
+1066 return self._apply_func_to_corr(np.cos)
View Source
- 1064 def tan(self):
-1065 return self._apply_func_to_corr(np.tan)
+ 1068 def tan(self):
+1069 return self._apply_func_to_corr(np.tan)
View Source
- 1067 def sinh(self):
-1068 return self._apply_func_to_corr(np.sinh)
+ 1071 def sinh(self):
+1072 return self._apply_func_to_corr(np.sinh)
View Source
- 1070 def cosh(self):
-1071 return self._apply_func_to_corr(np.cosh)
+ 1074 def cosh(self):
+1075 return self._apply_func_to_corr(np.cosh)
View Source
- 1073 def tanh(self):
-1074 return self._apply_func_to_corr(np.tanh)
+ 1077 def tanh(self):
+1078 return self._apply_func_to_corr(np.tanh)
View Source
- 1076 def arcsin(self):
-1077 return self._apply_func_to_corr(np.arcsin)
+ 1080 def arcsin(self):
+1081 return self._apply_func_to_corr(np.arcsin)
View Source
- 1079 def arccos(self):
-1080 return self._apply_func_to_corr(np.arccos)
+ 1083 def arccos(self):
+1084 return self._apply_func_to_corr(np.arccos)
View Source
- 1082 def arctan(self):
-1083 return self._apply_func_to_corr(np.arctan)
+ 1086 def arctan(self):
+1087 return self._apply_func_to_corr(np.arctan)
View Source
- 1085 def arcsinh(self):
-1086 return self._apply_func_to_corr(np.arcsinh)
+ 1089 def arcsinh(self):
+1090 return self._apply_func_to_corr(np.arcsinh)
View Source
- 1088 def arccosh(self):
-1089 return self._apply_func_to_corr(np.arccosh)
+ 1092 def arccosh(self):
+1093 return self._apply_func_to_corr(np.arccosh)
View Source
- 1091 def arctanh(self):
-1092 return self._apply_func_to_corr(np.arctanh)
+ 1095 def arctanh(self):
+1096 return self._apply_func_to_corr(np.arctanh)
View Source
- 1127 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1128 r''' Project large correlation matrix to lowest states
-1129
-1130 This method can be used to reduce the size of an (N x N) correlation matrix
-1131 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1132 is still small.
+ 1131 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1132 r''' Project large correlation matrix to lowest states
1133
-1134 Parameters
-1135 ----------
-1136 Ntrunc: int
-1137 Rank of the target matrix.
-1138 tproj: int
-1139 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1140 The default value is 3.
-1141 t0proj: int
-1142 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1143 discouraged for O(a) improved theories, since the correctness of the procedure
-1144 cannot be granted in this case. The default value is 2.
-1145 basematrix : Corr
-1146 Correlation matrix that is used to determine the eigenvectors of the
-1147 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1148 is is not specified.
-1149
-1150 Notes
-1151 -----
-1152 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1153 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
-1154 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1155 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1156 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1157 correlation matrix and to remove some noise that is added by irrelevant operators.
-1158 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1159 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1160 '''
-1161
-1162 if self.N == 1:
-1163 raise Exception('Method cannot be applied to one-dimensional correlators.')
-1164 if basematrix is None:
-1165 basematrix = self
-1166 if Ntrunc >= basematrix.N:
-1167 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
-1168 if basematrix.N != self.N:
-1169 raise Exception('basematrix and targetmatrix have to be of the same size.')
-1170
-1171 evecs = []
-1172 for i in range(Ntrunc):
-1173 evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
+1134 This method can be used to reduce the size of an (N x N) correlation matrix
+1135 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1136 is still small.
+1137
+1138 Parameters
+1139 ----------
+1140 Ntrunc: int
+1141 Rank of the target matrix.
+1142 tproj: int
+1143 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1144 The default value is 3.
+1145 t0proj: int
+1146 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1147 discouraged for O(a) improved theories, since the correctness of the procedure
+1148 cannot be granted in this case. The default value is 2.
+1149 basematrix : Corr
+1150 Correlation matrix that is used to determine the eigenvectors of the
+1151 lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1152 is is not specified.
+1153
+1154 Notes
+1155 -----
+1156 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1157 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
+1158 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1159 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1160 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1161 correlation matrix and to remove some noise that is added by irrelevant operators.
+1162 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1163 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1164 '''
+1165
+1166 if self.N == 1:
+1167 raise Exception('Method cannot be applied to one-dimensional correlators.')
+1168 if basematrix is None:
+1169 basematrix = self
+1170 if Ntrunc >= basematrix.N:
+1171 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
+1172 if basematrix.N != self.N:
+1173 raise Exception('basematrix and targetmatrix have to be of the same size.')
1174
-1175 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1176 rmat = []
-1177 for t in range(basematrix.T):
-1178 for i in range(Ntrunc):
-1179 for j in range(Ntrunc):
-1180 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1181 rmat.append(np.copy(tmpmat))
-1182
-1183 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
-1184 return Corr(newcontent)
+1175 evecs = []
+1176 for i in range(Ntrunc):
+1177 evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None))
+1178
+1179 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1180 rmat = []
+1181 for t in range(basematrix.T):
+1182 for i in range(Ntrunc):
+1183 for j in range(Ntrunc):
+1184 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1185 rmat.append(np.copy(tmpmat))
+1186
+1187 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
+1188 return Corr(newcontent)