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