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)
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)
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)))
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])
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)
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)
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)
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
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.")
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.")
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.')
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
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)
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
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.")
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()
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))
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))
View Source
- 1040 def sqrt(self):
-1041 return self**0.5
+ 1042 def sqrt(self):
+1043 return self**0.5
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)