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