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