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, title=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 title : string 761 Optional title of the figure. 762 """ 763 if self.N != 1: 764 raise Exception("Correlator must be projected before plotting") 765 766 if auto_gamma: 767 self.gamma_method() 768 769 if x_range is None: 770 x_range = [0, self.T - 1] 771 772 fig = plt.figure() 773 ax1 = fig.add_subplot(111) 774 775 x, y, y_err = self.plottable() 776 if hide_sigma: 777 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 778 else: 779 hide_from = None 780 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) 781 if logscale: 782 ax1.set_yscale('log') 783 else: 784 if y_range is None: 785 try: 786 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)]) 787 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)]) 788 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 789 except Exception: 790 pass 791 else: 792 ax1.set_ylim(y_range) 793 if comp: 794 if isinstance(comp, (Corr, list)): 795 for corr in comp if isinstance(comp, list) else [comp]: 796 if auto_gamma: 797 corr.gamma_method() 798 x, y, y_err = corr.plottable() 799 if hide_sigma: 800 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 801 else: 802 hide_from = None 803 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 804 else: 805 raise Exception("'comp' must be a correlator or a list of correlators.") 806 807 if plateau: 808 if isinstance(plateau, Obs): 809 if auto_gamma: 810 plateau.gamma_method() 811 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 812 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 813 else: 814 raise Exception("'plateau' must be an Obs") 815 816 if references: 817 if isinstance(references, list): 818 for ref in references: 819 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') 820 else: 821 raise Exception("'references' must be a list of floating pint values.") 822 823 if self.prange: 824 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 825 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 826 827 if fit_res: 828 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 829 ax1.plot(x_samples, 830 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 831 ls='-', marker=',', lw=2) 832 833 ax1.set_xlabel(r'$x_0 / a$') 834 if ylabel: 835 ax1.set_ylabel(ylabel) 836 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 837 838 handles, labels = ax1.get_legend_handles_labels() 839 if labels: 840 ax1.legend() 841 842 if title: 843 plt.title(title) 844 845 plt.draw() 846 847 if save: 848 if isinstance(save, str): 849 fig.savefig(save, bbox_inches='tight') 850 else: 851 raise Exception("'save' has to be a string.") 852 853 def spaghetti_plot(self, logscale=True): 854 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 855 856 Parameters 857 ---------- 858 logscale : bool 859 Determines whether the scale of the y-axis is logarithmic or standard. 860 """ 861 if self.N != 1: 862 raise Exception("Correlator needs to be projected first.") 863 864 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])) 865 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 866 867 for name in mc_names: 868 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 869 870 fig = plt.figure() 871 ax = fig.add_subplot(111) 872 for dat in data: 873 ax.plot(x0_vals, dat, ls='-', marker='') 874 875 if logscale is True: 876 ax.set_yscale('log') 877 878 ax.set_xlabel(r'$x_0 / a$') 879 plt.title(name) 880 plt.draw() 881 882 def dump(self, filename, datatype="json.gz", **kwargs): 883 """Dumps the Corr into a file of chosen type 884 Parameters 885 ---------- 886 filename : str 887 Name of the file to be saved. 888 datatype : str 889 Format of the exported file. Supported formats include 890 "json.gz" and "pickle" 891 path : str 892 specifies a custom path for the file (default '.') 893 """ 894 if datatype == "json.gz": 895 from .input.json import dump_to_json 896 if 'path' in kwargs: 897 file_name = kwargs.get('path') + '/' + filename 898 else: 899 file_name = filename 900 dump_to_json(self, file_name) 901 elif datatype == "pickle": 902 dump_object(self, filename, **kwargs) 903 else: 904 raise Exception("Unknown datatype " + str(datatype)) 905 906 def print(self, print_range=None): 907 print(self.__repr__(print_range)) 908 909 def __repr__(self, print_range=None): 910 if print_range is None: 911 print_range = [0, None] 912 913 content_string = "" 914 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 915 916 if self.tag is not None: 917 content_string += "Description: " + self.tag + "\n" 918 if self.N != 1: 919 return content_string 920 921 if print_range[1]: 922 print_range[1] += 1 923 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' 924 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): 925 if sub_corr is None: 926 content_string += str(i + print_range[0]) + '\n' 927 else: 928 content_string += str(i + print_range[0]) 929 for element in sub_corr: 930 content_string += '\t' + ' ' * int(element >= 0) + str(element) 931 content_string += '\n' 932 return content_string 933 934 def __str__(self): 935 return self.__repr__() 936 937 # We define the basic operations, that can be performed with correlators. 938 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. 939 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. 940 # One could try and tell Obs to check if the y in __mul__ is a Corr and 941 942 def __add__(self, y): 943 if isinstance(y, Corr): 944 if ((self.N != y.N) or (self.T != y.T)): 945 raise Exception("Addition of Corrs with different shape") 946 newcontent = [] 947 for t in range(self.T): 948 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): 949 newcontent.append(None) 950 else: 951 newcontent.append(self.content[t] + y.content[t]) 952 return Corr(newcontent) 953 954 elif isinstance(y, (Obs, int, float, CObs)): 955 newcontent = [] 956 for t in range(self.T): 957 if _check_for_none(self, self.content[t]): 958 newcontent.append(None) 959 else: 960 newcontent.append(self.content[t] + y) 961 return Corr(newcontent, prange=self.prange) 962 elif isinstance(y, np.ndarray): 963 if y.shape == (self.T,): 964 return Corr(list((np.array(self.content).T + y).T)) 965 else: 966 raise ValueError("operands could not be broadcast together") 967 else: 968 raise TypeError("Corr + wrong type") 969 970 def __mul__(self, y): 971 if isinstance(y, Corr): 972 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 973 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 974 newcontent = [] 975 for t in range(self.T): 976 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): 977 newcontent.append(None) 978 else: 979 newcontent.append(self.content[t] * y.content[t]) 980 return Corr(newcontent) 981 982 elif isinstance(y, (Obs, int, float, CObs)): 983 newcontent = [] 984 for t in range(self.T): 985 if _check_for_none(self, self.content[t]): 986 newcontent.append(None) 987 else: 988 newcontent.append(self.content[t] * y) 989 return Corr(newcontent, prange=self.prange) 990 elif isinstance(y, np.ndarray): 991 if y.shape == (self.T,): 992 return Corr(list((np.array(self.content).T * y).T)) 993 else: 994 raise ValueError("operands could not be broadcast together") 995 else: 996 raise TypeError("Corr * wrong type") 997 998 def __truediv__(self, y): 999 if isinstance(y, Corr): 1000 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 1001 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 1002 newcontent = [] 1003 for t in range(self.T): 1004 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): 1005 newcontent.append(None) 1006 else: 1007 newcontent.append(self.content[t] / y.content[t]) 1008 for t in range(self.T): 1009 if _check_for_none(self, newcontent[t]): 1010 continue 1011 if np.isnan(np.sum(newcontent[t]).value): 1012 newcontent[t] = None 1013 1014 if all([item is None for item in newcontent]): 1015 raise Exception("Division returns completely undefined correlator") 1016 return Corr(newcontent) 1017 1018 elif isinstance(y, (Obs, CObs)): 1019 if isinstance(y, Obs): 1020 if y.value == 0: 1021 raise Exception('Division by zero will return undefined correlator') 1022 if isinstance(y, CObs): 1023 if y.is_zero(): 1024 raise Exception('Division by zero will return undefined correlator') 1025 1026 newcontent = [] 1027 for t in range(self.T): 1028 if _check_for_none(self, self.content[t]): 1029 newcontent.append(None) 1030 else: 1031 newcontent.append(self.content[t] / y) 1032 return Corr(newcontent, prange=self.prange) 1033 1034 elif isinstance(y, (int, float)): 1035 if y == 0: 1036 raise Exception('Division by zero will return undefined correlator') 1037 newcontent = [] 1038 for t in range(self.T): 1039 if _check_for_none(self, self.content[t]): 1040 newcontent.append(None) 1041 else: 1042 newcontent.append(self.content[t] / y) 1043 return Corr(newcontent, prange=self.prange) 1044 elif isinstance(y, np.ndarray): 1045 if y.shape == (self.T,): 1046 return Corr(list((np.array(self.content).T / y).T)) 1047 else: 1048 raise ValueError("operands could not be broadcast together") 1049 else: 1050 raise TypeError('Corr / wrong type') 1051 1052 def __neg__(self): 1053 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] 1054 return Corr(newcontent, prange=self.prange) 1055 1056 def __sub__(self, y): 1057 return self + (-y) 1058 1059 def __pow__(self, y): 1060 if isinstance(y, (Obs, int, float, CObs)): 1061 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] 1062 return Corr(newcontent, prange=self.prange) 1063 else: 1064 raise TypeError('Type of exponent not supported') 1065 1066 def __abs__(self): 1067 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] 1068 return Corr(newcontent, prange=self.prange) 1069 1070 # The numpy functions: 1071 def sqrt(self): 1072 return self ** 0.5 1073 1074 def log(self): 1075 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] 1076 return Corr(newcontent, prange=self.prange) 1077 1078 def exp(self): 1079 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] 1080 return Corr(newcontent, prange=self.prange) 1081 1082 def _apply_func_to_corr(self, func): 1083 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] 1084 for t in range(self.T): 1085 if _check_for_none(self, newcontent[t]): 1086 continue 1087 if np.isnan(np.sum(newcontent[t]).value): 1088 newcontent[t] = None 1089 if all([item is None for item in newcontent]): 1090 raise Exception('Operation returns undefined correlator') 1091 return Corr(newcontent) 1092 1093 def sin(self): 1094 return self._apply_func_to_corr(np.sin) 1095 1096 def cos(self): 1097 return self._apply_func_to_corr(np.cos) 1098 1099 def tan(self): 1100 return self._apply_func_to_corr(np.tan) 1101 1102 def sinh(self): 1103 return self._apply_func_to_corr(np.sinh) 1104 1105 def cosh(self): 1106 return self._apply_func_to_corr(np.cosh) 1107 1108 def tanh(self): 1109 return self._apply_func_to_corr(np.tanh) 1110 1111 def arcsin(self): 1112 return self._apply_func_to_corr(np.arcsin) 1113 1114 def arccos(self): 1115 return self._apply_func_to_corr(np.arccos) 1116 1117 def arctan(self): 1118 return self._apply_func_to_corr(np.arctan) 1119 1120 def arcsinh(self): 1121 return self._apply_func_to_corr(np.arcsinh) 1122 1123 def arccosh(self): 1124 return self._apply_func_to_corr(np.arccosh) 1125 1126 def arctanh(self): 1127 return self._apply_func_to_corr(np.arctanh) 1128 1129 # Right hand side operations (require tweak in main module to work) 1130 def __radd__(self, y): 1131 return self + y 1132 1133 def __rsub__(self, y): 1134 return -self + y 1135 1136 def __rmul__(self, y): 1137 return self * y 1138 1139 def __rtruediv__(self, y): 1140 return (self / y) ** (-1) 1141 1142 @property 1143 def real(self): 1144 def return_real(obs_OR_cobs): 1145 if isinstance(obs_OR_cobs, CObs): 1146 return obs_OR_cobs.real 1147 else: 1148 return obs_OR_cobs 1149 1150 return self._apply_func_to_corr(return_real) 1151 1152 @property 1153 def imag(self): 1154 def return_imag(obs_OR_cobs): 1155 if isinstance(obs_OR_cobs, CObs): 1156 return obs_OR_cobs.imag 1157 else: 1158 return obs_OR_cobs * 0 # So it stays the right type 1159 1160 return self._apply_func_to_corr(return_imag) 1161 1162 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): 1163 r''' Project large correlation matrix to lowest states 1164 1165 This method can be used to reduce the size of an (N x N) correlation matrix 1166 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise 1167 is still small. 1168 1169 Parameters 1170 ---------- 1171 Ntrunc: int 1172 Rank of the target matrix. 1173 tproj: int 1174 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. 1175 The default value is 3. 1176 t0proj: int 1177 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly 1178 discouraged for O(a) improved theories, since the correctness of the procedure 1179 cannot be granted in this case. The default value is 2. 1180 basematrix : Corr 1181 Correlation matrix that is used to determine the eigenvectors of the 1182 lowest states based on a GEVP. basematrix is taken to be the Corr itself if 1183 is is not specified. 1184 1185 Notes 1186 ----- 1187 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving 1188 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}$ 1189 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the 1190 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via 1191 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large 1192 correlation matrix and to remove some noise that is added by irrelevant operators. 1193 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated 1194 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. 1195 ''' 1196 1197 if self.N == 1: 1198 raise Exception('Method cannot be applied to one-dimensional correlators.') 1199 if basematrix is None: 1200 basematrix = self 1201 if Ntrunc >= basematrix.N: 1202 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) 1203 if basematrix.N != self.N: 1204 raise Exception('basematrix and targetmatrix have to be of the same size.') 1205 1206 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] 1207 1208 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) 1209 rmat = [] 1210 for t in range(basematrix.T): 1211 for i in range(Ntrunc): 1212 for j in range(Ntrunc): 1213 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] 1214 rmat.append(np.copy(tmpmat)) 1215 1216 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] 1217 return Corr(newcontent) 1218 1219 1220def _sort_vectors(vec_set, ts): 1221 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" 1222 reference_sorting = np.array(vec_set[ts]) 1223 N = reference_sorting.shape[0] 1224 sorted_vec_set = [] 1225 for t in range(len(vec_set)): 1226 if vec_set[t] is None: 1227 sorted_vec_set.append(None) 1228 elif not t == ts: 1229 perms = [list(o) for o in permutations([i for i in range(N)], N)] 1230 best_score = 0 1231 for perm in perms: 1232 current_score = 1 1233 for k in range(N): 1234 new_sorting = reference_sorting.copy() 1235 new_sorting[perm[k], :] = vec_set[t][k] 1236 current_score *= abs(np.linalg.det(new_sorting)) 1237 if current_score > best_score: 1238 best_score = current_score 1239 best_perm = perm 1240 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) 1241 else: 1242 sorted_vec_set.append(vec_set[t]) 1243 1244 return sorted_vec_set 1245 1246 1247def _check_for_none(corr, entry): 1248 """Checks if entry for correlator corr is None""" 1249 return len(list(filter(None, np.asarray(entry).flatten()))) < corr.N ** 2 1250 1251 1252def _GEVP_solver(Gt, G0): 1253 """Helper function for solving the GEVP and sorting the eigenvectors. 1254 1255 The helper function assumes that both provided matrices are symmetric and 1256 only processes the lower triangular part of both matrices. In case the matrices 1257 are not symmetric the upper triangular parts are effectively discarded.""" 1258 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, title=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 title : string 762 Optional title of the figure. 763 """ 764 if self.N != 1: 765 raise Exception("Correlator must be projected before plotting") 766 767 if auto_gamma: 768 self.gamma_method() 769 770 if x_range is None: 771 x_range = [0, self.T - 1] 772 773 fig = plt.figure() 774 ax1 = fig.add_subplot(111) 775 776 x, y, y_err = self.plottable() 777 if hide_sigma: 778 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 779 else: 780 hide_from = None 781 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) 782 if logscale: 783 ax1.set_yscale('log') 784 else: 785 if y_range is None: 786 try: 787 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)]) 788 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)]) 789 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 790 except Exception: 791 pass 792 else: 793 ax1.set_ylim(y_range) 794 if comp: 795 if isinstance(comp, (Corr, list)): 796 for corr in comp if isinstance(comp, list) else [comp]: 797 if auto_gamma: 798 corr.gamma_method() 799 x, y, y_err = corr.plottable() 800 if hide_sigma: 801 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 802 else: 803 hide_from = None 804 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 805 else: 806 raise Exception("'comp' must be a correlator or a list of correlators.") 807 808 if plateau: 809 if isinstance(plateau, Obs): 810 if auto_gamma: 811 plateau.gamma_method() 812 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 813 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 814 else: 815 raise Exception("'plateau' must be an Obs") 816 817 if references: 818 if isinstance(references, list): 819 for ref in references: 820 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') 821 else: 822 raise Exception("'references' must be a list of floating pint values.") 823 824 if self.prange: 825 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 826 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 827 828 if fit_res: 829 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 830 ax1.plot(x_samples, 831 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 832 ls='-', marker=',', lw=2) 833 834 ax1.set_xlabel(r'$x_0 / a$') 835 if ylabel: 836 ax1.set_ylabel(ylabel) 837 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 838 839 handles, labels = ax1.get_legend_handles_labels() 840 if labels: 841 ax1.legend() 842 843 if title: 844 plt.title(title) 845 846 plt.draw() 847 848 if save: 849 if isinstance(save, str): 850 fig.savefig(save, bbox_inches='tight') 851 else: 852 raise Exception("'save' has to be a string.") 853 854 def spaghetti_plot(self, logscale=True): 855 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 856 857 Parameters 858 ---------- 859 logscale : bool 860 Determines whether the scale of the y-axis is logarithmic or standard. 861 """ 862 if self.N != 1: 863 raise Exception("Correlator needs to be projected first.") 864 865 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])) 866 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 867 868 for name in mc_names: 869 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 870 871 fig = plt.figure() 872 ax = fig.add_subplot(111) 873 for dat in data: 874 ax.plot(x0_vals, dat, ls='-', marker='') 875 876 if logscale is True: 877 ax.set_yscale('log') 878 879 ax.set_xlabel(r'$x_0 / a$') 880 plt.title(name) 881 plt.draw() 882 883 def dump(self, filename, datatype="json.gz", **kwargs): 884 """Dumps the Corr into a file of chosen type 885 Parameters 886 ---------- 887 filename : str 888 Name of the file to be saved. 889 datatype : str 890 Format of the exported file. Supported formats include 891 "json.gz" and "pickle" 892 path : str 893 specifies a custom path for the file (default '.') 894 """ 895 if datatype == "json.gz": 896 from .input.json import dump_to_json 897 if 'path' in kwargs: 898 file_name = kwargs.get('path') + '/' + filename 899 else: 900 file_name = filename 901 dump_to_json(self, file_name) 902 elif datatype == "pickle": 903 dump_object(self, filename, **kwargs) 904 else: 905 raise Exception("Unknown datatype " + str(datatype)) 906 907 def print(self, print_range=None): 908 print(self.__repr__(print_range)) 909 910 def __repr__(self, print_range=None): 911 if print_range is None: 912 print_range = [0, None] 913 914 content_string = "" 915 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 916 917 if self.tag is not None: 918 content_string += "Description: " + self.tag + "\n" 919 if self.N != 1: 920 return content_string 921 922 if print_range[1]: 923 print_range[1] += 1 924 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' 925 for i, sub_corr in enumerate(self.content[print_range[0]:print_range[1]]): 926 if sub_corr is None: 927 content_string += str(i + print_range[0]) + '\n' 928 else: 929 content_string += str(i + print_range[0]) 930 for element in sub_corr: 931 content_string += '\t' + ' ' * int(element >= 0) + str(element) 932 content_string += '\n' 933 return content_string 934 935 def __str__(self): 936 return self.__repr__() 937 938 # We define the basic operations, that can be performed with correlators. 939 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. 940 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. 941 # One could try and tell Obs to check if the y in __mul__ is a Corr and 942 943 def __add__(self, y): 944 if isinstance(y, Corr): 945 if ((self.N != y.N) or (self.T != y.T)): 946 raise Exception("Addition of Corrs with different shape") 947 newcontent = [] 948 for t in range(self.T): 949 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): 950 newcontent.append(None) 951 else: 952 newcontent.append(self.content[t] + y.content[t]) 953 return Corr(newcontent) 954 955 elif isinstance(y, (Obs, int, float, CObs)): 956 newcontent = [] 957 for t in range(self.T): 958 if _check_for_none(self, self.content[t]): 959 newcontent.append(None) 960 else: 961 newcontent.append(self.content[t] + y) 962 return Corr(newcontent, prange=self.prange) 963 elif isinstance(y, np.ndarray): 964 if y.shape == (self.T,): 965 return Corr(list((np.array(self.content).T + y).T)) 966 else: 967 raise ValueError("operands could not be broadcast together") 968 else: 969 raise TypeError("Corr + wrong type") 970 971 def __mul__(self, y): 972 if isinstance(y, Corr): 973 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 974 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 975 newcontent = [] 976 for t in range(self.T): 977 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): 978 newcontent.append(None) 979 else: 980 newcontent.append(self.content[t] * y.content[t]) 981 return Corr(newcontent) 982 983 elif isinstance(y, (Obs, int, float, CObs)): 984 newcontent = [] 985 for t in range(self.T): 986 if _check_for_none(self, self.content[t]): 987 newcontent.append(None) 988 else: 989 newcontent.append(self.content[t] * y) 990 return Corr(newcontent, prange=self.prange) 991 elif isinstance(y, np.ndarray): 992 if y.shape == (self.T,): 993 return Corr(list((np.array(self.content).T * y).T)) 994 else: 995 raise ValueError("operands could not be broadcast together") 996 else: 997 raise TypeError("Corr * wrong type") 998 999 def __truediv__(self, y): 1000 if isinstance(y, Corr): 1001 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 1002 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 1003 newcontent = [] 1004 for t in range(self.T): 1005 if _check_for_none(self, self.content[t]) or _check_for_none(y, y.content[t]): 1006 newcontent.append(None) 1007 else: 1008 newcontent.append(self.content[t] / y.content[t]) 1009 for t in range(self.T): 1010 if _check_for_none(self, newcontent[t]): 1011 continue 1012 if np.isnan(np.sum(newcontent[t]).value): 1013 newcontent[t] = None 1014 1015 if all([item is None for item in newcontent]): 1016 raise Exception("Division returns completely undefined correlator") 1017 return Corr(newcontent) 1018 1019 elif isinstance(y, (Obs, CObs)): 1020 if isinstance(y, Obs): 1021 if y.value == 0: 1022 raise Exception('Division by zero will return undefined correlator') 1023 if isinstance(y, CObs): 1024 if y.is_zero(): 1025 raise Exception('Division by zero will return undefined correlator') 1026 1027 newcontent = [] 1028 for t in range(self.T): 1029 if _check_for_none(self, self.content[t]): 1030 newcontent.append(None) 1031 else: 1032 newcontent.append(self.content[t] / y) 1033 return Corr(newcontent, prange=self.prange) 1034 1035 elif isinstance(y, (int, float)): 1036 if y == 0: 1037 raise Exception('Division by zero will return undefined correlator') 1038 newcontent = [] 1039 for t in range(self.T): 1040 if _check_for_none(self, self.content[t]): 1041 newcontent.append(None) 1042 else: 1043 newcontent.append(self.content[t] / y) 1044 return Corr(newcontent, prange=self.prange) 1045 elif isinstance(y, np.ndarray): 1046 if y.shape == (self.T,): 1047 return Corr(list((np.array(self.content).T / y).T)) 1048 else: 1049 raise ValueError("operands could not be broadcast together") 1050 else: 1051 raise TypeError('Corr / wrong type') 1052 1053 def __neg__(self): 1054 newcontent = [None if _check_for_none(self, item) else -1. * item for item in self.content] 1055 return Corr(newcontent, prange=self.prange) 1056 1057 def __sub__(self, y): 1058 return self + (-y) 1059 1060 def __pow__(self, y): 1061 if isinstance(y, (Obs, int, float, CObs)): 1062 newcontent = [None if _check_for_none(self, item) else item**y for item in self.content] 1063 return Corr(newcontent, prange=self.prange) 1064 else: 1065 raise TypeError('Type of exponent not supported') 1066 1067 def __abs__(self): 1068 newcontent = [None if _check_for_none(self, item) else np.abs(item) for item in self.content] 1069 return Corr(newcontent, prange=self.prange) 1070 1071 # The numpy functions: 1072 def sqrt(self): 1073 return self ** 0.5 1074 1075 def log(self): 1076 newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content] 1077 return Corr(newcontent, prange=self.prange) 1078 1079 def exp(self): 1080 newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content] 1081 return Corr(newcontent, prange=self.prange) 1082 1083 def _apply_func_to_corr(self, func): 1084 newcontent = [None if _check_for_none(self, item) else func(item) for item in self.content] 1085 for t in range(self.T): 1086 if _check_for_none(self, newcontent[t]): 1087 continue 1088 if np.isnan(np.sum(newcontent[t]).value): 1089 newcontent[t] = None 1090 if all([item is None for item in newcontent]): 1091 raise Exception('Operation returns undefined correlator') 1092 return Corr(newcontent) 1093 1094 def sin(self): 1095 return self._apply_func_to_corr(np.sin) 1096 1097 def cos(self): 1098 return self._apply_func_to_corr(np.cos) 1099 1100 def tan(self): 1101 return self._apply_func_to_corr(np.tan) 1102 1103 def sinh(self): 1104 return self._apply_func_to_corr(np.sinh) 1105 1106 def cosh(self): 1107 return self._apply_func_to_corr(np.cosh) 1108 1109 def tanh(self): 1110 return self._apply_func_to_corr(np.tanh) 1111 1112 def arcsin(self): 1113 return self._apply_func_to_corr(np.arcsin) 1114 1115 def arccos(self): 1116 return self._apply_func_to_corr(np.arccos) 1117 1118 def arctan(self): 1119 return self._apply_func_to_corr(np.arctan) 1120 1121 def arcsinh(self): 1122 return self._apply_func_to_corr(np.arcsinh) 1123 1124 def arccosh(self): 1125 return self._apply_func_to_corr(np.arccosh) 1126 1127 def arctanh(self): 1128 return self._apply_func_to_corr(np.arctanh) 1129 1130 # Right hand side operations (require tweak in main module to work) 1131 def __radd__(self, y): 1132 return self + y 1133 1134 def __rsub__(self, y): 1135 return -self + y 1136 1137 def __rmul__(self, y): 1138 return self * y 1139 1140 def __rtruediv__(self, y): 1141 return (self / y) ** (-1) 1142 1143 @property 1144 def real(self): 1145 def return_real(obs_OR_cobs): 1146 if isinstance(obs_OR_cobs, CObs): 1147 return obs_OR_cobs.real 1148 else: 1149 return obs_OR_cobs 1150 1151 return self._apply_func_to_corr(return_real) 1152 1153 @property 1154 def imag(self): 1155 def return_imag(obs_OR_cobs): 1156 if isinstance(obs_OR_cobs, CObs): 1157 return obs_OR_cobs.imag 1158 else: 1159 return obs_OR_cobs * 0 # So it stays the right type 1160 1161 return self._apply_func_to_corr(return_imag) 1162 1163 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): 1164 r''' Project large correlation matrix to lowest states 1165 1166 This method can be used to reduce the size of an (N x N) correlation matrix 1167 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise 1168 is still small. 1169 1170 Parameters 1171 ---------- 1172 Ntrunc: int 1173 Rank of the target matrix. 1174 tproj: int 1175 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. 1176 The default value is 3. 1177 t0proj: int 1178 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly 1179 discouraged for O(a) improved theories, since the correctness of the procedure 1180 cannot be granted in this case. The default value is 2. 1181 basematrix : Corr 1182 Correlation matrix that is used to determine the eigenvectors of the 1183 lowest states based on a GEVP. basematrix is taken to be the Corr itself if 1184 is is not specified. 1185 1186 Notes 1187 ----- 1188 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving 1189 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}$ 1190 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the 1191 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via 1192 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large 1193 correlation matrix and to remove some noise that is added by irrelevant operators. 1194 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated 1195 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. 1196 ''' 1197 1198 if self.N == 1: 1199 raise Exception('Method cannot be applied to one-dimensional correlators.') 1200 if basematrix is None: 1201 basematrix = self 1202 if Ntrunc >= basematrix.N: 1203 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) 1204 if basematrix.N != self.N: 1205 raise Exception('basematrix and targetmatrix have to be of the same size.') 1206 1207 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] 1208 1209 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) 1210 rmat = [] 1211 for t in range(basematrix.T): 1212 for i in range(Ntrunc): 1213 for j in range(Ntrunc): 1214 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] 1215 rmat.append(np.copy(tmpmat)) 1216 1217 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] 1218 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, title=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 title : string 762 Optional title of the figure. 763 """ 764 if self.N != 1: 765 raise Exception("Correlator must be projected before plotting") 766 767 if auto_gamma: 768 self.gamma_method() 769 770 if x_range is None: 771 x_range = [0, self.T - 1] 772 773 fig = plt.figure() 774 ax1 = fig.add_subplot(111) 775 776 x, y, y_err = self.plottable() 777 if hide_sigma: 778 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 779 else: 780 hide_from = None 781 ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag) 782 if logscale: 783 ax1.set_yscale('log') 784 else: 785 if y_range is None: 786 try: 787 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)]) 788 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)]) 789 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 790 except Exception: 791 pass 792 else: 793 ax1.set_ylim(y_range) 794 if comp: 795 if isinstance(comp, (Corr, list)): 796 for corr in comp if isinstance(comp, list) else [comp]: 797 if auto_gamma: 798 corr.gamma_method() 799 x, y, y_err = corr.plottable() 800 if hide_sigma: 801 hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1 802 else: 803 hide_from = None 804 plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 805 else: 806 raise Exception("'comp' must be a correlator or a list of correlators.") 807 808 if plateau: 809 if isinstance(plateau, Obs): 810 if auto_gamma: 811 plateau.gamma_method() 812 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 813 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 814 else: 815 raise Exception("'plateau' must be an Obs") 816 817 if references: 818 if isinstance(references, list): 819 for ref in references: 820 ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--') 821 else: 822 raise Exception("'references' must be a list of floating pint values.") 823 824 if self.prange: 825 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 826 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 827 828 if fit_res: 829 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 830 ax1.plot(x_samples, 831 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 832 ls='-', marker=',', lw=2) 833 834 ax1.set_xlabel(r'$x_0 / a$') 835 if ylabel: 836 ax1.set_ylabel(ylabel) 837 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 838 839 handles, labels = ax1.get_legend_handles_labels() 840 if labels: 841 ax1.legend() 842 843 if title: 844 plt.title(title) 845 846 plt.draw() 847 848 if save: 849 if isinstance(save, str): 850 fig.savefig(save, bbox_inches='tight') 851 else: 852 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.
- title (string): Optional title of the figure.
854 def spaghetti_plot(self, logscale=True): 855 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 856 857 Parameters 858 ---------- 859 logscale : bool 860 Determines whether the scale of the y-axis is logarithmic or standard. 861 """ 862 if self.N != 1: 863 raise Exception("Correlator needs to be projected first.") 864 865 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])) 866 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 867 868 for name in mc_names: 869 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 870 871 fig = plt.figure() 872 ax = fig.add_subplot(111) 873 for dat in data: 874 ax.plot(x0_vals, dat, ls='-', marker='') 875 876 if logscale is True: 877 ax.set_yscale('log') 878 879 ax.set_xlabel(r'$x_0 / a$') 880 plt.title(name) 881 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.
883 def dump(self, filename, datatype="json.gz", **kwargs): 884 """Dumps the Corr into a file of chosen type 885 Parameters 886 ---------- 887 filename : str 888 Name of the file to be saved. 889 datatype : str 890 Format of the exported file. Supported formats include 891 "json.gz" and "pickle" 892 path : str 893 specifies a custom path for the file (default '.') 894 """ 895 if datatype == "json.gz": 896 from .input.json import dump_to_json 897 if 'path' in kwargs: 898 file_name = kwargs.get('path') + '/' + filename 899 else: 900 file_name = filename 901 dump_to_json(self, file_name) 902 elif datatype == "pickle": 903 dump_object(self, filename, **kwargs) 904 else: 905 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 '.')
1163 def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None): 1164 r''' Project large correlation matrix to lowest states 1165 1166 This method can be used to reduce the size of an (N x N) correlation matrix 1167 to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise 1168 is still small. 1169 1170 Parameters 1171 ---------- 1172 Ntrunc: int 1173 Rank of the target matrix. 1174 tproj: int 1175 Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method. 1176 The default value is 3. 1177 t0proj: int 1178 Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly 1179 discouraged for O(a) improved theories, since the correctness of the procedure 1180 cannot be granted in this case. The default value is 2. 1181 basematrix : Corr 1182 Correlation matrix that is used to determine the eigenvectors of the 1183 lowest states based on a GEVP. basematrix is taken to be the Corr itself if 1184 is is not specified. 1185 1186 Notes 1187 ----- 1188 We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving 1189 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}$ 1190 and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the 1191 resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via 1192 $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large 1193 correlation matrix and to remove some noise that is added by irrelevant operators. 1194 This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated 1195 bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$. 1196 ''' 1197 1198 if self.N == 1: 1199 raise Exception('Method cannot be applied to one-dimensional correlators.') 1200 if basematrix is None: 1201 basematrix = self 1202 if Ntrunc >= basematrix.N: 1203 raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N)) 1204 if basematrix.N != self.N: 1205 raise Exception('basematrix and targetmatrix have to be of the same size.') 1206 1207 evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc] 1208 1209 tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object) 1210 rmat = [] 1211 for t in range(basematrix.T): 1212 for i in range(Ntrunc): 1213 for j in range(Ntrunc): 1214 tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j] 1215 rmat.append(np.copy(tmpmat)) 1216 1217 newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)] 1218 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)$.