pyerrors.correlators
1import warnings 2from itertools import permutations 3import numpy as np 4import autograd.numpy as anp 5import matplotlib.pyplot as plt 6import scipy.linalg 7from .obs import Obs, reweight, correlate, CObs 8from .misc import dump_object, _assert_equal_properties 9from .fits import least_squares 10from .roots import find_root 11 12 13class Corr: 14 """The class for a correlator (time dependent sequence of pe.Obs). 15 16 Everything, this class does, can be achieved using lists or arrays of Obs. 17 But it is simply more convenient to have a dedicated object for correlators. 18 One often wants to add or multiply correlators of the same length at every timeslice and it is inconvenient 19 to iterate over all timeslices for every operation. This is especially true, when dealing with matrices. 20 21 The correlator can have two types of content: An Obs at every timeslice OR a GEVP 22 matrix at every timeslice. Other dependency (eg. spatial) are not supported. 23 24 """ 25 26 def __init__(self, data_input, padding=[0, 0], prange=None): 27 """ Initialize a Corr object. 28 29 Parameters 30 ---------- 31 data_input : list or array 32 list of Obs or list of arrays of Obs or array of Corrs 33 padding : list, optional 34 List with two entries where the first labels the padding 35 at the front of the correlator and the second the padding 36 at the back. 37 prange : list, optional 38 List containing the first and last timeslice of the plateau 39 region indentified for this correlator. 40 """ 41 42 if isinstance(data_input, np.ndarray): 43 44 # This only works, if the array fulfills the conditions below 45 if not len(data_input.shape) == 2 and data_input.shape[0] == data_input.shape[1]: 46 raise Exception("Incompatible array shape") 47 if not all([isinstance(item, Corr) for item in data_input.flatten()]): 48 raise Exception("If the input is an array, its elements must be of type pe.Corr") 49 if not all([item.N == 1 for item in data_input.flatten()]): 50 raise Exception("Can only construct matrix correlator from single valued correlators") 51 if not len(set([item.T for item in data_input.flatten()])) == 1: 52 raise Exception("All input Correlators must be defined over the same timeslices.") 53 54 T = data_input[0, 0].T 55 N = data_input.shape[0] 56 input_as_list = [] 57 for t in range(T): 58 if any([(item.content[t] is None) for item in data_input.flatten()]): 59 if not all([(item.content[t] is None) for item in data_input.flatten()]): 60 warnings.warn("Input ill-defined at different timeslices. Conversion leads to data loss!", RuntimeWarning) 61 input_as_list.append(None) 62 else: 63 array_at_timeslace = np.empty([N, N], dtype="object") 64 for i in range(N): 65 for j in range(N): 66 array_at_timeslace[i, j] = data_input[i, j][t] 67 input_as_list.append(array_at_timeslace) 68 data_input = input_as_list 69 70 if isinstance(data_input, list): 71 72 if all([isinstance(item, (Obs, CObs)) for item in data_input]): 73 _assert_equal_properties(data_input) 74 self.content = [np.asarray([item]) for item in data_input] 75 if all([isinstance(item, (Obs, CObs)) or item is None for item in data_input]): 76 _assert_equal_properties([o for o in data_input if o is not None]) 77 self.content = [np.asarray([item]) if item is not None else None for item in data_input] 78 self.N = 1 79 80 elif all([isinstance(item, np.ndarray) or item is None for item in data_input]) and any([isinstance(item, np.ndarray) for item in data_input]): 81 self.content = data_input 82 noNull = [a for a in self.content if not (a is None)] # To check if the matrices are correct for all undefined elements 83 self.N = noNull[0].shape[0] 84 if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]: 85 raise Exception("Smearing matrices are not NxN") 86 if (not all([item.shape == noNull[0].shape for item in noNull])): 87 raise Exception("Items in data_input are not of identical shape." + str(noNull)) 88 else: 89 raise Exception("data_input contains item of wrong type") 90 else: 91 raise Exception("Data input was not given as list or correct array") 92 93 self.tag = None 94 95 # An undefined timeslice is represented by the None object 96 self.content = [None] * padding[0] + self.content + [None] * padding[1] 97 self.T = len(self.content) 98 self.prange = prange 99 100 def __getitem__(self, idx): 101 """Return the content of timeslice idx""" 102 if self.content[idx] is None: 103 return None 104 elif len(self.content[idx]) == 1: 105 return self.content[idx][0] 106 else: 107 return self.content[idx] 108 109 @property 110 def reweighted(self): 111 bool_array = np.array([list(map(lambda x: x.reweighted, o)) for o in [x for x in self.content if x is not None]]) 112 if np.all(bool_array == 1): 113 return True 114 elif np.all(bool_array == 0): 115 return False 116 else: 117 raise Exception("Reweighting status of correlator corrupted.") 118 119 def gamma_method(self, **kwargs): 120 """Apply the gamma method to the content of the Corr.""" 121 for item in self.content: 122 if not(item is None): 123 if self.N == 1: 124 item[0].gamma_method(**kwargs) 125 else: 126 for i in range(self.N): 127 for j in range(self.N): 128 item[i, j].gamma_method(**kwargs) 129 130 def projected(self, vector_l=None, vector_r=None, normalize=False): 131 """We need to project the Correlator with a Vector to get a single value at each timeslice. 132 133 The method can use one or two vectors. 134 If two are specified it returns v1@G@v2 (the order might be very important.) 135 By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to 136 """ 137 if self.N == 1: 138 raise Exception("Trying to project a Corr, that already has N=1.") 139 140 if vector_l is None: 141 vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.]) 142 elif(vector_r is None): 143 vector_r = vector_l 144 if isinstance(vector_l, list) and not isinstance(vector_r, list): 145 if len(vector_l) != self.T: 146 raise Exception("Length of vector list must be equal to T") 147 vector_r = [vector_r] * self.T 148 if isinstance(vector_r, list) and not isinstance(vector_l, list): 149 if len(vector_r) != self.T: 150 raise Exception("Length of vector list must be equal to T") 151 vector_l = [vector_l] * self.T 152 153 if not isinstance(vector_l, list): 154 if not vector_l.shape == vector_r.shape == (self.N,): 155 raise Exception("Vectors are of wrong shape!") 156 if normalize: 157 vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r) 158 newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content] 159 160 else: 161 # There are no checks here yet. There are so many possible scenarios, where this can go wrong. 162 if normalize: 163 for t in range(self.T): 164 vector_l[t], vector_r[t] = vector_l[t] / np.sqrt((vector_l[t] @ vector_l[t])), vector_r[t] / np.sqrt(vector_r[t] @ vector_r[t]) 165 166 newcontent = [None if (self.content[t] is None or vector_l[t] is None or vector_r[t] is None) else np.asarray([vector_l[t].T @ self.content[t] @ vector_r[t]]) for t in range(self.T)] 167 return Corr(newcontent) 168 169 def item(self, i, j): 170 """Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice. 171 172 Parameters 173 ---------- 174 i : int 175 First index to be picked. 176 j : int 177 Second index to be picked. 178 """ 179 if self.N == 1: 180 raise Exception("Trying to pick item from projected Corr") 181 newcontent = [None if(item is None) else item[i, j] for item in self.content] 182 return Corr(newcontent) 183 184 def plottable(self): 185 """Outputs the correlator in a plotable format. 186 187 Outputs three lists containing the timeslice index, the value on each 188 timeslice and the error on each timeslice. 189 """ 190 if self.N != 1: 191 raise Exception("Can only make Corr[N=1] plottable") 192 x_list = [x for x in range(self.T) if not self.content[x] is None] 193 y_list = [y[0].value for y in self.content if y is not None] 194 y_err_list = [y[0].dvalue for y in self.content if y is not None] 195 196 return x_list, y_list, y_err_list 197 198 def symmetric(self): 199 """ Symmetrize the correlator around x0=0.""" 200 if self.T % 2 != 0: 201 raise Exception("Can not symmetrize odd T") 202 203 if np.argmax(np.abs(self.content)) != 0: 204 warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning) 205 206 newcontent = [self.content[0]] 207 for t in range(1, self.T): 208 if (self.content[t] is None) or (self.content[self.T - t] is None): 209 newcontent.append(None) 210 else: 211 newcontent.append(0.5 * (self.content[t] + self.content[self.T - t])) 212 if(all([x is None for x in newcontent])): 213 raise Exception("Corr could not be symmetrized: No redundant values") 214 return Corr(newcontent, prange=self.prange) 215 216 def anti_symmetric(self): 217 """Anti-symmetrize the correlator around x0=0.""" 218 if self.T % 2 != 0: 219 raise Exception("Can not symmetrize odd T") 220 221 test = 1 * self 222 test.gamma_method() 223 if not all([o.is_zero_within_error(3) for o in test.content[0]]): 224 warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning) 225 226 newcontent = [self.content[0]] 227 for t in range(1, self.T): 228 if (self.content[t] is None) or (self.content[self.T - t] is None): 229 newcontent.append(None) 230 else: 231 newcontent.append(0.5 * (self.content[t] - self.content[self.T - t])) 232 if(all([x is None for x in newcontent])): 233 raise Exception("Corr could not be symmetrized: No redundant values") 234 return Corr(newcontent, prange=self.prange) 235 236 def matrix_symmetric(self): 237 """Symmetrizes the correlator matrices on every timeslice.""" 238 if self.N > 1: 239 transposed = [None if (G is None) else G.T for G in self.content] 240 return 0.5 * (Corr(transposed) + self) 241 if self.N == 1: 242 raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.") 243 244 def GEVP(self, t0, ts=None, state=0, sorted_list="Eigenvalue"): 245 """Solve the general eigenvalue problem on the current correlator 246 247 Parameters 248 ---------- 249 t0 : int 250 The time t0 for G(t)v= lambda G(t_0)v 251 ts : int 252 fixed time G(t_s)v= lambda G(t_0)v if return_list=False 253 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. 254 state : int 255 The state one is interested in, ordered by energy. The lowest state is zero. 256 sorted_list : string 257 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. 258 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. 259 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. 260 The reference state is identified by its eigenvalue at t=ts 261 """ 262 263 if self.N == 1: 264 raise Exception("GEVP methods only works on correlator matrices and not single correlators.") 265 266 symmetric_corr = self.matrix_symmetric() 267 if sorted_list is None: 268 if (ts is None): 269 raise Exception("ts is required if sorted_list=None.") 270 if (ts <= t0): 271 raise Exception("ts has to be larger than t0.") 272 if (self.content[t0] is None) or (self.content[ts] is None): 273 raise Exception("Corr not defined at t0/ts.") 274 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") 275 for i in range(self.N): 276 for j in range(self.N): 277 G0[i, j] = symmetric_corr[t0][i, j].value 278 Gt[i, j] = symmetric_corr[ts][i, j].value 279 280 sp_vecs = _GEVP_solver(Gt, G0) 281 sp_vec = sp_vecs[state] 282 return sp_vec 283 elif sorted_list in ["Eigenvalue", "Eigenvector"]: 284 if sorted_list == "Eigenvalue" and ts is not None: 285 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) 286 all_vecs = [None] * (t0 + 1) 287 for t in range(t0 + 1, self.T): 288 try: 289 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") 290 for i in range(self.N): 291 for j in range(self.N): 292 G0[i, j] = symmetric_corr[t0][i, j].value 293 Gt[i, j] = symmetric_corr[t][i, j].value 294 295 sp_vecs = _GEVP_solver(Gt, G0) 296 if sorted_list == "Eigenvalue": 297 sp_vec = sp_vecs[state] 298 all_vecs.append(sp_vec) 299 else: 300 all_vecs.append(sp_vecs) 301 except Exception: 302 all_vecs.append(None) 303 if sorted_list == "Eigenvector": 304 if (ts is None): 305 raise Exception("ts is required for the Eigenvector sorting method.") 306 all_vecs = _sort_vectors(all_vecs, ts) 307 all_vecs = [a[state] if a is not None else None for a in all_vecs] 308 else: 309 raise Exception("Unkown value for 'sorted_list'.") 310 311 return all_vecs 312 313 def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None): 314 """Determines the eigenvalue of the GEVP by solving and projecting the correlator 315 316 Parameters 317 ---------- 318 t0 : int 319 The time t0 for G(t)v= lambda G(t_0)v 320 ts : int 321 fixed time G(t_s)v= lambda G(t_0)v if return_list=False 322 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. 323 state : int 324 The state one is interested in ordered by energy. The lowest state is zero. 325 sorted_list : string 326 if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. 327 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. 328 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. 329 The reference state is identified by its eigenvalue at t=ts 330 """ 331 vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list) 332 return self.projected(vec) 333 334 def Hankel(self, N, periodic=False): 335 """Constructs an NxN Hankel matrix 336 337 C(t) c(t+1) ... c(t+n-1) 338 C(t+1) c(t+2) ... c(t+n) 339 ................. 340 C(t+(n-1)) c(t+n) ... c(t+2(n-1)) 341 342 Parameters 343 ---------- 344 N : int 345 Dimension of the Hankel matrix 346 periodic : bool, optional 347 determines whether the matrix is extended periodically 348 """ 349 350 if self.N != 1: 351 raise Exception("Multi-operator Prony not implemented!") 352 353 array = np.empty([N, N], dtype="object") 354 new_content = [] 355 for t in range(self.T): 356 new_content.append(array.copy()) 357 358 def wrap(i): 359 while i >= self.T: 360 i -= self.T 361 return i 362 363 for t in range(self.T): 364 for i in range(N): 365 for j in range(N): 366 if periodic: 367 new_content[t][i, j] = self.content[wrap(t + i + j)][0] 368 elif (t + i + j) >= self.T: 369 new_content[t] = None 370 else: 371 new_content[t][i, j] = self.content[t + i + j][0] 372 373 return Corr(new_content) 374 375 def roll(self, dt): 376 """Periodically shift the correlator by dt timeslices 377 378 Parameters 379 ---------- 380 dt : int 381 number of timeslices 382 """ 383 return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) 384 385 def reverse(self): 386 """Reverse the time ordering of the Corr""" 387 return Corr(self.content[:: -1]) 388 389 def thin(self, spacing=2, offset=0): 390 """Thin out a correlator to suppress correlations 391 392 Parameters 393 ---------- 394 spacing : int 395 Keep only every 'spacing'th entry of the correlator 396 offset : int 397 Offset the equal spacing 398 """ 399 new_content = [] 400 for t in range(self.T): 401 if (offset + t) % spacing != 0: 402 new_content.append(None) 403 else: 404 new_content.append(self.content[t]) 405 return Corr(new_content) 406 407 def correlate(self, partner): 408 """Correlate the correlator with another correlator or Obs 409 410 Parameters 411 ---------- 412 partner : Obs or Corr 413 partner to correlate the correlator with. 414 Can either be an Obs which is correlated with all entries of the 415 correlator or a Corr of same length. 416 """ 417 new_content = [] 418 for x0, t_slice in enumerate(self.content): 419 if t_slice is None: 420 new_content.append(None) 421 else: 422 if isinstance(partner, Corr): 423 if partner.content[x0] is None: 424 new_content.append(None) 425 else: 426 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) 427 elif isinstance(partner, Obs): # Should this include CObs? 428 new_content.append(np.array([correlate(o, partner) for o in t_slice])) 429 else: 430 raise Exception("Can only correlate with an Obs or a Corr.") 431 432 return Corr(new_content) 433 434 def reweight(self, weight, **kwargs): 435 """Reweight the correlator. 436 437 Parameters 438 ---------- 439 weight : Obs 440 Reweighting factor. An Observable that has to be defined on a superset of the 441 configurations in obs[i].idl for all i. 442 all_configs : bool 443 if True, the reweighted observables are normalized by the average of 444 the reweighting factor on all configurations in weight.idl and not 445 on the configurations in obs[i].idl. 446 """ 447 new_content = [] 448 for t_slice in self.content: 449 if t_slice is None: 450 new_content.append(None) 451 else: 452 new_content.append(np.array(reweight(weight, t_slice, **kwargs))) 453 return Corr(new_content) 454 455 def T_symmetry(self, partner, parity=+1): 456 """Return the time symmetry average of the correlator and its partner 457 458 Parameters 459 ---------- 460 partner : Corr 461 Time symmetry partner of the Corr 462 partity : int 463 Parity quantum number of the correlator, can be +1 or -1 464 """ 465 if not isinstance(partner, Corr): 466 raise Exception("T partner has to be a Corr object.") 467 if parity not in [+1, -1]: 468 raise Exception("Parity has to be +1 or -1.") 469 T_partner = parity * partner.reverse() 470 471 t_slices = [] 472 test = (self - T_partner) 473 test.gamma_method() 474 for x0, t_slice in enumerate(test.content): 475 if t_slice is not None: 476 if not t_slice[0].is_zero_within_error(5): 477 t_slices.append(x0) 478 if t_slices: 479 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning) 480 481 return (self + T_partner) / 2 482 483 def deriv(self, variant="symmetric"): 484 """Return the first derivative of the correlator with respect to x0. 485 486 Parameters 487 ---------- 488 variant : str 489 decides which definition of the finite differences derivative is used. 490 Available choice: symmetric, forward, backward, improved, default: symmetric 491 """ 492 if variant == "symmetric": 493 newcontent = [] 494 for t in range(1, self.T - 1): 495 if (self.content[t - 1] is None) or (self.content[t + 1] is None): 496 newcontent.append(None) 497 else: 498 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) 499 if(all([x is None for x in newcontent])): 500 raise Exception('Derivative is undefined at all timeslices') 501 return Corr(newcontent, padding=[1, 1]) 502 elif variant == "forward": 503 newcontent = [] 504 for t in range(self.T - 1): 505 if (self.content[t] is None) or (self.content[t + 1] is None): 506 newcontent.append(None) 507 else: 508 newcontent.append(self.content[t + 1] - self.content[t]) 509 if(all([x is None for x in newcontent])): 510 raise Exception("Derivative is undefined at all timeslices") 511 return Corr(newcontent, padding=[0, 1]) 512 elif variant == "backward": 513 newcontent = [] 514 for t in range(1, self.T): 515 if (self.content[t - 1] is None) or (self.content[t] is None): 516 newcontent.append(None) 517 else: 518 newcontent.append(self.content[t] - self.content[t - 1]) 519 if(all([x is None for x in newcontent])): 520 raise Exception("Derivative is undefined at all timeslices") 521 return Corr(newcontent, padding=[1, 0]) 522 elif variant == "improved": 523 newcontent = [] 524 for t in range(2, self.T - 2): 525 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None): 526 newcontent.append(None) 527 else: 528 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2])) 529 if(all([x is None for x in newcontent])): 530 raise Exception('Derivative is undefined at all timeslices') 531 return Corr(newcontent, padding=[2, 2]) 532 else: 533 raise Exception("Unknown variant.") 534 535 def second_deriv(self, variant="symmetric"): 536 """Return the second derivative of the correlator with respect to x0. 537 538 Parameters 539 ---------- 540 variant : str 541 decides which definition of the finite differences derivative is used. 542 Available choice: symmetric, improved, default: symmetric 543 """ 544 if variant == "symmetric": 545 newcontent = [] 546 for t in range(1, self.T - 1): 547 if (self.content[t - 1] is None) or (self.content[t + 1] is None): 548 newcontent.append(None) 549 else: 550 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) 551 if(all([x is None for x in newcontent])): 552 raise Exception("Derivative is undefined at all timeslices") 553 return Corr(newcontent, padding=[1, 1]) 554 elif variant == "improved": 555 newcontent = [] 556 for t in range(2, self.T - 2): 557 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None): 558 newcontent.append(None) 559 else: 560 newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2])) 561 if(all([x is None for x in newcontent])): 562 raise Exception("Derivative is undefined at all timeslices") 563 return Corr(newcontent, padding=[2, 2]) 564 else: 565 raise Exception("Unknown variant.") 566 567 def m_eff(self, variant='log', guess=1.0): 568 """Returns the effective mass of the correlator as correlator object 569 570 Parameters 571 ---------- 572 variant : str 573 log : uses the standard effective mass log(C(t) / C(t+1)) 574 cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. 575 sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m. 576 See, e.g., arXiv:1205.5380 577 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) 578 guess : float 579 guess for the root finder, only relevant for the root variant 580 """ 581 if self.N != 1: 582 raise Exception('Correlator must be projected before getting m_eff') 583 if variant == 'log': 584 newcontent = [] 585 for t in range(self.T - 1): 586 if (self.content[t] is None) or (self.content[t + 1] is None): 587 newcontent.append(None) 588 else: 589 newcontent.append(self.content[t] / self.content[t + 1]) 590 if(all([x is None for x in newcontent])): 591 raise Exception('m_eff is undefined at all timeslices') 592 593 return np.log(Corr(newcontent, padding=[0, 1])) 594 595 elif variant in ['periodic', 'cosh', 'sinh']: 596 if variant in ['periodic', 'cosh']: 597 func = anp.cosh 598 else: 599 func = anp.sinh 600 601 def root_function(x, d): 602 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d 603 604 newcontent = [] 605 for t in range(self.T - 1): 606 if (self.content[t] is None) or (self.content[t + 1] is None): 607 newcontent.append(None) 608 # Fill the two timeslices in the middle of the lattice with their predecessors 609 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: 610 newcontent.append(newcontent[-1]) 611 else: 612 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) 613 if(all([x is None for x in newcontent])): 614 raise Exception('m_eff is undefined at all timeslices') 615 616 return Corr(newcontent, padding=[0, 1]) 617 618 elif variant == 'arccosh': 619 newcontent = [] 620 for t in range(1, self.T - 1): 621 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None): 622 newcontent.append(None) 623 else: 624 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) 625 if(all([x is None for x in newcontent])): 626 raise Exception("m_eff is undefined at all timeslices") 627 return np.arccosh(Corr(newcontent, padding=[1, 1])) 628 629 else: 630 raise Exception('Unknown variant.') 631 632 def fit(self, function, fitrange=None, silent=False, **kwargs): 633 r'''Fits function to the data 634 635 Parameters 636 ---------- 637 function : obj 638 function to fit to the data. See fits.least_squares for details. 639 fitrange : list 640 Two element list containing the timeslices on which the fit is supposed to start and stop. 641 Caution: This range is inclusive as opposed to standard python indexing. 642 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. 643 If not specified, self.prange or all timeslices are used. 644 silent : bool 645 Decides whether output is printed to the standard output. 646 ''' 647 if self.N != 1: 648 raise Exception("Correlator must be projected before fitting") 649 650 if fitrange is None: 651 if self.prange: 652 fitrange = self.prange 653 else: 654 fitrange = [0, self.T - 1] 655 else: 656 if not isinstance(fitrange, list): 657 raise Exception("fitrange has to be a list with two elements") 658 if len(fitrange) != 2: 659 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") 660 661 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] 662 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] 663 result = least_squares(xs, ys, function, silent=silent, **kwargs) 664 return result 665 666 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): 667 """ Extract a plateau value from a Corr object 668 669 Parameters 670 ---------- 671 plateau_range : list 672 list with two entries, indicating the first and the last timeslice 673 of the plateau region. 674 method : str 675 method to extract the plateau. 676 'fit' fits a constant to the plateau region 677 'avg', 'average' or 'mean' just average over the given timeslices. 678 auto_gamma : bool 679 apply gamma_method with default parameters to the Corr. Defaults to None 680 """ 681 if not plateau_range: 682 if self.prange: 683 plateau_range = self.prange 684 else: 685 raise Exception("no plateau range provided") 686 if self.N != 1: 687 raise Exception("Correlator must be projected before getting a plateau.") 688 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): 689 raise Exception("plateau is undefined at all timeslices in plateaurange.") 690 if auto_gamma: 691 self.gamma_method() 692 if method == "fit": 693 def const_func(a, t): 694 return a[0] 695 return self.fit(const_func, plateau_range)[0] 696 elif method in ["avg", "average", "mean"]: 697 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) 698 return returnvalue 699 700 else: 701 raise Exception("Unsupported plateau method: " + method) 702 703 def set_prange(self, prange): 704 """Sets the attribute prange of the Corr object.""" 705 if not len(prange) == 2: 706 raise Exception("prange must be a list or array with two values") 707 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): 708 raise Exception("Start and end point must be integers") 709 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): 710 raise Exception("Start and end point must define a range in the interval 0,T") 711 712 self.prange = prange 713 return 714 715 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None): 716 """Plots the correlator using the tag of the correlator as label if available. 717 718 Parameters 719 ---------- 720 x_range : list 721 list of two values, determining the range of the x-axis e.g. [4, 8] 722 comp : Corr or list of Corr 723 Correlator or list of correlators which are plotted for comparison. 724 The tags of these correlators are used as labels if available. 725 logscale : bool 726 Sets y-axis to logscale 727 plateau : Obs 728 Plateau value to be visualized in the figure 729 fit_res : Fit_result 730 Fit_result object to be visualized 731 ylabel : str 732 Label for the y-axis 733 save : str 734 path to file in which the figure should be saved 735 auto_gamma : bool 736 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. 737 hide_sigma : float 738 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. 739 references : list 740 List of floating point values that are displayed as horizontal lines for reference. 741 """ 742 if self.N != 1: 743 raise Exception("Correlator must be projected before plotting") 744 745 if auto_gamma: 746 self.gamma_method() 747 748 if x_range is None: 749 x_range = [0, self.T - 1] 750 751 fig = plt.figure() 752 ax1 = fig.add_subplot(111) 753 754 x, y, y_err = self.plottable() 755 if hide_sigma: 756 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1 757 else: 758 hide_from = None 759 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) 760 if logscale: 761 ax1.set_yscale('log') 762 else: 763 if y_range is None: 764 try: 765 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)]) 766 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)]) 767 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 768 except Exception: 769 pass 770 else: 771 ax1.set_ylim(y_range) 772 if comp: 773 if isinstance(comp, (Corr, list)): 774 for corr in comp if isinstance(comp, list) else [comp]: 775 if auto_gamma: 776 corr.gamma_method() 777 x, y, y_err = corr.plottable() 778 if hide_sigma: 779 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1 780 else: 781 hide_from = None 782 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 783 else: 784 raise Exception("'comp' must be a correlator or a list of correlators.") 785 786 if plateau: 787 if isinstance(plateau, Obs): 788 if auto_gamma: 789 plateau.gamma_method() 790 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 791 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 792 else: 793 raise Exception("'plateau' must be an Obs") 794 795 if references: 796 if isinstance(references, list): 797 for ref in references: 798 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') 799 else: 800 raise Exception("'references' must be a list of floating pint values.") 801 802 if self.prange: 803 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 804 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 805 806 if fit_res: 807 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 808 ax1.plot(x_samples, 809 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 810 ls='-', marker=',', lw=2) 811 812 ax1.set_xlabel(r'$x_0 / a$') 813 if ylabel: 814 ax1.set_ylabel(ylabel) 815 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 816 817 handles, labels = ax1.get_legend_handles_labels() 818 if labels: 819 ax1.legend() 820 plt.draw() 821 822 if save: 823 if isinstance(save, str): 824 fig.savefig(save) 825 else: 826 raise Exception("'save' has to be a string.") 827 828 def spaghetti_plot(self, logscale=True): 829 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 830 831 Parameters 832 ---------- 833 logscale : bool 834 Determines whether the scale of the y-axis is logarithmic or standard. 835 """ 836 if self.N != 1: 837 raise Exception("Correlator needs to be projected first.") 838 839 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist])) 840 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 841 842 for name in mc_names: 843 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 844 845 fig = plt.figure() 846 ax = fig.add_subplot(111) 847 for dat in data: 848 ax.plot(x0_vals, dat, ls='-', marker='') 849 850 if logscale is True: 851 ax.set_yscale('log') 852 853 ax.set_xlabel(r'$x_0 / a$') 854 plt.title(name) 855 plt.draw() 856 857 def dump(self, filename, datatype="json.gz", **kwargs): 858 """Dumps the Corr into a file of chosen type 859 Parameters 860 ---------- 861 filename : str 862 Name of the file to be saved. 863 datatype : str 864 Format of the exported file. Supported formats include 865 "json.gz" and "pickle" 866 path : str 867 specifies a custom path for the file (default '.') 868 """ 869 if datatype == "json.gz": 870 from .input.json import dump_to_json 871 if 'path' in kwargs: 872 file_name = kwargs.get('path') + '/' + filename 873 else: 874 file_name = filename 875 dump_to_json(self, file_name) 876 elif datatype == "pickle": 877 dump_object(self, filename, **kwargs) 878 else: 879 raise Exception("Unknown datatype " + str(datatype)) 880 881 def print(self, range=[0, None]): 882 print(self.__repr__(range)) 883 884 def __repr__(self, range=[0, None]): 885 content_string = "" 886 887 content_string += "Corr T=" + str(self.T) + " N=" + str(self.N) + "\n" # +" filled with"+ str(type(self.content[0][0])) there should be a good solution here 888 889 if self.tag is not None: 890 content_string += "Description: " + self.tag + "\n" 891 if self.N != 1: 892 return content_string 893 894 if range[1]: 895 range[1] += 1 896 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' 897 for i, sub_corr in enumerate(self.content[range[0]:range[1]]): 898 if sub_corr is None: 899 content_string += str(i + range[0]) + '\n' 900 else: 901 content_string += str(i + range[0]) 902 for element in sub_corr: 903 content_string += '\t' + ' ' * int(element >= 0) + str(element) 904 content_string += '\n' 905 return content_string 906 907 def __str__(self): 908 return self.__repr__() 909 910 # We define the basic operations, that can be performed with correlators. 911 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. 912 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. 913 # One could try and tell Obs to check if the y in __mul__ is a Corr and 914 915 def __add__(self, y): 916 if isinstance(y, Corr): 917 if ((self.N != y.N) or (self.T != y.T)): 918 raise Exception("Addition of Corrs with different shape") 919 newcontent = [] 920 for t in range(self.T): 921 if (self.content[t] is None) or (y.content[t] is None): 922 newcontent.append(None) 923 else: 924 newcontent.append(self.content[t] + y.content[t]) 925 return Corr(newcontent) 926 927 elif isinstance(y, (Obs, int, float, CObs)): 928 newcontent = [] 929 for t in range(self.T): 930 if (self.content[t] is None): 931 newcontent.append(None) 932 else: 933 newcontent.append(self.content[t] + y) 934 return Corr(newcontent, prange=self.prange) 935 elif isinstance(y, np.ndarray): 936 if y.shape == (self.T,): 937 return Corr(list((np.array(self.content).T + y).T)) 938 else: 939 raise ValueError("operands could not be broadcast together") 940 else: 941 raise TypeError("Corr + wrong type") 942 943 def __mul__(self, y): 944 if isinstance(y, Corr): 945 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 946 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 947 newcontent = [] 948 for t in range(self.T): 949 if (self.content[t] is None) or (y.content[t] is None): 950 newcontent.append(None) 951 else: 952 newcontent.append(self.content[t] * y.content[t]) 953 return Corr(newcontent) 954 955 elif isinstance(y, (Obs, int, float, CObs)): 956 newcontent = [] 957 for t in range(self.T): 958 if (self.content[t] is None): 959 newcontent.append(None) 960 else: 961 newcontent.append(self.content[t] * y) 962 return Corr(newcontent, prange=self.prange) 963 elif isinstance(y, np.ndarray): 964 if y.shape == (self.T,): 965 return Corr(list((np.array(self.content).T * y).T)) 966 else: 967 raise ValueError("operands could not be broadcast together") 968 else: 969 raise TypeError("Corr * wrong type") 970 971 def __truediv__(self, y): 972 if isinstance(y, Corr): 973 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 974 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 975 newcontent = [] 976 for t in range(self.T): 977 if (self.content[t] is None) or (y.content[t] is None): 978 newcontent.append(None) 979 else: 980 newcontent.append(self.content[t] / y.content[t]) 981 for t in range(self.T): 982 if newcontent[t] is None: 983 continue 984 if np.isnan(np.sum(newcontent[t]).value): 985 newcontent[t] = None 986 987 if all([item is None for item in newcontent]): 988 raise Exception("Division returns completely undefined correlator") 989 return Corr(newcontent) 990 991 elif isinstance(y, (Obs, CObs)): 992 if isinstance(y, Obs): 993 if y.value == 0: 994 raise Exception('Division by zero will return undefined correlator') 995 if isinstance(y, CObs): 996 if y.is_zero(): 997 raise Exception('Division by zero will return undefined correlator') 998 999 newcontent = [] 1000 for t in range(self.T): 1001 if (self.content[t] is None): 1002 newcontent.append(None) 1003 else: 1004 newcontent.append(self.content[t] / y) 1005 return Corr(newcontent, prange=self.prange) 1006 1007 elif isinstance(y, (int, float)): 1008 if y == 0: 1009 raise Exception('Division by zero will return undefined correlator') 1010 newcontent = [] 1011 for t in range(self.T): 1012 if (self.content[t] is None): 1013 newcontent.append(None) 1014 else: 1015 newcontent.append(self.content[t] / y) 1016 return Corr(newcontent, prange=self.prange) 1017 elif isinstance(y, np.ndarray): 1018 if y.shape == (self.T,): 1019 return Corr(list((np.array(self.content).T / y).T)) 1020 else: 1021 raise ValueError("operands could not be broadcast together") 1022 else: 1023 raise TypeError('Corr / wrong type') 1024 1025 def __neg__(self): 1026 newcontent = [None if (item is None) else -1. * item for item in self.content] 1027 return Corr(newcontent, prange=self.prange) 1028 1029 def __sub__(self, y): 1030 return self + (-y) 1031 1032 def __pow__(self, y): 1033 if isinstance(y, (Obs, int, float, CObs)): 1034 newcontent = [None if (item is None) else item**y for item in self.content] 1035 return Corr(newcontent, prange=self.prange) 1036 else: 1037 raise TypeError('Type of exponent not supported') 1038 1039 def __abs__(self): 1040 newcontent = [None if (item is None) else np.abs(item) for item in self.content] 1041 return Corr(newcontent, prange=self.prange) 1042 1043 # The numpy functions: 1044 def sqrt(self): 1045 return self**0.5 1046 1047 def log(self): 1048 newcontent = [None if (item is None) else np.log(item) for item in self.content] 1049 return Corr(newcontent, prange=self.prange) 1050 1051 def exp(self): 1052 newcontent = [None if (item is None) else np.exp(item) for item in self.content] 1053 return Corr(newcontent, prange=self.prange) 1054 1055 def _apply_func_to_corr(self, func): 1056 newcontent = [None if (item is None) else func(item) for item in self.content] 1057 for t in range(self.T): 1058 if newcontent[t] is None: 1059 continue 1060 if np.isnan(np.sum(newcontent[t]).value): 1061 newcontent[t] = None 1062 if all([item is None for item in newcontent]): 1063 raise Exception('Operation returns undefined correlator') 1064 return Corr(newcontent) 1065 1066 def sin(self): 1067 return self._apply_func_to_corr(np.sin) 1068 1069 def cos(self): 1070 return self._apply_func_to_corr(np.cos) 1071 1072 def tan(self): 1073 return self._apply_func_to_corr(np.tan) 1074 1075 def sinh(self): 1076 return self._apply_func_to_corr(np.sinh) 1077 1078 def cosh(self): 1079 return self._apply_func_to_corr(np.cosh) 1080 1081 def tanh(self): 1082 return self._apply_func_to_corr(np.tanh) 1083 1084 def arcsin(self): 1085 return self._apply_func_to_corr(np.arcsin) 1086 1087 def arccos(self): 1088 return self._apply_func_to_corr(np.arccos) 1089 1090 def arctan(self): 1091 return self._apply_func_to_corr(np.arctan) 1092 1093 def arcsinh(self): 1094 return self._apply_func_to_corr(np.arcsinh) 1095 1096 def arccosh(self): 1097 return self._apply_func_to_corr(np.arccosh) 1098 1099 def arctanh(self): 1100 return self._apply_func_to_corr(np.arctanh) 1101 1102 # Right hand side operations (require tweak in main module to work) 1103 def __radd__(self, y): 1104 return self + y 1105 1106 def __rsub__(self, y): 1107 return -self + y 1108 1109 def __rmul__(self, y): 1110 return self * y 1111 1112 def __rtruediv__(self, y): 1113 return (self / y) ** (-1) 1114 1115 @property 1116 def real(self): 1117 def return_real(obs_OR_cobs): 1118 if isinstance(obs_OR_cobs, CObs): 1119 return obs_OR_cobs.real 1120 else: 1121 return obs_OR_cobs 1122 1123 return self._apply_func_to_corr(return_real) 1124 1125 @property 1126 def imag(self): 1127 def return_imag(obs_OR_cobs): 1128 if isinstance(obs_OR_cobs, CObs): 1129 return obs_OR_cobs.imag 1130 else: 1131 return obs_OR_cobs * 0 # So it stays the right type 1132 1133 return self._apply_func_to_corr(return_imag) 1134 1135 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): 1136 r''' Project large correlation matrix to lowest states 1137 1138 This method can be used to reduce the size of an (N x N) correlation matrix 1139 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise 1140 is still small. 1141 1142 Parameters 1143 ---------- 1144 Ntrunc: int 1145 Rank of the target matrix. 1146 tproj: int 1147 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. 1148 The default value is 3. 1149 t0proj: int 1150 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly 1151 discouraged for O(a) improved theories, since the correctness of the procedure 1152 cannot be granted in this case. The default value is 2. 1153 basematrix : Corr 1154 Correlation matrix that is used to determine the eigenvectors of the 1155 lowest states based on a GEVP. basematrix is taken to be the Corr itself if 1156 is is not specified. 1157 1158 Notes 1159 ----- 1160 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving 1161 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$ 1162 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the 1163 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via 1164 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large 1165 correlation matrix and to remove some noise that is added by irrelevant operators. 1166 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated 1167 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. 1168 ''' 1169 1170 if self.N == 1: 1171 raise Exception('Method cannot be applied to one-dimensional correlators.') 1172 if basematrix is None: 1173 basematrix = self 1174 if Ntrunc >= basematrix.N: 1175 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) 1176 if basematrix.N != self.N: 1177 raise Exception('basematrix and targetmatrix have to be of the same size.') 1178 1179 evecs = [] 1180 for i in range(Ntrunc): 1181 evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None)) 1182 1183 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) 1184 rmat = [] 1185 for t in range(basematrix.T): 1186 for i in range(Ntrunc): 1187 for j in range(Ntrunc): 1188 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] 1189 rmat.append(np.copy(tmpmat)) 1190 1191 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] 1192 return Corr(newcontent) 1193 1194 1195def _sort_vectors(vec_set, ts): 1196 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" 1197 reference_sorting = np.array(vec_set[ts]) 1198 N = reference_sorting.shape[0] 1199 sorted_vec_set = [] 1200 for t in range(len(vec_set)): 1201 if vec_set[t] is None: 1202 sorted_vec_set.append(None) 1203 elif not t == ts: 1204 perms = [list(o) for o in permutations([i for i in range(N)], N)] 1205 best_score = 0 1206 for perm in perms: 1207 current_score = 1 1208 for k in range(N): 1209 new_sorting = reference_sorting.copy() 1210 new_sorting[perm[k], :] = vec_set[t][k] 1211 current_score *= abs(np.linalg.det(new_sorting)) 1212 if current_score > best_score: 1213 best_score = current_score 1214 best_perm = perm 1215 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) 1216 else: 1217 sorted_vec_set.append(vec_set[t]) 1218 1219 return sorted_vec_set 1220 1221 1222def _GEVP_solver(Gt, G0): # Just so normalization an sorting does not need to be repeated. Here we could later put in some checks 1223 sp_val, sp_vecs = scipy.linalg.eigh(Gt, G0) 1224 sp_vecs = [sp_vecs[:, np.argsort(sp_val)[-i]] for i in range(1, sp_vecs.shape[0] + 1)] 1225 sp_vecs = [v / np.sqrt((v.T @ G0 @ v)) for v in sp_vecs] 1226 return sp_vecs
14class Corr: 15 """The class for a correlator (time dependent sequence of pe.Obs). 16 17 Everything, this class does, can be achieved using lists or arrays of Obs. 18 But it is simply more convenient to have a dedicated object for correlators. 19 One often wants to add or multiply correlators of the same length at every timeslice and it is inconvenient 20 to iterate over all timeslices for every operation. This is especially true, when dealing with matrices. 21 22 The correlator can have two types of content: An Obs at every timeslice OR a GEVP 23 matrix at every timeslice. Other dependency (eg. spatial) are not supported. 24 25 """ 26 27 def __init__(self, data_input, padding=[0, 0], prange=None): 28 """ Initialize a Corr object. 29 30 Parameters 31 ---------- 32 data_input : list or array 33 list of Obs or list of arrays of Obs or array of Corrs 34 padding : list, optional 35 List with two entries where the first labels the padding 36 at the front of the correlator and the second the padding 37 at the back. 38 prange : list, optional 39 List containing the first and last timeslice of the plateau 40 region indentified for this correlator. 41 """ 42 43 if isinstance(data_input, np.ndarray): 44 45 # This only works, if the array fulfills the conditions below 46 if not len(data_input.shape) == 2 and data_input.shape[0] == data_input.shape[1]: 47 raise Exception("Incompatible array shape") 48 if not all([isinstance(item, Corr) for item in data_input.flatten()]): 49 raise Exception("If the input is an array, its elements must be of type pe.Corr") 50 if not all([item.N == 1 for item in data_input.flatten()]): 51 raise Exception("Can only construct matrix correlator from single valued correlators") 52 if not len(set([item.T for item in data_input.flatten()])) == 1: 53 raise Exception("All input Correlators must be defined over the same timeslices.") 54 55 T = data_input[0, 0].T 56 N = data_input.shape[0] 57 input_as_list = [] 58 for t in range(T): 59 if any([(item.content[t] is None) for item in data_input.flatten()]): 60 if not all([(item.content[t] is None) for item in data_input.flatten()]): 61 warnings.warn("Input ill-defined at different timeslices. Conversion leads to data loss!", RuntimeWarning) 62 input_as_list.append(None) 63 else: 64 array_at_timeslace = np.empty([N, N], dtype="object") 65 for i in range(N): 66 for j in range(N): 67 array_at_timeslace[i, j] = data_input[i, j][t] 68 input_as_list.append(array_at_timeslace) 69 data_input = input_as_list 70 71 if isinstance(data_input, list): 72 73 if all([isinstance(item, (Obs, CObs)) for item in data_input]): 74 _assert_equal_properties(data_input) 75 self.content = [np.asarray([item]) for item in data_input] 76 if all([isinstance(item, (Obs, CObs)) or item is None for item in data_input]): 77 _assert_equal_properties([o for o in data_input if o is not None]) 78 self.content = [np.asarray([item]) if item is not None else None for item in data_input] 79 self.N = 1 80 81 elif all([isinstance(item, np.ndarray) or item is None for item in data_input]) and any([isinstance(item, np.ndarray) for item in data_input]): 82 self.content = data_input 83 noNull = [a for a in self.content if not (a is None)] # To check if the matrices are correct for all undefined elements 84 self.N = noNull[0].shape[0] 85 if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]: 86 raise Exception("Smearing matrices are not NxN") 87 if (not all([item.shape == noNull[0].shape for item in noNull])): 88 raise Exception("Items in data_input are not of identical shape." + str(noNull)) 89 else: 90 raise Exception("data_input contains item of wrong type") 91 else: 92 raise Exception("Data input was not given as list or correct array") 93 94 self.tag = None 95 96 # An undefined timeslice is represented by the None object 97 self.content = [None] * padding[0] + self.content + [None] * padding[1] 98 self.T = len(self.content) 99 self.prange = prange 100 101 def __getitem__(self, idx): 102 """Return the content of timeslice idx""" 103 if self.content[idx] is None: 104 return None 105 elif len(self.content[idx]) == 1: 106 return self.content[idx][0] 107 else: 108 return self.content[idx] 109 110 @property 111 def reweighted(self): 112 bool_array = np.array([list(map(lambda x: x.reweighted, o)) for o in [x for x in self.content if x is not None]]) 113 if np.all(bool_array == 1): 114 return True 115 elif np.all(bool_array == 0): 116 return False 117 else: 118 raise Exception("Reweighting status of correlator corrupted.") 119 120 def gamma_method(self, **kwargs): 121 """Apply the gamma method to the content of the Corr.""" 122 for item in self.content: 123 if not(item is None): 124 if self.N == 1: 125 item[0].gamma_method(**kwargs) 126 else: 127 for i in range(self.N): 128 for j in range(self.N): 129 item[i, j].gamma_method(**kwargs) 130 131 def projected(self, vector_l=None, vector_r=None, normalize=False): 132 """We need to project the Correlator with a Vector to get a single value at each timeslice. 133 134 The method can use one or two vectors. 135 If two are specified it returns v1@G@v2 (the order might be very important.) 136 By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to 137 """ 138 if self.N == 1: 139 raise Exception("Trying to project a Corr, that already has N=1.") 140 141 if vector_l is None: 142 vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.]) 143 elif(vector_r is None): 144 vector_r = vector_l 145 if isinstance(vector_l, list) and not isinstance(vector_r, list): 146 if len(vector_l) != self.T: 147 raise Exception("Length of vector list must be equal to T") 148 vector_r = [vector_r] * self.T 149 if isinstance(vector_r, list) and not isinstance(vector_l, list): 150 if len(vector_r) != self.T: 151 raise Exception("Length of vector list must be equal to T") 152 vector_l = [vector_l] * self.T 153 154 if not isinstance(vector_l, list): 155 if not vector_l.shape == vector_r.shape == (self.N,): 156 raise Exception("Vectors are of wrong shape!") 157 if normalize: 158 vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r) 159 newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content] 160 161 else: 162 # There are no checks here yet. There are so many possible scenarios, where this can go wrong. 163 if normalize: 164 for t in range(self.T): 165 vector_l[t], vector_r[t] = vector_l[t] / np.sqrt((vector_l[t] @ vector_l[t])), vector_r[t] / np.sqrt(vector_r[t] @ vector_r[t]) 166 167 newcontent = [None if (self.content[t] is None or vector_l[t] is None or vector_r[t] is None) else np.asarray([vector_l[t].T @ self.content[t] @ vector_r[t]]) for t in range(self.T)] 168 return Corr(newcontent) 169 170 def item(self, i, j): 171 """Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice. 172 173 Parameters 174 ---------- 175 i : int 176 First index to be picked. 177 j : int 178 Second index to be picked. 179 """ 180 if self.N == 1: 181 raise Exception("Trying to pick item from projected Corr") 182 newcontent = [None if(item is None) else item[i, j] for item in self.content] 183 return Corr(newcontent) 184 185 def plottable(self): 186 """Outputs the correlator in a plotable format. 187 188 Outputs three lists containing the timeslice index, the value on each 189 timeslice and the error on each timeslice. 190 """ 191 if self.N != 1: 192 raise Exception("Can only make Corr[N=1] plottable") 193 x_list = [x for x in range(self.T) if not self.content[x] is None] 194 y_list = [y[0].value for y in self.content if y is not None] 195 y_err_list = [y[0].dvalue for y in self.content if y is not None] 196 197 return x_list, y_list, y_err_list 198 199 def symmetric(self): 200 """ Symmetrize the correlator around x0=0.""" 201 if self.T % 2 != 0: 202 raise Exception("Can not symmetrize odd T") 203 204 if np.argmax(np.abs(self.content)) != 0: 205 warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning) 206 207 newcontent = [self.content[0]] 208 for t in range(1, self.T): 209 if (self.content[t] is None) or (self.content[self.T - t] is None): 210 newcontent.append(None) 211 else: 212 newcontent.append(0.5 * (self.content[t] + self.content[self.T - t])) 213 if(all([x is None for x in newcontent])): 214 raise Exception("Corr could not be symmetrized: No redundant values") 215 return Corr(newcontent, prange=self.prange) 216 217 def anti_symmetric(self): 218 """Anti-symmetrize the correlator around x0=0.""" 219 if self.T % 2 != 0: 220 raise Exception("Can not symmetrize odd T") 221 222 test = 1 * self 223 test.gamma_method() 224 if not all([o.is_zero_within_error(3) for o in test.content[0]]): 225 warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning) 226 227 newcontent = [self.content[0]] 228 for t in range(1, self.T): 229 if (self.content[t] is None) or (self.content[self.T - t] is None): 230 newcontent.append(None) 231 else: 232 newcontent.append(0.5 * (self.content[t] - self.content[self.T - t])) 233 if(all([x is None for x in newcontent])): 234 raise Exception("Corr could not be symmetrized: No redundant values") 235 return Corr(newcontent, prange=self.prange) 236 237 def matrix_symmetric(self): 238 """Symmetrizes the correlator matrices on every timeslice.""" 239 if self.N > 1: 240 transposed = [None if (G is None) else G.T for G in self.content] 241 return 0.5 * (Corr(transposed) + self) 242 if self.N == 1: 243 raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.") 244 245 def GEVP(self, t0, ts=None, state=0, sorted_list="Eigenvalue"): 246 """Solve the general eigenvalue problem on the current correlator 247 248 Parameters 249 ---------- 250 t0 : int 251 The time t0 for G(t)v= lambda G(t_0)v 252 ts : int 253 fixed time G(t_s)v= lambda G(t_0)v if return_list=False 254 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. 255 state : int 256 The state one is interested in, ordered by energy. The lowest state is zero. 257 sorted_list : string 258 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. 259 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. 260 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. 261 The reference state is identified by its eigenvalue at t=ts 262 """ 263 264 if self.N == 1: 265 raise Exception("GEVP methods only works on correlator matrices and not single correlators.") 266 267 symmetric_corr = self.matrix_symmetric() 268 if sorted_list is None: 269 if (ts is None): 270 raise Exception("ts is required if sorted_list=None.") 271 if (ts <= t0): 272 raise Exception("ts has to be larger than t0.") 273 if (self.content[t0] is None) or (self.content[ts] is None): 274 raise Exception("Corr not defined at t0/ts.") 275 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") 276 for i in range(self.N): 277 for j in range(self.N): 278 G0[i, j] = symmetric_corr[t0][i, j].value 279 Gt[i, j] = symmetric_corr[ts][i, j].value 280 281 sp_vecs = _GEVP_solver(Gt, G0) 282 sp_vec = sp_vecs[state] 283 return sp_vec 284 elif sorted_list in ["Eigenvalue", "Eigenvector"]: 285 if sorted_list == "Eigenvalue" and ts is not None: 286 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) 287 all_vecs = [None] * (t0 + 1) 288 for t in range(t0 + 1, self.T): 289 try: 290 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") 291 for i in range(self.N): 292 for j in range(self.N): 293 G0[i, j] = symmetric_corr[t0][i, j].value 294 Gt[i, j] = symmetric_corr[t][i, j].value 295 296 sp_vecs = _GEVP_solver(Gt, G0) 297 if sorted_list == "Eigenvalue": 298 sp_vec = sp_vecs[state] 299 all_vecs.append(sp_vec) 300 else: 301 all_vecs.append(sp_vecs) 302 except Exception: 303 all_vecs.append(None) 304 if sorted_list == "Eigenvector": 305 if (ts is None): 306 raise Exception("ts is required for the Eigenvector sorting method.") 307 all_vecs = _sort_vectors(all_vecs, ts) 308 all_vecs = [a[state] if a is not None else None for a in all_vecs] 309 else: 310 raise Exception("Unkown value for 'sorted_list'.") 311 312 return all_vecs 313 314 def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None): 315 """Determines the eigenvalue of the GEVP by solving and projecting the correlator 316 317 Parameters 318 ---------- 319 t0 : int 320 The time t0 for G(t)v= lambda G(t_0)v 321 ts : int 322 fixed time G(t_s)v= lambda G(t_0)v if return_list=False 323 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. 324 state : int 325 The state one is interested in ordered by energy. The lowest state is zero. 326 sorted_list : string 327 if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. 328 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. 329 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. 330 The reference state is identified by its eigenvalue at t=ts 331 """ 332 vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list) 333 return self.projected(vec) 334 335 def Hankel(self, N, periodic=False): 336 """Constructs an NxN Hankel matrix 337 338 C(t) c(t+1) ... c(t+n-1) 339 C(t+1) c(t+2) ... c(t+n) 340 ................. 341 C(t+(n-1)) c(t+n) ... c(t+2(n-1)) 342 343 Parameters 344 ---------- 345 N : int 346 Dimension of the Hankel matrix 347 periodic : bool, optional 348 determines whether the matrix is extended periodically 349 """ 350 351 if self.N != 1: 352 raise Exception("Multi-operator Prony not implemented!") 353 354 array = np.empty([N, N], dtype="object") 355 new_content = [] 356 for t in range(self.T): 357 new_content.append(array.copy()) 358 359 def wrap(i): 360 while i >= self.T: 361 i -= self.T 362 return i 363 364 for t in range(self.T): 365 for i in range(N): 366 for j in range(N): 367 if periodic: 368 new_content[t][i, j] = self.content[wrap(t + i + j)][0] 369 elif (t + i + j) >= self.T: 370 new_content[t] = None 371 else: 372 new_content[t][i, j] = self.content[t + i + j][0] 373 374 return Corr(new_content) 375 376 def roll(self, dt): 377 """Periodically shift the correlator by dt timeslices 378 379 Parameters 380 ---------- 381 dt : int 382 number of timeslices 383 """ 384 return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) 385 386 def reverse(self): 387 """Reverse the time ordering of the Corr""" 388 return Corr(self.content[:: -1]) 389 390 def thin(self, spacing=2, offset=0): 391 """Thin out a correlator to suppress correlations 392 393 Parameters 394 ---------- 395 spacing : int 396 Keep only every 'spacing'th entry of the correlator 397 offset : int 398 Offset the equal spacing 399 """ 400 new_content = [] 401 for t in range(self.T): 402 if (offset + t) % spacing != 0: 403 new_content.append(None) 404 else: 405 new_content.append(self.content[t]) 406 return Corr(new_content) 407 408 def correlate(self, partner): 409 """Correlate the correlator with another correlator or Obs 410 411 Parameters 412 ---------- 413 partner : Obs or Corr 414 partner to correlate the correlator with. 415 Can either be an Obs which is correlated with all entries of the 416 correlator or a Corr of same length. 417 """ 418 new_content = [] 419 for x0, t_slice in enumerate(self.content): 420 if t_slice is None: 421 new_content.append(None) 422 else: 423 if isinstance(partner, Corr): 424 if partner.content[x0] is None: 425 new_content.append(None) 426 else: 427 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) 428 elif isinstance(partner, Obs): # Should this include CObs? 429 new_content.append(np.array([correlate(o, partner) for o in t_slice])) 430 else: 431 raise Exception("Can only correlate with an Obs or a Corr.") 432 433 return Corr(new_content) 434 435 def reweight(self, weight, **kwargs): 436 """Reweight the correlator. 437 438 Parameters 439 ---------- 440 weight : Obs 441 Reweighting factor. An Observable that has to be defined on a superset of the 442 configurations in obs[i].idl for all i. 443 all_configs : bool 444 if True, the reweighted observables are normalized by the average of 445 the reweighting factor on all configurations in weight.idl and not 446 on the configurations in obs[i].idl. 447 """ 448 new_content = [] 449 for t_slice in self.content: 450 if t_slice is None: 451 new_content.append(None) 452 else: 453 new_content.append(np.array(reweight(weight, t_slice, **kwargs))) 454 return Corr(new_content) 455 456 def T_symmetry(self, partner, parity=+1): 457 """Return the time symmetry average of the correlator and its partner 458 459 Parameters 460 ---------- 461 partner : Corr 462 Time symmetry partner of the Corr 463 partity : int 464 Parity quantum number of the correlator, can be +1 or -1 465 """ 466 if not isinstance(partner, Corr): 467 raise Exception("T partner has to be a Corr object.") 468 if parity not in [+1, -1]: 469 raise Exception("Parity has to be +1 or -1.") 470 T_partner = parity * partner.reverse() 471 472 t_slices = [] 473 test = (self - T_partner) 474 test.gamma_method() 475 for x0, t_slice in enumerate(test.content): 476 if t_slice is not None: 477 if not t_slice[0].is_zero_within_error(5): 478 t_slices.append(x0) 479 if t_slices: 480 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning) 481 482 return (self + T_partner) / 2 483 484 def deriv(self, variant="symmetric"): 485 """Return the first derivative of the correlator with respect to x0. 486 487 Parameters 488 ---------- 489 variant : str 490 decides which definition of the finite differences derivative is used. 491 Available choice: symmetric, forward, backward, improved, default: symmetric 492 """ 493 if variant == "symmetric": 494 newcontent = [] 495 for t in range(1, self.T - 1): 496 if (self.content[t - 1] is None) or (self.content[t + 1] is None): 497 newcontent.append(None) 498 else: 499 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) 500 if(all([x is None for x in newcontent])): 501 raise Exception('Derivative is undefined at all timeslices') 502 return Corr(newcontent, padding=[1, 1]) 503 elif variant == "forward": 504 newcontent = [] 505 for t in range(self.T - 1): 506 if (self.content[t] is None) or (self.content[t + 1] is None): 507 newcontent.append(None) 508 else: 509 newcontent.append(self.content[t + 1] - self.content[t]) 510 if(all([x is None for x in newcontent])): 511 raise Exception("Derivative is undefined at all timeslices") 512 return Corr(newcontent, padding=[0, 1]) 513 elif variant == "backward": 514 newcontent = [] 515 for t in range(1, self.T): 516 if (self.content[t - 1] is None) or (self.content[t] is None): 517 newcontent.append(None) 518 else: 519 newcontent.append(self.content[t] - self.content[t - 1]) 520 if(all([x is None for x in newcontent])): 521 raise Exception("Derivative is undefined at all timeslices") 522 return Corr(newcontent, padding=[1, 0]) 523 elif variant == "improved": 524 newcontent = [] 525 for t in range(2, self.T - 2): 526 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None): 527 newcontent.append(None) 528 else: 529 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2])) 530 if(all([x is None for x in newcontent])): 531 raise Exception('Derivative is undefined at all timeslices') 532 return Corr(newcontent, padding=[2, 2]) 533 else: 534 raise Exception("Unknown variant.") 535 536 def second_deriv(self, variant="symmetric"): 537 """Return the second derivative of the correlator with respect to x0. 538 539 Parameters 540 ---------- 541 variant : str 542 decides which definition of the finite differences derivative is used. 543 Available choice: symmetric, improved, default: symmetric 544 """ 545 if variant == "symmetric": 546 newcontent = [] 547 for t in range(1, self.T - 1): 548 if (self.content[t - 1] is None) or (self.content[t + 1] is None): 549 newcontent.append(None) 550 else: 551 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) 552 if(all([x is None for x in newcontent])): 553 raise Exception("Derivative is undefined at all timeslices") 554 return Corr(newcontent, padding=[1, 1]) 555 elif variant == "improved": 556 newcontent = [] 557 for t in range(2, self.T - 2): 558 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None): 559 newcontent.append(None) 560 else: 561 newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2])) 562 if(all([x is None for x in newcontent])): 563 raise Exception("Derivative is undefined at all timeslices") 564 return Corr(newcontent, padding=[2, 2]) 565 else: 566 raise Exception("Unknown variant.") 567 568 def m_eff(self, variant='log', guess=1.0): 569 """Returns the effective mass of the correlator as correlator object 570 571 Parameters 572 ---------- 573 variant : str 574 log : uses the standard effective mass log(C(t) / C(t+1)) 575 cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. 576 sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m. 577 See, e.g., arXiv:1205.5380 578 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) 579 guess : float 580 guess for the root finder, only relevant for the root variant 581 """ 582 if self.N != 1: 583 raise Exception('Correlator must be projected before getting m_eff') 584 if variant == 'log': 585 newcontent = [] 586 for t in range(self.T - 1): 587 if (self.content[t] is None) or (self.content[t + 1] is None): 588 newcontent.append(None) 589 else: 590 newcontent.append(self.content[t] / self.content[t + 1]) 591 if(all([x is None for x in newcontent])): 592 raise Exception('m_eff is undefined at all timeslices') 593 594 return np.log(Corr(newcontent, padding=[0, 1])) 595 596 elif variant in ['periodic', 'cosh', 'sinh']: 597 if variant in ['periodic', 'cosh']: 598 func = anp.cosh 599 else: 600 func = anp.sinh 601 602 def root_function(x, d): 603 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d 604 605 newcontent = [] 606 for t in range(self.T - 1): 607 if (self.content[t] is None) or (self.content[t + 1] is None): 608 newcontent.append(None) 609 # Fill the two timeslices in the middle of the lattice with their predecessors 610 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: 611 newcontent.append(newcontent[-1]) 612 else: 613 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) 614 if(all([x is None for x in newcontent])): 615 raise Exception('m_eff is undefined at all timeslices') 616 617 return Corr(newcontent, padding=[0, 1]) 618 619 elif variant == 'arccosh': 620 newcontent = [] 621 for t in range(1, self.T - 1): 622 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None): 623 newcontent.append(None) 624 else: 625 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) 626 if(all([x is None for x in newcontent])): 627 raise Exception("m_eff is undefined at all timeslices") 628 return np.arccosh(Corr(newcontent, padding=[1, 1])) 629 630 else: 631 raise Exception('Unknown variant.') 632 633 def fit(self, function, fitrange=None, silent=False, **kwargs): 634 r'''Fits function to the data 635 636 Parameters 637 ---------- 638 function : obj 639 function to fit to the data. See fits.least_squares for details. 640 fitrange : list 641 Two element list containing the timeslices on which the fit is supposed to start and stop. 642 Caution: This range is inclusive as opposed to standard python indexing. 643 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. 644 If not specified, self.prange or all timeslices are used. 645 silent : bool 646 Decides whether output is printed to the standard output. 647 ''' 648 if self.N != 1: 649 raise Exception("Correlator must be projected before fitting") 650 651 if fitrange is None: 652 if self.prange: 653 fitrange = self.prange 654 else: 655 fitrange = [0, self.T - 1] 656 else: 657 if not isinstance(fitrange, list): 658 raise Exception("fitrange has to be a list with two elements") 659 if len(fitrange) != 2: 660 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") 661 662 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] 663 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] 664 result = least_squares(xs, ys, function, silent=silent, **kwargs) 665 return result 666 667 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): 668 """ Extract a plateau value from a Corr object 669 670 Parameters 671 ---------- 672 plateau_range : list 673 list with two entries, indicating the first and the last timeslice 674 of the plateau region. 675 method : str 676 method to extract the plateau. 677 'fit' fits a constant to the plateau region 678 'avg', 'average' or 'mean' just average over the given timeslices. 679 auto_gamma : bool 680 apply gamma_method with default parameters to the Corr. Defaults to None 681 """ 682 if not plateau_range: 683 if self.prange: 684 plateau_range = self.prange 685 else: 686 raise Exception("no plateau range provided") 687 if self.N != 1: 688 raise Exception("Correlator must be projected before getting a plateau.") 689 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): 690 raise Exception("plateau is undefined at all timeslices in plateaurange.") 691 if auto_gamma: 692 self.gamma_method() 693 if method == "fit": 694 def const_func(a, t): 695 return a[0] 696 return self.fit(const_func, plateau_range)[0] 697 elif method in ["avg", "average", "mean"]: 698 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) 699 return returnvalue 700 701 else: 702 raise Exception("Unsupported plateau method: " + method) 703 704 def set_prange(self, prange): 705 """Sets the attribute prange of the Corr object.""" 706 if not len(prange) == 2: 707 raise Exception("prange must be a list or array with two values") 708 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): 709 raise Exception("Start and end point must be integers") 710 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): 711 raise Exception("Start and end point must define a range in the interval 0,T") 712 713 self.prange = prange 714 return 715 716 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None): 717 """Plots the correlator using the tag of the correlator as label if available. 718 719 Parameters 720 ---------- 721 x_range : list 722 list of two values, determining the range of the x-axis e.g. [4, 8] 723 comp : Corr or list of Corr 724 Correlator or list of correlators which are plotted for comparison. 725 The tags of these correlators are used as labels if available. 726 logscale : bool 727 Sets y-axis to logscale 728 plateau : Obs 729 Plateau value to be visualized in the figure 730 fit_res : Fit_result 731 Fit_result object to be visualized 732 ylabel : str 733 Label for the y-axis 734 save : str 735 path to file in which the figure should be saved 736 auto_gamma : bool 737 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. 738 hide_sigma : float 739 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. 740 references : list 741 List of floating point values that are displayed as horizontal lines for reference. 742 """ 743 if self.N != 1: 744 raise Exception("Correlator must be projected before plotting") 745 746 if auto_gamma: 747 self.gamma_method() 748 749 if x_range is None: 750 x_range = [0, self.T - 1] 751 752 fig = plt.figure() 753 ax1 = fig.add_subplot(111) 754 755 x, y, y_err = self.plottable() 756 if hide_sigma: 757 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1 758 else: 759 hide_from = None 760 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) 761 if logscale: 762 ax1.set_yscale('log') 763 else: 764 if y_range is None: 765 try: 766 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)]) 767 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)]) 768 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 769 except Exception: 770 pass 771 else: 772 ax1.set_ylim(y_range) 773 if comp: 774 if isinstance(comp, (Corr, list)): 775 for corr in comp if isinstance(comp, list) else [comp]: 776 if auto_gamma: 777 corr.gamma_method() 778 x, y, y_err = corr.plottable() 779 if hide_sigma: 780 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1 781 else: 782 hide_from = None 783 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 784 else: 785 raise Exception("'comp' must be a correlator or a list of correlators.") 786 787 if plateau: 788 if isinstance(plateau, Obs): 789 if auto_gamma: 790 plateau.gamma_method() 791 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 792 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 793 else: 794 raise Exception("'plateau' must be an Obs") 795 796 if references: 797 if isinstance(references, list): 798 for ref in references: 799 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') 800 else: 801 raise Exception("'references' must be a list of floating pint values.") 802 803 if self.prange: 804 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 805 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 806 807 if fit_res: 808 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 809 ax1.plot(x_samples, 810 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 811 ls='-', marker=',', lw=2) 812 813 ax1.set_xlabel(r'$x_0 / a$') 814 if ylabel: 815 ax1.set_ylabel(ylabel) 816 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 817 818 handles, labels = ax1.get_legend_handles_labels() 819 if labels: 820 ax1.legend() 821 plt.draw() 822 823 if save: 824 if isinstance(save, str): 825 fig.savefig(save) 826 else: 827 raise Exception("'save' has to be a string.") 828 829 def spaghetti_plot(self, logscale=True): 830 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 831 832 Parameters 833 ---------- 834 logscale : bool 835 Determines whether the scale of the y-axis is logarithmic or standard. 836 """ 837 if self.N != 1: 838 raise Exception("Correlator needs to be projected first.") 839 840 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist])) 841 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 842 843 for name in mc_names: 844 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 845 846 fig = plt.figure() 847 ax = fig.add_subplot(111) 848 for dat in data: 849 ax.plot(x0_vals, dat, ls='-', marker='') 850 851 if logscale is True: 852 ax.set_yscale('log') 853 854 ax.set_xlabel(r'$x_0 / a$') 855 plt.title(name) 856 plt.draw() 857 858 def dump(self, filename, datatype="json.gz", **kwargs): 859 """Dumps the Corr into a file of chosen type 860 Parameters 861 ---------- 862 filename : str 863 Name of the file to be saved. 864 datatype : str 865 Format of the exported file. Supported formats include 866 "json.gz" and "pickle" 867 path : str 868 specifies a custom path for the file (default '.') 869 """ 870 if datatype == "json.gz": 871 from .input.json import dump_to_json 872 if 'path' in kwargs: 873 file_name = kwargs.get('path') + '/' + filename 874 else: 875 file_name = filename 876 dump_to_json(self, file_name) 877 elif datatype == "pickle": 878 dump_object(self, filename, **kwargs) 879 else: 880 raise Exception("Unknown datatype " + str(datatype)) 881 882 def print(self, range=[0, None]): 883 print(self.__repr__(range)) 884 885 def __repr__(self, range=[0, None]): 886 content_string = "" 887 888 content_string += "Corr T=" + str(self.T) + " N=" + str(self.N) + "\n" # +" filled with"+ str(type(self.content[0][0])) there should be a good solution here 889 890 if self.tag is not None: 891 content_string += "Description: " + self.tag + "\n" 892 if self.N != 1: 893 return content_string 894 895 if range[1]: 896 range[1] += 1 897 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' 898 for i, sub_corr in enumerate(self.content[range[0]:range[1]]): 899 if sub_corr is None: 900 content_string += str(i + range[0]) + '\n' 901 else: 902 content_string += str(i + range[0]) 903 for element in sub_corr: 904 content_string += '\t' + ' ' * int(element >= 0) + str(element) 905 content_string += '\n' 906 return content_string 907 908 def __str__(self): 909 return self.__repr__() 910 911 # We define the basic operations, that can be performed with correlators. 912 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. 913 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. 914 # One could try and tell Obs to check if the y in __mul__ is a Corr and 915 916 def __add__(self, y): 917 if isinstance(y, Corr): 918 if ((self.N != y.N) or (self.T != y.T)): 919 raise Exception("Addition of Corrs with different shape") 920 newcontent = [] 921 for t in range(self.T): 922 if (self.content[t] is None) or (y.content[t] is None): 923 newcontent.append(None) 924 else: 925 newcontent.append(self.content[t] + y.content[t]) 926 return Corr(newcontent) 927 928 elif isinstance(y, (Obs, int, float, CObs)): 929 newcontent = [] 930 for t in range(self.T): 931 if (self.content[t] is None): 932 newcontent.append(None) 933 else: 934 newcontent.append(self.content[t] + y) 935 return Corr(newcontent, prange=self.prange) 936 elif isinstance(y, np.ndarray): 937 if y.shape == (self.T,): 938 return Corr(list((np.array(self.content).T + y).T)) 939 else: 940 raise ValueError("operands could not be broadcast together") 941 else: 942 raise TypeError("Corr + wrong type") 943 944 def __mul__(self, y): 945 if isinstance(y, Corr): 946 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 947 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 948 newcontent = [] 949 for t in range(self.T): 950 if (self.content[t] is None) or (y.content[t] is None): 951 newcontent.append(None) 952 else: 953 newcontent.append(self.content[t] * y.content[t]) 954 return Corr(newcontent) 955 956 elif isinstance(y, (Obs, int, float, CObs)): 957 newcontent = [] 958 for t in range(self.T): 959 if (self.content[t] is None): 960 newcontent.append(None) 961 else: 962 newcontent.append(self.content[t] * y) 963 return Corr(newcontent, prange=self.prange) 964 elif isinstance(y, np.ndarray): 965 if y.shape == (self.T,): 966 return Corr(list((np.array(self.content).T * y).T)) 967 else: 968 raise ValueError("operands could not be broadcast together") 969 else: 970 raise TypeError("Corr * wrong type") 971 972 def __truediv__(self, y): 973 if isinstance(y, Corr): 974 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 975 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 976 newcontent = [] 977 for t in range(self.T): 978 if (self.content[t] is None) or (y.content[t] is None): 979 newcontent.append(None) 980 else: 981 newcontent.append(self.content[t] / y.content[t]) 982 for t in range(self.T): 983 if newcontent[t] is None: 984 continue 985 if np.isnan(np.sum(newcontent[t]).value): 986 newcontent[t] = None 987 988 if all([item is None for item in newcontent]): 989 raise Exception("Division returns completely undefined correlator") 990 return Corr(newcontent) 991 992 elif isinstance(y, (Obs, CObs)): 993 if isinstance(y, Obs): 994 if y.value == 0: 995 raise Exception('Division by zero will return undefined correlator') 996 if isinstance(y, CObs): 997 if y.is_zero(): 998 raise Exception('Division by zero will return undefined correlator') 999 1000 newcontent = [] 1001 for t in range(self.T): 1002 if (self.content[t] is None): 1003 newcontent.append(None) 1004 else: 1005 newcontent.append(self.content[t] / y) 1006 return Corr(newcontent, prange=self.prange) 1007 1008 elif isinstance(y, (int, float)): 1009 if y == 0: 1010 raise Exception('Division by zero will return undefined correlator') 1011 newcontent = [] 1012 for t in range(self.T): 1013 if (self.content[t] is None): 1014 newcontent.append(None) 1015 else: 1016 newcontent.append(self.content[t] / y) 1017 return Corr(newcontent, prange=self.prange) 1018 elif isinstance(y, np.ndarray): 1019 if y.shape == (self.T,): 1020 return Corr(list((np.array(self.content).T / y).T)) 1021 else: 1022 raise ValueError("operands could not be broadcast together") 1023 else: 1024 raise TypeError('Corr / wrong type') 1025 1026 def __neg__(self): 1027 newcontent = [None if (item is None) else -1. * item for item in self.content] 1028 return Corr(newcontent, prange=self.prange) 1029 1030 def __sub__(self, y): 1031 return self + (-y) 1032 1033 def __pow__(self, y): 1034 if isinstance(y, (Obs, int, float, CObs)): 1035 newcontent = [None if (item is None) else item**y for item in self.content] 1036 return Corr(newcontent, prange=self.prange) 1037 else: 1038 raise TypeError('Type of exponent not supported') 1039 1040 def __abs__(self): 1041 newcontent = [None if (item is None) else np.abs(item) for item in self.content] 1042 return Corr(newcontent, prange=self.prange) 1043 1044 # The numpy functions: 1045 def sqrt(self): 1046 return self**0.5 1047 1048 def log(self): 1049 newcontent = [None if (item is None) else np.log(item) for item in self.content] 1050 return Corr(newcontent, prange=self.prange) 1051 1052 def exp(self): 1053 newcontent = [None if (item is None) else np.exp(item) for item in self.content] 1054 return Corr(newcontent, prange=self.prange) 1055 1056 def _apply_func_to_corr(self, func): 1057 newcontent = [None if (item is None) else func(item) for item in self.content] 1058 for t in range(self.T): 1059 if newcontent[t] is None: 1060 continue 1061 if np.isnan(np.sum(newcontent[t]).value): 1062 newcontent[t] = None 1063 if all([item is None for item in newcontent]): 1064 raise Exception('Operation returns undefined correlator') 1065 return Corr(newcontent) 1066 1067 def sin(self): 1068 return self._apply_func_to_corr(np.sin) 1069 1070 def cos(self): 1071 return self._apply_func_to_corr(np.cos) 1072 1073 def tan(self): 1074 return self._apply_func_to_corr(np.tan) 1075 1076 def sinh(self): 1077 return self._apply_func_to_corr(np.sinh) 1078 1079 def cosh(self): 1080 return self._apply_func_to_corr(np.cosh) 1081 1082 def tanh(self): 1083 return self._apply_func_to_corr(np.tanh) 1084 1085 def arcsin(self): 1086 return self._apply_func_to_corr(np.arcsin) 1087 1088 def arccos(self): 1089 return self._apply_func_to_corr(np.arccos) 1090 1091 def arctan(self): 1092 return self._apply_func_to_corr(np.arctan) 1093 1094 def arcsinh(self): 1095 return self._apply_func_to_corr(np.arcsinh) 1096 1097 def arccosh(self): 1098 return self._apply_func_to_corr(np.arccosh) 1099 1100 def arctanh(self): 1101 return self._apply_func_to_corr(np.arctanh) 1102 1103 # Right hand side operations (require tweak in main module to work) 1104 def __radd__(self, y): 1105 return self + y 1106 1107 def __rsub__(self, y): 1108 return -self + y 1109 1110 def __rmul__(self, y): 1111 return self * y 1112 1113 def __rtruediv__(self, y): 1114 return (self / y) ** (-1) 1115 1116 @property 1117 def real(self): 1118 def return_real(obs_OR_cobs): 1119 if isinstance(obs_OR_cobs, CObs): 1120 return obs_OR_cobs.real 1121 else: 1122 return obs_OR_cobs 1123 1124 return self._apply_func_to_corr(return_real) 1125 1126 @property 1127 def imag(self): 1128 def return_imag(obs_OR_cobs): 1129 if isinstance(obs_OR_cobs, CObs): 1130 return obs_OR_cobs.imag 1131 else: 1132 return obs_OR_cobs * 0 # So it stays the right type 1133 1134 return self._apply_func_to_corr(return_imag) 1135 1136 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): 1137 r''' Project large correlation matrix to lowest states 1138 1139 This method can be used to reduce the size of an (N x N) correlation matrix 1140 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise 1141 is still small. 1142 1143 Parameters 1144 ---------- 1145 Ntrunc: int 1146 Rank of the target matrix. 1147 tproj: int 1148 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. 1149 The default value is 3. 1150 t0proj: int 1151 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly 1152 discouraged for O(a) improved theories, since the correctness of the procedure 1153 cannot be granted in this case. The default value is 2. 1154 basematrix : Corr 1155 Correlation matrix that is used to determine the eigenvectors of the 1156 lowest states based on a GEVP. basematrix is taken to be the Corr itself if 1157 is is not specified. 1158 1159 Notes 1160 ----- 1161 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving 1162 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$ 1163 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the 1164 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via 1165 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large 1166 correlation matrix and to remove some noise that is added by irrelevant operators. 1167 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated 1168 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. 1169 ''' 1170 1171 if self.N == 1: 1172 raise Exception('Method cannot be applied to one-dimensional correlators.') 1173 if basematrix is None: 1174 basematrix = self 1175 if Ntrunc >= basematrix.N: 1176 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) 1177 if basematrix.N != self.N: 1178 raise Exception('basematrix and targetmatrix have to be of the same size.') 1179 1180 evecs = [] 1181 for i in range(Ntrunc): 1182 evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None)) 1183 1184 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) 1185 rmat = [] 1186 for t in range(basematrix.T): 1187 for i in range(Ntrunc): 1188 for j in range(Ntrunc): 1189 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] 1190 rmat.append(np.copy(tmpmat)) 1191 1192 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] 1193 return Corr(newcontent)
The class for a correlator (time dependent sequence of pe.Obs).
Everything, this class does, can be achieved using lists or arrays of Obs. But it is simply more convenient to have a dedicated object for correlators. One often wants to add or multiply correlators of the same length at every timeslice and it is inconvenient to iterate over all timeslices for every operation. This is especially true, when dealing with matrices.
The correlator can have two types of content: An Obs at every timeslice OR a GEVP matrix at every timeslice. Other dependency (eg. spatial) are not supported.
27 def __init__(self, data_input, padding=[0, 0], prange=None): 28 """ Initialize a Corr object. 29 30 Parameters 31 ---------- 32 data_input : list or array 33 list of Obs or list of arrays of Obs or array of Corrs 34 padding : list, optional 35 List with two entries where the first labels the padding 36 at the front of the correlator and the second the padding 37 at the back. 38 prange : list, optional 39 List containing the first and last timeslice of the plateau 40 region indentified for this correlator. 41 """ 42 43 if isinstance(data_input, np.ndarray): 44 45 # This only works, if the array fulfills the conditions below 46 if not len(data_input.shape) == 2 and data_input.shape[0] == data_input.shape[1]: 47 raise Exception("Incompatible array shape") 48 if not all([isinstance(item, Corr) for item in data_input.flatten()]): 49 raise Exception("If the input is an array, its elements must be of type pe.Corr") 50 if not all([item.N == 1 for item in data_input.flatten()]): 51 raise Exception("Can only construct matrix correlator from single valued correlators") 52 if not len(set([item.T for item in data_input.flatten()])) == 1: 53 raise Exception("All input Correlators must be defined over the same timeslices.") 54 55 T = data_input[0, 0].T 56 N = data_input.shape[0] 57 input_as_list = [] 58 for t in range(T): 59 if any([(item.content[t] is None) for item in data_input.flatten()]): 60 if not all([(item.content[t] is None) for item in data_input.flatten()]): 61 warnings.warn("Input ill-defined at different timeslices. Conversion leads to data loss!", RuntimeWarning) 62 input_as_list.append(None) 63 else: 64 array_at_timeslace = np.empty([N, N], dtype="object") 65 for i in range(N): 66 for j in range(N): 67 array_at_timeslace[i, j] = data_input[i, j][t] 68 input_as_list.append(array_at_timeslace) 69 data_input = input_as_list 70 71 if isinstance(data_input, list): 72 73 if all([isinstance(item, (Obs, CObs)) for item in data_input]): 74 _assert_equal_properties(data_input) 75 self.content = [np.asarray([item]) for item in data_input] 76 if all([isinstance(item, (Obs, CObs)) or item is None for item in data_input]): 77 _assert_equal_properties([o for o in data_input if o is not None]) 78 self.content = [np.asarray([item]) if item is not None else None for item in data_input] 79 self.N = 1 80 81 elif all([isinstance(item, np.ndarray) or item is None for item in data_input]) and any([isinstance(item, np.ndarray) for item in data_input]): 82 self.content = data_input 83 noNull = [a for a in self.content if not (a is None)] # To check if the matrices are correct for all undefined elements 84 self.N = noNull[0].shape[0] 85 if self.N > 1 and noNull[0].shape[0] != noNull[0].shape[1]: 86 raise Exception("Smearing matrices are not NxN") 87 if (not all([item.shape == noNull[0].shape for item in noNull])): 88 raise Exception("Items in data_input are not of identical shape." + str(noNull)) 89 else: 90 raise Exception("data_input contains item of wrong type") 91 else: 92 raise Exception("Data input was not given as list or correct array") 93 94 self.tag = None 95 96 # An undefined timeslice is represented by the None object 97 self.content = [None] * padding[0] + self.content + [None] * padding[1] 98 self.T = len(self.content) 99 self.prange = prange
Initialize a Corr object.
Parameters
- data_input (list or array): list of Obs or list of arrays of Obs or array of Corrs
- padding (list, optional): List with two entries where the first labels the padding at the front of the correlator and the second the padding at the back.
- prange (list, optional): List containing the first and last timeslice of the plateau region indentified for this correlator.
120 def gamma_method(self, **kwargs): 121 """Apply the gamma method to the content of the Corr.""" 122 for item in self.content: 123 if not(item is None): 124 if self.N == 1: 125 item[0].gamma_method(**kwargs) 126 else: 127 for i in range(self.N): 128 for j in range(self.N): 129 item[i, j].gamma_method(**kwargs)
Apply the gamma method to the content of the Corr.
131 def projected(self, vector_l=None, vector_r=None, normalize=False): 132 """We need to project the Correlator with a Vector to get a single value at each timeslice. 133 134 The method can use one or two vectors. 135 If two are specified it returns v1@G@v2 (the order might be very important.) 136 By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to 137 """ 138 if self.N == 1: 139 raise Exception("Trying to project a Corr, that already has N=1.") 140 141 if vector_l is None: 142 vector_l, vector_r = np.asarray([1.] + (self.N - 1) * [0.]), np.asarray([1.] + (self.N - 1) * [0.]) 143 elif(vector_r is None): 144 vector_r = vector_l 145 if isinstance(vector_l, list) and not isinstance(vector_r, list): 146 if len(vector_l) != self.T: 147 raise Exception("Length of vector list must be equal to T") 148 vector_r = [vector_r] * self.T 149 if isinstance(vector_r, list) and not isinstance(vector_l, list): 150 if len(vector_r) != self.T: 151 raise Exception("Length of vector list must be equal to T") 152 vector_l = [vector_l] * self.T 153 154 if not isinstance(vector_l, list): 155 if not vector_l.shape == vector_r.shape == (self.N,): 156 raise Exception("Vectors are of wrong shape!") 157 if normalize: 158 vector_l, vector_r = vector_l / np.sqrt((vector_l @ vector_l)), vector_r / np.sqrt(vector_r @ vector_r) 159 newcontent = [None if (item is None) else np.asarray([vector_l.T @ item @ vector_r]) for item in self.content] 160 161 else: 162 # There are no checks here yet. There are so many possible scenarios, where this can go wrong. 163 if normalize: 164 for t in range(self.T): 165 vector_l[t], vector_r[t] = vector_l[t] / np.sqrt((vector_l[t] @ vector_l[t])), vector_r[t] / np.sqrt(vector_r[t] @ vector_r[t]) 166 167 newcontent = [None if (self.content[t] is None or vector_l[t] is None or vector_r[t] is None) else np.asarray([vector_l[t].T @ self.content[t] @ vector_r[t]]) for t in range(self.T)] 168 return Corr(newcontent)
We need to project the Correlator with a Vector to get a single value at each timeslice.
The method can use one or two vectors. If two are specified it returns v1@G@v2 (the order might be very important.) By default it will return the lowest source, which usually means unsmeared-unsmeared (0,0), but it does not have to
170 def item(self, i, j): 171 """Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice. 172 173 Parameters 174 ---------- 175 i : int 176 First index to be picked. 177 j : int 178 Second index to be picked. 179 """ 180 if self.N == 1: 181 raise Exception("Trying to pick item from projected Corr") 182 newcontent = [None if(item is None) else item[i, j] for item in self.content] 183 return Corr(newcontent)
Picks the element [i,j] from every matrix and returns a correlator containing one Obs per timeslice.
Parameters
- i (int): First index to be picked.
- j (int): Second index to be picked.
185 def plottable(self): 186 """Outputs the correlator in a plotable format. 187 188 Outputs three lists containing the timeslice index, the value on each 189 timeslice and the error on each timeslice. 190 """ 191 if self.N != 1: 192 raise Exception("Can only make Corr[N=1] plottable") 193 x_list = [x for x in range(self.T) if not self.content[x] is None] 194 y_list = [y[0].value for y in self.content if y is not None] 195 y_err_list = [y[0].dvalue for y in self.content if y is not None] 196 197 return x_list, y_list, y_err_list
Outputs the correlator in a plotable format.
Outputs three lists containing the timeslice index, the value on each timeslice and the error on each timeslice.
199 def symmetric(self): 200 """ Symmetrize the correlator around x0=0.""" 201 if self.T % 2 != 0: 202 raise Exception("Can not symmetrize odd T") 203 204 if np.argmax(np.abs(self.content)) != 0: 205 warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning) 206 207 newcontent = [self.content[0]] 208 for t in range(1, self.T): 209 if (self.content[t] is None) or (self.content[self.T - t] is None): 210 newcontent.append(None) 211 else: 212 newcontent.append(0.5 * (self.content[t] + self.content[self.T - t])) 213 if(all([x is None for x in newcontent])): 214 raise Exception("Corr could not be symmetrized: No redundant values") 215 return Corr(newcontent, prange=self.prange)
Symmetrize the correlator around x0=0.
217 def anti_symmetric(self): 218 """Anti-symmetrize the correlator around x0=0.""" 219 if self.T % 2 != 0: 220 raise Exception("Can not symmetrize odd T") 221 222 test = 1 * self 223 test.gamma_method() 224 if not all([o.is_zero_within_error(3) for o in test.content[0]]): 225 warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning) 226 227 newcontent = [self.content[0]] 228 for t in range(1, self.T): 229 if (self.content[t] is None) or (self.content[self.T - t] is None): 230 newcontent.append(None) 231 else: 232 newcontent.append(0.5 * (self.content[t] - self.content[self.T - t])) 233 if(all([x is None for x in newcontent])): 234 raise Exception("Corr could not be symmetrized: No redundant values") 235 return Corr(newcontent, prange=self.prange)
Anti-symmetrize the correlator around x0=0.
237 def matrix_symmetric(self): 238 """Symmetrizes the correlator matrices on every timeslice.""" 239 if self.N > 1: 240 transposed = [None if (G is None) else G.T for G in self.content] 241 return 0.5 * (Corr(transposed) + self) 242 if self.N == 1: 243 raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
Symmetrizes the correlator matrices on every timeslice.
245 def GEVP(self, t0, ts=None, state=0, sorted_list="Eigenvalue"): 246 """Solve the general eigenvalue problem on the current correlator 247 248 Parameters 249 ---------- 250 t0 : int 251 The time t0 for G(t)v= lambda G(t_0)v 252 ts : int 253 fixed time G(t_s)v= lambda G(t_0)v if return_list=False 254 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. 255 state : int 256 The state one is interested in, ordered by energy. The lowest state is zero. 257 sorted_list : string 258 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. 259 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. 260 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. 261 The reference state is identified by its eigenvalue at t=ts 262 """ 263 264 if self.N == 1: 265 raise Exception("GEVP methods only works on correlator matrices and not single correlators.") 266 267 symmetric_corr = self.matrix_symmetric() 268 if sorted_list is None: 269 if (ts is None): 270 raise Exception("ts is required if sorted_list=None.") 271 if (ts <= t0): 272 raise Exception("ts has to be larger than t0.") 273 if (self.content[t0] is None) or (self.content[ts] is None): 274 raise Exception("Corr not defined at t0/ts.") 275 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") 276 for i in range(self.N): 277 for j in range(self.N): 278 G0[i, j] = symmetric_corr[t0][i, j].value 279 Gt[i, j] = symmetric_corr[ts][i, j].value 280 281 sp_vecs = _GEVP_solver(Gt, G0) 282 sp_vec = sp_vecs[state] 283 return sp_vec 284 elif sorted_list in ["Eigenvalue", "Eigenvector"]: 285 if sorted_list == "Eigenvalue" and ts is not None: 286 warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning) 287 all_vecs = [None] * (t0 + 1) 288 for t in range(t0 + 1, self.T): 289 try: 290 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") 291 for i in range(self.N): 292 for j in range(self.N): 293 G0[i, j] = symmetric_corr[t0][i, j].value 294 Gt[i, j] = symmetric_corr[t][i, j].value 295 296 sp_vecs = _GEVP_solver(Gt, G0) 297 if sorted_list == "Eigenvalue": 298 sp_vec = sp_vecs[state] 299 all_vecs.append(sp_vec) 300 else: 301 all_vecs.append(sp_vecs) 302 except Exception: 303 all_vecs.append(None) 304 if sorted_list == "Eigenvector": 305 if (ts is None): 306 raise Exception("ts is required for the Eigenvector sorting method.") 307 all_vecs = _sort_vectors(all_vecs, ts) 308 all_vecs = [a[state] if a is not None else None for a in all_vecs] 309 else: 310 raise Exception("Unkown value for 'sorted_list'.") 311 312 return all_vecs
Solve the general eigenvalue problem on the current correlator
Parameters
- t0 (int): The time t0 for G(t)v= lambda G(t_0)v
- ts (int): fixed time G(t_s)v= lambda G(t_0)v if return_list=False If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
- state (int): The state one is interested in, ordered by energy. The lowest state is zero.
- sorted_list (string): 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. "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at t=ts
314 def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None): 315 """Determines the eigenvalue of the GEVP by solving and projecting the correlator 316 317 Parameters 318 ---------- 319 t0 : int 320 The time t0 for G(t)v= lambda G(t_0)v 321 ts : int 322 fixed time G(t_s)v= lambda G(t_0)v if return_list=False 323 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. 324 state : int 325 The state one is interested in ordered by energy. The lowest state is zero. 326 sorted_list : string 327 if this argument is set, a list of vectors (len=self.T) is returned. If it is left as None, only one vector is returned. 328 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. 329 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. 330 The reference state is identified by its eigenvalue at t=ts 331 """ 332 vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list) 333 return self.projected(vec)
Determines the eigenvalue of the GEVP by solving and projecting the correlator
Parameters
- t0 (int): The time t0 for G(t)v= lambda G(t_0)v
- ts (int): fixed time G(t_s)v= lambda G(t_0)v if return_list=False If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method.
- state (int): The state one is interested in ordered by energy. The lowest state is zero.
- sorted_list (string): 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. "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. The reference state is identified by its eigenvalue at t=ts
335 def Hankel(self, N, periodic=False): 336 """Constructs an NxN Hankel matrix 337 338 C(t) c(t+1) ... c(t+n-1) 339 C(t+1) c(t+2) ... c(t+n) 340 ................. 341 C(t+(n-1)) c(t+n) ... c(t+2(n-1)) 342 343 Parameters 344 ---------- 345 N : int 346 Dimension of the Hankel matrix 347 periodic : bool, optional 348 determines whether the matrix is extended periodically 349 """ 350 351 if self.N != 1: 352 raise Exception("Multi-operator Prony not implemented!") 353 354 array = np.empty([N, N], dtype="object") 355 new_content = [] 356 for t in range(self.T): 357 new_content.append(array.copy()) 358 359 def wrap(i): 360 while i >= self.T: 361 i -= self.T 362 return i 363 364 for t in range(self.T): 365 for i in range(N): 366 for j in range(N): 367 if periodic: 368 new_content[t][i, j] = self.content[wrap(t + i + j)][0] 369 elif (t + i + j) >= self.T: 370 new_content[t] = None 371 else: 372 new_content[t][i, j] = self.content[t + i + j][0] 373 374 return Corr(new_content)
Constructs an NxN Hankel matrix
C(t) c(t+1) ... c(t+n-1) C(t+1) c(t+2) ... c(t+n) ................. C(t+(n-1)) c(t+n) ... c(t+2(n-1))
Parameters
- N (int): Dimension of the Hankel matrix
- periodic (bool, optional): determines whether the matrix is extended periodically
376 def roll(self, dt): 377 """Periodically shift the correlator by dt timeslices 378 379 Parameters 380 ---------- 381 dt : int 382 number of timeslices 383 """ 384 return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
Periodically shift the correlator by dt timeslices
Parameters
- dt (int): number of timeslices
386 def reverse(self): 387 """Reverse the time ordering of the Corr""" 388 return Corr(self.content[:: -1])
Reverse the time ordering of the Corr
390 def thin(self, spacing=2, offset=0): 391 """Thin out a correlator to suppress correlations 392 393 Parameters 394 ---------- 395 spacing : int 396 Keep only every 'spacing'th entry of the correlator 397 offset : int 398 Offset the equal spacing 399 """ 400 new_content = [] 401 for t in range(self.T): 402 if (offset + t) % spacing != 0: 403 new_content.append(None) 404 else: 405 new_content.append(self.content[t]) 406 return Corr(new_content)
Thin out a correlator to suppress correlations
Parameters
- spacing (int): Keep only every 'spacing'th entry of the correlator
- offset (int): Offset the equal spacing
408 def correlate(self, partner): 409 """Correlate the correlator with another correlator or Obs 410 411 Parameters 412 ---------- 413 partner : Obs or Corr 414 partner to correlate the correlator with. 415 Can either be an Obs which is correlated with all entries of the 416 correlator or a Corr of same length. 417 """ 418 new_content = [] 419 for x0, t_slice in enumerate(self.content): 420 if t_slice is None: 421 new_content.append(None) 422 else: 423 if isinstance(partner, Corr): 424 if partner.content[x0] is None: 425 new_content.append(None) 426 else: 427 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) 428 elif isinstance(partner, Obs): # Should this include CObs? 429 new_content.append(np.array([correlate(o, partner) for o in t_slice])) 430 else: 431 raise Exception("Can only correlate with an Obs or a Corr.") 432 433 return Corr(new_content)
Correlate the correlator with another correlator or Obs
Parameters
- partner (Obs or Corr): partner to correlate the correlator with. Can either be an Obs which is correlated with all entries of the correlator or a Corr of same length.
435 def reweight(self, weight, **kwargs): 436 """Reweight the correlator. 437 438 Parameters 439 ---------- 440 weight : Obs 441 Reweighting factor. An Observable that has to be defined on a superset of the 442 configurations in obs[i].idl for all i. 443 all_configs : bool 444 if True, the reweighted observables are normalized by the average of 445 the reweighting factor on all configurations in weight.idl and not 446 on the configurations in obs[i].idl. 447 """ 448 new_content = [] 449 for t_slice in self.content: 450 if t_slice is None: 451 new_content.append(None) 452 else: 453 new_content.append(np.array(reweight(weight, t_slice, **kwargs))) 454 return Corr(new_content)
Reweight the correlator.
Parameters
- weight (Obs): Reweighting factor. An Observable that has to be defined on a superset of the configurations in obs[i].idl for all i.
- all_configs (bool): if True, the reweighted observables are normalized by the average of the reweighting factor on all configurations in weight.idl and not on the configurations in obs[i].idl.
456 def T_symmetry(self, partner, parity=+1): 457 """Return the time symmetry average of the correlator and its partner 458 459 Parameters 460 ---------- 461 partner : Corr 462 Time symmetry partner of the Corr 463 partity : int 464 Parity quantum number of the correlator, can be +1 or -1 465 """ 466 if not isinstance(partner, Corr): 467 raise Exception("T partner has to be a Corr object.") 468 if parity not in [+1, -1]: 469 raise Exception("Parity has to be +1 or -1.") 470 T_partner = parity * partner.reverse() 471 472 t_slices = [] 473 test = (self - T_partner) 474 test.gamma_method() 475 for x0, t_slice in enumerate(test.content): 476 if t_slice is not None: 477 if not t_slice[0].is_zero_within_error(5): 478 t_slices.append(x0) 479 if t_slices: 480 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning) 481 482 return (self + T_partner) / 2
Return the time symmetry average of the correlator and its partner
Parameters
- partner (Corr): Time symmetry partner of the Corr
- partity (int): Parity quantum number of the correlator, can be +1 or -1
484 def deriv(self, variant="symmetric"): 485 """Return the first derivative of the correlator with respect to x0. 486 487 Parameters 488 ---------- 489 variant : str 490 decides which definition of the finite differences derivative is used. 491 Available choice: symmetric, forward, backward, improved, default: symmetric 492 """ 493 if variant == "symmetric": 494 newcontent = [] 495 for t in range(1, self.T - 1): 496 if (self.content[t - 1] is None) or (self.content[t + 1] is None): 497 newcontent.append(None) 498 else: 499 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) 500 if(all([x is None for x in newcontent])): 501 raise Exception('Derivative is undefined at all timeslices') 502 return Corr(newcontent, padding=[1, 1]) 503 elif variant == "forward": 504 newcontent = [] 505 for t in range(self.T - 1): 506 if (self.content[t] is None) or (self.content[t + 1] is None): 507 newcontent.append(None) 508 else: 509 newcontent.append(self.content[t + 1] - self.content[t]) 510 if(all([x is None for x in newcontent])): 511 raise Exception("Derivative is undefined at all timeslices") 512 return Corr(newcontent, padding=[0, 1]) 513 elif variant == "backward": 514 newcontent = [] 515 for t in range(1, self.T): 516 if (self.content[t - 1] is None) or (self.content[t] is None): 517 newcontent.append(None) 518 else: 519 newcontent.append(self.content[t] - self.content[t - 1]) 520 if(all([x is None for x in newcontent])): 521 raise Exception("Derivative is undefined at all timeslices") 522 return Corr(newcontent, padding=[1, 0]) 523 elif variant == "improved": 524 newcontent = [] 525 for t in range(2, self.T - 2): 526 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None): 527 newcontent.append(None) 528 else: 529 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2])) 530 if(all([x is None for x in newcontent])): 531 raise Exception('Derivative is undefined at all timeslices') 532 return Corr(newcontent, padding=[2, 2]) 533 else: 534 raise Exception("Unknown variant.")
Return the first derivative of the correlator with respect to x0.
Parameters
- variant (str): decides which definition of the finite differences derivative is used. Available choice: symmetric, forward, backward, improved, default: symmetric
536 def second_deriv(self, variant="symmetric"): 537 """Return the second derivative of the correlator with respect to x0. 538 539 Parameters 540 ---------- 541 variant : str 542 decides which definition of the finite differences derivative is used. 543 Available choice: symmetric, improved, default: symmetric 544 """ 545 if variant == "symmetric": 546 newcontent = [] 547 for t in range(1, self.T - 1): 548 if (self.content[t - 1] is None) or (self.content[t + 1] is None): 549 newcontent.append(None) 550 else: 551 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) 552 if(all([x is None for x in newcontent])): 553 raise Exception("Derivative is undefined at all timeslices") 554 return Corr(newcontent, padding=[1, 1]) 555 elif variant == "improved": 556 newcontent = [] 557 for t in range(2, self.T - 2): 558 if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None): 559 newcontent.append(None) 560 else: 561 newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2])) 562 if(all([x is None for x in newcontent])): 563 raise Exception("Derivative is undefined at all timeslices") 564 return Corr(newcontent, padding=[2, 2]) 565 else: 566 raise Exception("Unknown variant.")
Return the second derivative of the correlator with respect to x0.
Parameters
- variant (str): decides which definition of the finite differences derivative is used. Available choice: symmetric, improved, default: symmetric
568 def m_eff(self, variant='log', guess=1.0): 569 """Returns the effective mass of the correlator as correlator object 570 571 Parameters 572 ---------- 573 variant : str 574 log : uses the standard effective mass log(C(t) / C(t+1)) 575 cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m. 576 sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m. 577 See, e.g., arXiv:1205.5380 578 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) 579 guess : float 580 guess for the root finder, only relevant for the root variant 581 """ 582 if self.N != 1: 583 raise Exception('Correlator must be projected before getting m_eff') 584 if variant == 'log': 585 newcontent = [] 586 for t in range(self.T - 1): 587 if (self.content[t] is None) or (self.content[t + 1] is None): 588 newcontent.append(None) 589 else: 590 newcontent.append(self.content[t] / self.content[t + 1]) 591 if(all([x is None for x in newcontent])): 592 raise Exception('m_eff is undefined at all timeslices') 593 594 return np.log(Corr(newcontent, padding=[0, 1])) 595 596 elif variant in ['periodic', 'cosh', 'sinh']: 597 if variant in ['periodic', 'cosh']: 598 func = anp.cosh 599 else: 600 func = anp.sinh 601 602 def root_function(x, d): 603 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d 604 605 newcontent = [] 606 for t in range(self.T - 1): 607 if (self.content[t] is None) or (self.content[t + 1] is None): 608 newcontent.append(None) 609 # Fill the two timeslices in the middle of the lattice with their predecessors 610 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: 611 newcontent.append(newcontent[-1]) 612 else: 613 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) 614 if(all([x is None for x in newcontent])): 615 raise Exception('m_eff is undefined at all timeslices') 616 617 return Corr(newcontent, padding=[0, 1]) 618 619 elif variant == 'arccosh': 620 newcontent = [] 621 for t in range(1, self.T - 1): 622 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None): 623 newcontent.append(None) 624 else: 625 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) 626 if(all([x is None for x in newcontent])): 627 raise Exception("m_eff is undefined at all timeslices") 628 return np.arccosh(Corr(newcontent, padding=[1, 1])) 629 630 else: 631 raise Exception('Unknown variant.')
Returns the effective mass of the correlator as correlator object
Parameters
- variant (str): log : uses the standard effective mass log(C(t) / C(t+1)) 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. 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. See, e.g., arXiv:1205.5380 arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
- guess (float): guess for the root finder, only relevant for the root variant
633 def fit(self, function, fitrange=None, silent=False, **kwargs): 634 r'''Fits function to the data 635 636 Parameters 637 ---------- 638 function : obj 639 function to fit to the data. See fits.least_squares for details. 640 fitrange : list 641 Two element list containing the timeslices on which the fit is supposed to start and stop. 642 Caution: This range is inclusive as opposed to standard python indexing. 643 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. 644 If not specified, self.prange or all timeslices are used. 645 silent : bool 646 Decides whether output is printed to the standard output. 647 ''' 648 if self.N != 1: 649 raise Exception("Correlator must be projected before fitting") 650 651 if fitrange is None: 652 if self.prange: 653 fitrange = self.prange 654 else: 655 fitrange = [0, self.T - 1] 656 else: 657 if not isinstance(fitrange, list): 658 raise Exception("fitrange has to be a list with two elements") 659 if len(fitrange) != 2: 660 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") 661 662 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] 663 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] 664 result = least_squares(xs, ys, function, silent=silent, **kwargs) 665 return result
Fits function to the data
Parameters
- function (obj): function to fit to the data. See fits.least_squares for details.
- fitrange (list):
Two element list containing the timeslices on which the fit is supposed to start and stop.
Caution: This range is inclusive as opposed to standard python indexing.
fitrange=[4, 6]
corresponds to the three entries 4, 5 and 6. If not specified, self.prange or all timeslices are used. - silent (bool): Decides whether output is printed to the standard output.
667 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): 668 """ Extract a plateau value from a Corr object 669 670 Parameters 671 ---------- 672 plateau_range : list 673 list with two entries, indicating the first and the last timeslice 674 of the plateau region. 675 method : str 676 method to extract the plateau. 677 'fit' fits a constant to the plateau region 678 'avg', 'average' or 'mean' just average over the given timeslices. 679 auto_gamma : bool 680 apply gamma_method with default parameters to the Corr. Defaults to None 681 """ 682 if not plateau_range: 683 if self.prange: 684 plateau_range = self.prange 685 else: 686 raise Exception("no plateau range provided") 687 if self.N != 1: 688 raise Exception("Correlator must be projected before getting a plateau.") 689 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): 690 raise Exception("plateau is undefined at all timeslices in plateaurange.") 691 if auto_gamma: 692 self.gamma_method() 693 if method == "fit": 694 def const_func(a, t): 695 return a[0] 696 return self.fit(const_func, plateau_range)[0] 697 elif method in ["avg", "average", "mean"]: 698 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) 699 return returnvalue 700 701 else: 702 raise Exception("Unsupported plateau method: " + method)
Extract a plateau value from a Corr object
Parameters
- plateau_range (list): list with two entries, indicating the first and the last timeslice of the plateau region.
- method (str): method to extract the plateau. 'fit' fits a constant to the plateau region 'avg', 'average' or 'mean' just average over the given timeslices.
- auto_gamma (bool): apply gamma_method with default parameters to the Corr. Defaults to None
704 def set_prange(self, prange): 705 """Sets the attribute prange of the Corr object.""" 706 if not len(prange) == 2: 707 raise Exception("prange must be a list or array with two values") 708 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): 709 raise Exception("Start and end point must be integers") 710 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): 711 raise Exception("Start and end point must define a range in the interval 0,T") 712 713 self.prange = prange 714 return
Sets the attribute prange of the Corr object.
716 def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None): 717 """Plots the correlator using the tag of the correlator as label if available. 718 719 Parameters 720 ---------- 721 x_range : list 722 list of two values, determining the range of the x-axis e.g. [4, 8] 723 comp : Corr or list of Corr 724 Correlator or list of correlators which are plotted for comparison. 725 The tags of these correlators are used as labels if available. 726 logscale : bool 727 Sets y-axis to logscale 728 plateau : Obs 729 Plateau value to be visualized in the figure 730 fit_res : Fit_result 731 Fit_result object to be visualized 732 ylabel : str 733 Label for the y-axis 734 save : str 735 path to file in which the figure should be saved 736 auto_gamma : bool 737 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. 738 hide_sigma : float 739 Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors. 740 references : list 741 List of floating point values that are displayed as horizontal lines for reference. 742 """ 743 if self.N != 1: 744 raise Exception("Correlator must be projected before plotting") 745 746 if auto_gamma: 747 self.gamma_method() 748 749 if x_range is None: 750 x_range = [0, self.T - 1] 751 752 fig = plt.figure() 753 ax1 = fig.add_subplot(111) 754 755 x, y, y_err = self.plottable() 756 if hide_sigma: 757 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1 758 else: 759 hide_from = None 760 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) 761 if logscale: 762 ax1.set_yscale('log') 763 else: 764 if y_range is None: 765 try: 766 y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)]) 767 y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)]) 768 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 769 except Exception: 770 pass 771 else: 772 ax1.set_ylim(y_range) 773 if comp: 774 if isinstance(comp, (Corr, list)): 775 for corr in comp if isinstance(comp, list) else [comp]: 776 if auto_gamma: 777 corr.gamma_method() 778 x, y, y_err = corr.plottable() 779 if hide_sigma: 780 hide_from = np.argmax((hide_sigma * np.array(y_err)) > np.abs(y)) - 1 781 else: 782 hide_from = None 783 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 784 else: 785 raise Exception("'comp' must be a correlator or a list of correlators.") 786 787 if plateau: 788 if isinstance(plateau, Obs): 789 if auto_gamma: 790 plateau.gamma_method() 791 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 792 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 793 else: 794 raise Exception("'plateau' must be an Obs") 795 796 if references: 797 if isinstance(references, list): 798 for ref in references: 799 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') 800 else: 801 raise Exception("'references' must be a list of floating pint values.") 802 803 if self.prange: 804 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 805 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 806 807 if fit_res: 808 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 809 ax1.plot(x_samples, 810 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 811 ls='-', marker=',', lw=2) 812 813 ax1.set_xlabel(r'$x_0 / a$') 814 if ylabel: 815 ax1.set_ylabel(ylabel) 816 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 817 818 handles, labels = ax1.get_legend_handles_labels() 819 if labels: 820 ax1.legend() 821 plt.draw() 822 823 if save: 824 if isinstance(save, str): 825 fig.savefig(save) 826 else: 827 raise Exception("'save' has to be a string.")
Plots the correlator using the tag of the correlator as label if available.
Parameters
- x_range (list): list of two values, determining the range of the x-axis e.g. [4, 8]
- comp (Corr or list of Corr): Correlator or list of correlators which are plotted for comparison. The tags of these correlators are used as labels if available.
- logscale (bool): Sets y-axis to logscale
- plateau (Obs): Plateau value to be visualized in the figure
- fit_res (Fit_result): Fit_result object to be visualized
- ylabel (str): Label for the y-axis
- save (str): path to file in which the figure should be saved
- auto_gamma (bool): Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
- hide_sigma (float): Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
- references (list): List of floating point values that are displayed as horizontal lines for reference.
829 def spaghetti_plot(self, logscale=True): 830 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 831 832 Parameters 833 ---------- 834 logscale : bool 835 Determines whether the scale of the y-axis is logarithmic or standard. 836 """ 837 if self.N != 1: 838 raise Exception("Correlator needs to be projected first.") 839 840 mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist])) 841 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 842 843 for name in mc_names: 844 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 845 846 fig = plt.figure() 847 ax = fig.add_subplot(111) 848 for dat in data: 849 ax.plot(x0_vals, dat, ls='-', marker='') 850 851 if logscale is True: 852 ax.set_yscale('log') 853 854 ax.set_xlabel(r'$x_0 / a$') 855 plt.title(name) 856 plt.draw()
Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
Parameters
- logscale (bool): Determines whether the scale of the y-axis is logarithmic or standard.
858 def dump(self, filename, datatype="json.gz", **kwargs): 859 """Dumps the Corr into a file of chosen type 860 Parameters 861 ---------- 862 filename : str 863 Name of the file to be saved. 864 datatype : str 865 Format of the exported file. Supported formats include 866 "json.gz" and "pickle" 867 path : str 868 specifies a custom path for the file (default '.') 869 """ 870 if datatype == "json.gz": 871 from .input.json import dump_to_json 872 if 'path' in kwargs: 873 file_name = kwargs.get('path') + '/' + filename 874 else: 875 file_name = filename 876 dump_to_json(self, file_name) 877 elif datatype == "pickle": 878 dump_object(self, filename, **kwargs) 879 else: 880 raise Exception("Unknown datatype " + str(datatype))
Dumps the Corr into a file of chosen type
Parameters
- filename (str): Name of the file to be saved.
- datatype (str): Format of the exported file. Supported formats include "json.gz" and "pickle"
- path (str): specifies a custom path for the file (default '.')
1136 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): 1137 r''' Project large correlation matrix to lowest states 1138 1139 This method can be used to reduce the size of an (N x N) correlation matrix 1140 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise 1141 is still small. 1142 1143 Parameters 1144 ---------- 1145 Ntrunc: int 1146 Rank of the target matrix. 1147 tproj: int 1148 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. 1149 The default value is 3. 1150 t0proj: int 1151 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly 1152 discouraged for O(a) improved theories, since the correctness of the procedure 1153 cannot be granted in this case. The default value is 2. 1154 basematrix : Corr 1155 Correlation matrix that is used to determine the eigenvectors of the 1156 lowest states based on a GEVP. basematrix is taken to be the Corr itself if 1157 is is not specified. 1158 1159 Notes 1160 ----- 1161 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving 1162 the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$ 1163 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the 1164 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via 1165 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large 1166 correlation matrix and to remove some noise that is added by irrelevant operators. 1167 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated 1168 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. 1169 ''' 1170 1171 if self.N == 1: 1172 raise Exception('Method cannot be applied to one-dimensional correlators.') 1173 if basematrix is None: 1174 basematrix = self 1175 if Ntrunc >= basematrix.N: 1176 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) 1177 if basematrix.N != self.N: 1178 raise Exception('basematrix and targetmatrix have to be of the same size.') 1179 1180 evecs = [] 1181 for i in range(Ntrunc): 1182 evecs.append(basematrix.GEVP(t0proj, tproj, state=i, sorted_list=None)) 1183 1184 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) 1185 rmat = [] 1186 for t in range(basematrix.T): 1187 for i in range(Ntrunc): 1188 for j in range(Ntrunc): 1189 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] 1190 rmat.append(np.copy(tmpmat)) 1191 1192 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] 1193 return Corr(newcontent)
Project large correlation matrix to lowest states
This method can be used to reduce the size of an (N x N) correlation matrix to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise is still small.
Parameters
- Ntrunc (int): Rank of the target matrix.
- tproj (int): Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. The default value is 3.
- t0proj (int): Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly discouraged for O(a) improved theories, since the correctness of the procedure cannot be granted in this case. The default value is 2.
- basematrix (Corr): Correlation matrix that is used to determine the eigenvectors of the lowest states based on a GEVP. basematrix is taken to be the Corr itself if is is not specified.
Notes
We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving 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}$ and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large correlation matrix and to remove some noise that is added by irrelevant operators. This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.