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