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