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 else: 277 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 301 return all_vecs 302 303 def Eigenvalue(self, t0, ts=None, state=0, sorted_list=None): 304 """Determines the eigenvalue of the GEVP by solving and projecting the correlator 305 306 Parameters 307 ---------- 308 t0 : int 309 The time t0 for G(t)v= lambda G(t_0)v 310 ts : int 311 fixed time G(t_s)v= lambda G(t_0)v if return_list=False 312 If return_list=True and sorting=Eigenvector it gives a reference point for the sorting method. 313 state : int 314 The state one is interested in ordered by energy. The lowest state is zero. 315 sorted_list : string 316 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. 317 "Eigenvalue" - The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice. 318 "Eigenvector" - Use the method described in arXiv:2004.10472 [hep-lat] to find the set of v(t) belonging to the state. 319 The reference state is identified by its eigenvalue at t=ts 320 """ 321 vec = self.GEVP(t0, ts=ts, state=state, sorted_list=sorted_list) 322 return self.projected(vec) 323 324 def Hankel(self, N, periodic=False): 325 """Constructs an NxN Hankel matrix 326 327 C(t) c(t+1) ... c(t+n-1) 328 C(t+1) c(t+2) ... c(t+n) 329 ................. 330 C(t+(n-1)) c(t+n) ... c(t+2(n-1)) 331 332 Parameters 333 ---------- 334 N : int 335 Dimension of the Hankel matrix 336 periodic : bool, optional 337 determines whether the matrix is extended periodically 338 """ 339 340 if self.N != 1: 341 raise Exception("Multi-operator Prony not implemented!") 342 343 array = np.empty([N, N], dtype="object") 344 new_content = [] 345 for t in range(self.T): 346 new_content.append(array.copy()) 347 348 def wrap(i): 349 while i >= self.T: 350 i -= self.T 351 return i 352 353 for t in range(self.T): 354 for i in range(N): 355 for j in range(N): 356 if periodic: 357 new_content[t][i, j] = self.content[wrap(t + i + j)][0] 358 elif (t + i + j) >= self.T: 359 new_content[t] = None 360 else: 361 new_content[t][i, j] = self.content[t + i + j][0] 362 363 return Corr(new_content) 364 365 def roll(self, dt): 366 """Periodically shift the correlator by dt timeslices 367 368 Parameters 369 ---------- 370 dt : int 371 number of timeslices 372 """ 373 return Corr(list(np.roll(np.array(self.content, dtype=object), dt))) 374 375 def reverse(self): 376 """Reverse the time ordering of the Corr""" 377 return Corr(self.content[:: -1]) 378 379 def thin(self, spacing=2, offset=0): 380 """Thin out a correlator to suppress correlations 381 382 Parameters 383 ---------- 384 spacing : int 385 Keep only every 'spacing'th entry of the correlator 386 offset : int 387 Offset the equal spacing 388 """ 389 new_content = [] 390 for t in range(self.T): 391 if (offset + t) % spacing != 0: 392 new_content.append(None) 393 else: 394 new_content.append(self.content[t]) 395 return Corr(new_content) 396 397 def correlate(self, partner): 398 """Correlate the correlator with another correlator or Obs 399 400 Parameters 401 ---------- 402 partner : Obs or Corr 403 partner to correlate the correlator with. 404 Can either be an Obs which is correlated with all entries of the 405 correlator or a Corr of same length. 406 """ 407 new_content = [] 408 for x0, t_slice in enumerate(self.content): 409 if t_slice is None: 410 new_content.append(None) 411 else: 412 if isinstance(partner, Corr): 413 if partner.content[x0] is None: 414 new_content.append(None) 415 else: 416 new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice])) 417 elif isinstance(partner, Obs): # Should this include CObs? 418 new_content.append(np.array([correlate(o, partner) for o in t_slice])) 419 else: 420 raise Exception("Can only correlate with an Obs or a Corr.") 421 422 return Corr(new_content) 423 424 def reweight(self, weight, **kwargs): 425 """Reweight the correlator. 426 427 Parameters 428 ---------- 429 weight : Obs 430 Reweighting factor. An Observable that has to be defined on a superset of the 431 configurations in obs[i].idl for all i. 432 all_configs : bool 433 if True, the reweighted observables are normalized by the average of 434 the reweighting factor on all configurations in weight.idl and not 435 on the configurations in obs[i].idl. 436 """ 437 new_content = [] 438 for t_slice in self.content: 439 if t_slice is None: 440 new_content.append(None) 441 else: 442 new_content.append(np.array(reweight(weight, t_slice, **kwargs))) 443 return Corr(new_content) 444 445 def T_symmetry(self, partner, parity=+1): 446 """Return the time symmetry average of the correlator and its partner 447 448 Parameters 449 ---------- 450 partner : Corr 451 Time symmetry partner of the Corr 452 partity : int 453 Parity quantum number of the correlator, can be +1 or -1 454 """ 455 if not isinstance(partner, Corr): 456 raise Exception("T partner has to be a Corr object.") 457 if parity not in [+1, -1]: 458 raise Exception("Parity has to be +1 or -1.") 459 T_partner = parity * partner.reverse() 460 461 t_slices = [] 462 test = (self - T_partner) 463 test.gamma_method() 464 for x0, t_slice in enumerate(test.content): 465 if t_slice is not None: 466 if not t_slice[0].is_zero_within_error(5): 467 t_slices.append(x0) 468 if t_slices: 469 warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning) 470 471 return (self + T_partner) / 2 472 473 def deriv(self, variant="symmetric"): 474 """Return the first derivative of the correlator with respect to x0. 475 476 Parameters 477 ---------- 478 variant : str 479 decides which definition of the finite differences derivative is used. 480 Available choice: symmetric, forward, backward, improved, default: symmetric 481 """ 482 if variant == "symmetric": 483 newcontent = [] 484 for t in range(1, self.T - 1): 485 if (self.content[t - 1] is None) or (self.content[t + 1] is None): 486 newcontent.append(None) 487 else: 488 newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1])) 489 if(all([x is None for x in newcontent])): 490 raise Exception('Derivative is undefined at all timeslices') 491 return Corr(newcontent, padding=[1, 1]) 492 elif variant == "forward": 493 newcontent = [] 494 for t in range(self.T - 1): 495 if (self.content[t] is None) or (self.content[t + 1] is None): 496 newcontent.append(None) 497 else: 498 newcontent.append(self.content[t + 1] - self.content[t]) 499 if(all([x is None for x in newcontent])): 500 raise Exception("Derivative is undefined at all timeslices") 501 return Corr(newcontent, padding=[0, 1]) 502 elif variant == "backward": 503 newcontent = [] 504 for t in range(1, self.T): 505 if (self.content[t - 1] is None) or (self.content[t] is None): 506 newcontent.append(None) 507 else: 508 newcontent.append(self.content[t] - self.content[t - 1]) 509 if(all([x is None for x in newcontent])): 510 raise Exception("Derivative is undefined at all timeslices") 511 return Corr(newcontent, padding=[1, 0]) 512 elif variant == "improved": 513 newcontent = [] 514 for t in range(2, self.T - 2): 515 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): 516 newcontent.append(None) 517 else: 518 newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2])) 519 if(all([x is None for x in newcontent])): 520 raise Exception('Derivative is undefined at all timeslices') 521 return Corr(newcontent, padding=[2, 2]) 522 else: 523 raise Exception("Unknown variant.") 524 525 def second_deriv(self, variant="symmetric"): 526 """Return the second derivative of the correlator with respect to x0. 527 528 Parameters 529 ---------- 530 variant : str 531 decides which definition of the finite differences derivative is used. 532 Available choice: symmetric, improved, default: symmetric 533 """ 534 if variant == "symmetric": 535 newcontent = [] 536 for t in range(1, self.T - 1): 537 if (self.content[t - 1] is None) or (self.content[t + 1] is None): 538 newcontent.append(None) 539 else: 540 newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1])) 541 if(all([x is None for x in newcontent])): 542 raise Exception("Derivative is undefined at all timeslices") 543 return Corr(newcontent, padding=[1, 1]) 544 elif variant == "improved": 545 newcontent = [] 546 for t in range(2, self.T - 2): 547 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): 548 newcontent.append(None) 549 else: 550 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])) 551 if(all([x is None for x in newcontent])): 552 raise Exception("Derivative is undefined at all timeslices") 553 return Corr(newcontent, padding=[2, 2]) 554 else: 555 raise Exception("Unknown variant.") 556 557 def m_eff(self, variant='log', guess=1.0): 558 """Returns the effective mass of the correlator as correlator object 559 560 Parameters 561 ---------- 562 variant : str 563 log : uses the standard effective mass log(C(t) / C(t+1)) 564 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. 565 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. 566 See, e.g., arXiv:1205.5380 567 arccosh : Uses the explicit form of the symmetrized correlator (not recommended) 568 guess : float 569 guess for the root finder, only relevant for the root variant 570 """ 571 if self.N != 1: 572 raise Exception('Correlator must be projected before getting m_eff') 573 if variant == 'log': 574 newcontent = [] 575 for t in range(self.T - 1): 576 if (self.content[t] is None) or (self.content[t + 1] is None): 577 newcontent.append(None) 578 else: 579 newcontent.append(self.content[t] / self.content[t + 1]) 580 if(all([x is None for x in newcontent])): 581 raise Exception('m_eff is undefined at all timeslices') 582 583 return np.log(Corr(newcontent, padding=[0, 1])) 584 585 elif variant in ['periodic', 'cosh', 'sinh']: 586 if variant in ['periodic', 'cosh']: 587 func = anp.cosh 588 else: 589 func = anp.sinh 590 591 def root_function(x, d): 592 return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d 593 594 newcontent = [] 595 for t in range(self.T - 1): 596 if (self.content[t] is None) or (self.content[t + 1] is None): 597 newcontent.append(None) 598 # Fill the two timeslices in the middle of the lattice with their predecessors 599 elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]: 600 newcontent.append(newcontent[-1]) 601 else: 602 newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess))) 603 if(all([x is None for x in newcontent])): 604 raise Exception('m_eff is undefined at all timeslices') 605 606 return Corr(newcontent, padding=[0, 1]) 607 608 elif variant == 'arccosh': 609 newcontent = [] 610 for t in range(1, self.T - 1): 611 if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None): 612 newcontent.append(None) 613 else: 614 newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t])) 615 if(all([x is None for x in newcontent])): 616 raise Exception("m_eff is undefined at all timeslices") 617 return np.arccosh(Corr(newcontent, padding=[1, 1])) 618 619 else: 620 raise Exception('Unknown variant.') 621 622 def fit(self, function, fitrange=None, silent=False, **kwargs): 623 r'''Fits function to the data 624 625 Parameters 626 ---------- 627 function : obj 628 function to fit to the data. See fits.least_squares for details. 629 fitrange : list 630 Two element list containing the timeslices on which the fit is supposed to start and stop. 631 Caution: This range is inclusive as opposed to standard python indexing. 632 `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6. 633 If not specified, self.prange or all timeslices are used. 634 silent : bool 635 Decides whether output is printed to the standard output. 636 ''' 637 if self.N != 1: 638 raise Exception("Correlator must be projected before fitting") 639 640 if fitrange is None: 641 if self.prange: 642 fitrange = self.prange 643 else: 644 fitrange = [0, self.T - 1] 645 else: 646 if not isinstance(fitrange, list): 647 raise Exception("fitrange has to be a list with two elements") 648 if len(fitrange) != 2: 649 raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]") 650 651 xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] 652 ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None] 653 result = least_squares(xs, ys, function, silent=silent, **kwargs) 654 return result 655 656 def plateau(self, plateau_range=None, method="fit", auto_gamma=False): 657 """ Extract a plateau value from a Corr object 658 659 Parameters 660 ---------- 661 plateau_range : list 662 list with two entries, indicating the first and the last timeslice 663 of the plateau region. 664 method : str 665 method to extract the plateau. 666 'fit' fits a constant to the plateau region 667 'avg', 'average' or 'mean' just average over the given timeslices. 668 auto_gamma : bool 669 apply gamma_method with default parameters to the Corr. Defaults to None 670 """ 671 if not plateau_range: 672 if self.prange: 673 plateau_range = self.prange 674 else: 675 raise Exception("no plateau range provided") 676 if self.N != 1: 677 raise Exception("Correlator must be projected before getting a plateau.") 678 if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])): 679 raise Exception("plateau is undefined at all timeslices in plateaurange.") 680 if auto_gamma: 681 self.gamma_method() 682 if method == "fit": 683 def const_func(a, t): 684 return a[0] 685 return self.fit(const_func, plateau_range)[0] 686 elif method in ["avg", "average", "mean"]: 687 returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None]) 688 return returnvalue 689 690 else: 691 raise Exception("Unsupported plateau method: " + method) 692 693 def set_prange(self, prange): 694 """Sets the attribute prange of the Corr object.""" 695 if not len(prange) == 2: 696 raise Exception("prange must be a list or array with two values") 697 if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))): 698 raise Exception("Start and end point must be integers") 699 if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]): 700 raise Exception("Start and end point must define a range in the interval 0,T") 701 702 self.prange = prange 703 return 704 705 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): 706 """Plots the correlator using the tag of the correlator as label if available. 707 708 Parameters 709 ---------- 710 x_range : list 711 list of two values, determining the range of the x-axis e.g. [4, 8] 712 comp : Corr or list of Corr 713 Correlator or list of correlators which are plotted for comparison. 714 The tags of these correlators are used as labels if available. 715 logscale : bool 716 Sets y-axis to logscale 717 plateau : Obs 718 Plateau value to be visualized in the figure 719 fit_res : Fit_result 720 Fit_result object to be visualized 721 ylabel : str 722 Label for the y-axis 723 save : str 724 path to file in which the figure should be saved 725 auto_gamma : bool 726 Apply the gamma method with standard parameters to all correlators and plateau values before plotting. 727 """ 728 if self.N != 1: 729 raise Exception("Correlator must be projected before plotting") 730 731 if auto_gamma: 732 self.gamma_method() 733 734 if x_range is None: 735 x_range = [0, self.T - 1] 736 737 fig = plt.figure() 738 ax1 = fig.add_subplot(111) 739 740 x, y, y_err = self.plottable() 741 ax1.errorbar(x, y, y_err, label=self.tag) 742 if logscale: 743 ax1.set_yscale('log') 744 else: 745 if y_range is None: 746 try: 747 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)]) 748 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)]) 749 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 750 except Exception: 751 pass 752 else: 753 ax1.set_ylim(y_range) 754 if comp: 755 if isinstance(comp, (Corr, list)): 756 for corr in comp if isinstance(comp, list) else [comp]: 757 if auto_gamma: 758 corr.gamma_method() 759 x, y, y_err = corr.plottable() 760 plt.errorbar(x, y, y_err, label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 761 else: 762 raise Exception("'comp' must be a correlator or a list of correlators.") 763 764 if plateau: 765 if isinstance(plateau, Obs): 766 if auto_gamma: 767 plateau.gamma_method() 768 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 769 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 770 else: 771 raise Exception("'plateau' must be an Obs") 772 if self.prange: 773 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 774 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 775 776 if fit_res: 777 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 778 ax1.plot(x_samples, 779 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 780 ls='-', marker=',', lw=2) 781 782 ax1.set_xlabel(r'$x_0 / a$') 783 if ylabel: 784 ax1.set_ylabel(ylabel) 785 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 786 787 handles, labels = ax1.get_legend_handles_labels() 788 if labels: 789 ax1.legend() 790 plt.draw() 791 792 if save: 793 if isinstance(save, str): 794 fig.savefig(save) 795 else: 796 raise Exception("'save' has to be a string.") 797 798 def spaghetti_plot(self, logscale=True): 799 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 800 801 Parameters 802 ---------- 803 logscale : bool 804 Determines whether the scale of the y-axis is logarithmic or standard. 805 """ 806 if self.N != 1: 807 raise Exception("Correlator needs to be projected first.") 808 809 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])) 810 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 811 812 for name in mc_names: 813 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 814 815 fig = plt.figure() 816 ax = fig.add_subplot(111) 817 for dat in data: 818 ax.plot(x0_vals, dat, ls='-', marker='') 819 820 if logscale is True: 821 ax.set_yscale('log') 822 823 ax.set_xlabel(r'$x_0 / a$') 824 plt.title(name) 825 plt.draw() 826 827 def dump(self, filename, datatype="json.gz", **kwargs): 828 """Dumps the Corr into a file of chosen type 829 Parameters 830 ---------- 831 filename : str 832 Name of the file to be saved. 833 datatype : str 834 Format of the exported file. Supported formats include 835 "json.gz" and "pickle" 836 path : str 837 specifies a custom path for the file (default '.') 838 """ 839 if datatype == "json.gz": 840 from .input.json import dump_to_json 841 if 'path' in kwargs: 842 file_name = kwargs.get('path') + '/' + filename 843 else: 844 file_name = filename 845 dump_to_json(self, file_name) 846 elif datatype == "pickle": 847 dump_object(self, filename, **kwargs) 848 else: 849 raise Exception("Unknown datatype " + str(datatype)) 850 851 def print(self, range=[0, None]): 852 print(self.__repr__(range)) 853 854 def __repr__(self, range=[0, None]): 855 content_string = "" 856 857 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 858 859 if self.tag is not None: 860 content_string += "Description: " + self.tag + "\n" 861 if self.N != 1: 862 return content_string 863 864 if range[1]: 865 range[1] += 1 866 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' 867 for i, sub_corr in enumerate(self.content[range[0]:range[1]]): 868 if sub_corr is None: 869 content_string += str(i + range[0]) + '\n' 870 else: 871 content_string += str(i + range[0]) 872 for element in sub_corr: 873 content_string += '\t' + ' ' * int(element >= 0) + str(element) 874 content_string += '\n' 875 return content_string 876 877 def __str__(self): 878 return self.__repr__() 879 880 # We define the basic operations, that can be performed with correlators. 881 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. 882 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. 883 # One could try and tell Obs to check if the y in __mul__ is a Corr and 884 885 def __add__(self, y): 886 if isinstance(y, Corr): 887 if ((self.N != y.N) or (self.T != y.T)): 888 raise Exception("Addition of Corrs with different shape") 889 newcontent = [] 890 for t in range(self.T): 891 if (self.content[t] is None) or (y.content[t] is None): 892 newcontent.append(None) 893 else: 894 newcontent.append(self.content[t] + y.content[t]) 895 return Corr(newcontent) 896 897 elif isinstance(y, (Obs, int, float, CObs)): 898 newcontent = [] 899 for t in range(self.T): 900 if (self.content[t] is None): 901 newcontent.append(None) 902 else: 903 newcontent.append(self.content[t] + y) 904 return Corr(newcontent, prange=self.prange) 905 elif isinstance(y, np.ndarray): 906 if y.shape == (self.T,): 907 return Corr(list((np.array(self.content).T + y).T)) 908 else: 909 raise ValueError("operands could not be broadcast together") 910 else: 911 raise TypeError("Corr + wrong type") 912 913 def __mul__(self, y): 914 if isinstance(y, Corr): 915 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 916 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 917 newcontent = [] 918 for t in range(self.T): 919 if (self.content[t] is None) or (y.content[t] is None): 920 newcontent.append(None) 921 else: 922 newcontent.append(self.content[t] * y.content[t]) 923 return Corr(newcontent) 924 925 elif isinstance(y, (Obs, int, float, CObs)): 926 newcontent = [] 927 for t in range(self.T): 928 if (self.content[t] is None): 929 newcontent.append(None) 930 else: 931 newcontent.append(self.content[t] * y) 932 return Corr(newcontent, prange=self.prange) 933 elif isinstance(y, np.ndarray): 934 if y.shape == (self.T,): 935 return Corr(list((np.array(self.content).T * y).T)) 936 else: 937 raise ValueError("operands could not be broadcast together") 938 else: 939 raise TypeError("Corr * wrong type") 940 941 def __truediv__(self, y): 942 if isinstance(y, Corr): 943 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 944 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 945 newcontent = [] 946 for t in range(self.T): 947 if (self.content[t] is None) or (y.content[t] is None): 948 newcontent.append(None) 949 else: 950 newcontent.append(self.content[t] / y.content[t]) 951 for t in range(self.T): 952 if newcontent[t] is None: 953 continue 954 if np.isnan(np.sum(newcontent[t]).value): 955 newcontent[t] = None 956 957 if all([item is None for item in newcontent]): 958 raise Exception("Division returns completely undefined correlator") 959 return Corr(newcontent) 960 961 elif isinstance(y, (Obs, CObs)): 962 if isinstance(y, Obs): 963 if y.value == 0: 964 raise Exception('Division by zero will return undefined correlator') 965 if isinstance(y, CObs): 966 if y.is_zero(): 967 raise Exception('Division by zero will return undefined correlator') 968 969 newcontent = [] 970 for t in range(self.T): 971 if (self.content[t] is None): 972 newcontent.append(None) 973 else: 974 newcontent.append(self.content[t] / y) 975 return Corr(newcontent, prange=self.prange) 976 977 elif isinstance(y, (int, float)): 978 if y == 0: 979 raise Exception('Division by zero will return undefined correlator') 980 newcontent = [] 981 for t in range(self.T): 982 if (self.content[t] is None): 983 newcontent.append(None) 984 else: 985 newcontent.append(self.content[t] / y) 986 return Corr(newcontent, prange=self.prange) 987 elif isinstance(y, np.ndarray): 988 if y.shape == (self.T,): 989 return Corr(list((np.array(self.content).T / y).T)) 990 else: 991 raise ValueError("operands could not be broadcast together") 992 else: 993 raise TypeError('Corr / wrong type') 994 995 def __neg__(self): 996 newcontent = [None if (item is None) else -1. * item for item in self.content] 997 return Corr(newcontent, prange=self.prange) 998 999 def __sub__(self, y): 1000 return self + (-y) 1001 1002 def __pow__(self, y): 1003 if isinstance(y, (Obs, int, float, CObs)): 1004 newcontent = [None if (item is None) else item**y for item in self.content] 1005 return Corr(newcontent, prange=self.prange) 1006 else: 1007 raise TypeError('Type of exponent not supported') 1008 1009 def __abs__(self): 1010 newcontent = [None if (item is None) else np.abs(item) for item in self.content] 1011 return Corr(newcontent, prange=self.prange) 1012 1013 # The numpy functions: 1014 def sqrt(self): 1015 return self**0.5 1016 1017 def log(self): 1018 newcontent = [None if (item is None) else np.log(item) for item in self.content] 1019 return Corr(newcontent, prange=self.prange) 1020 1021 def exp(self): 1022 newcontent = [None if (item is None) else np.exp(item) for item in self.content] 1023 return Corr(newcontent, prange=self.prange) 1024 1025 def _apply_func_to_corr(self, func): 1026 newcontent = [None if (item is None) else func(item) for item in self.content] 1027 for t in range(self.T): 1028 if newcontent[t] is None: 1029 continue 1030 if np.isnan(np.sum(newcontent[t]).value): 1031 newcontent[t] = None 1032 if all([item is None for item in newcontent]): 1033 raise Exception('Operation returns undefined correlator') 1034 return Corr(newcontent) 1035 1036 def sin(self): 1037 return self._apply_func_to_corr(np.sin) 1038 1039 def cos(self): 1040 return self._apply_func_to_corr(np.cos) 1041 1042 def tan(self): 1043 return self._apply_func_to_corr(np.tan) 1044 1045 def sinh(self): 1046 return self._apply_func_to_corr(np.sinh) 1047 1048 def cosh(self): 1049 return self._apply_func_to_corr(np.cosh) 1050 1051 def tanh(self): 1052 return self._apply_func_to_corr(np.tanh) 1053 1054 def arcsin(self): 1055 return self._apply_func_to_corr(np.arcsin) 1056 1057 def arccos(self): 1058 return self._apply_func_to_corr(np.arccos) 1059 1060 def arctan(self): 1061 return self._apply_func_to_corr(np.arctan) 1062 1063 def arcsinh(self): 1064 return self._apply_func_to_corr(np.arcsinh) 1065 1066 def arccosh(self): 1067 return self._apply_func_to_corr(np.arccosh) 1068 1069 def arctanh(self): 1070 return self._apply_func_to_corr(np.arctanh) 1071 1072 # Right hand side operations (require tweak in main module to work) 1073 def __radd__(self, y): 1074 return self + y 1075 1076 def __rsub__(self, y): 1077 return -self + y 1078 1079 def __rmul__(self, y): 1080 return self * y 1081 1082 def __rtruediv__(self, y): 1083 return (self / y) ** (-1) 1084 1085 @property 1086 def real(self): 1087 def return_real(obs_OR_cobs): 1088 if isinstance(obs_OR_cobs, CObs): 1089 return obs_OR_cobs.real 1090 else: 1091 return obs_OR_cobs 1092 1093 return self._apply_func_to_corr(return_real) 1094 1095 @property 1096 def imag(self): 1097 def return_imag(obs_OR_cobs): 1098 if isinstance(obs_OR_cobs, CObs): 1099 return obs_OR_cobs.imag 1100 else: 1101 return obs_OR_cobs * 0 # So it stays the right type 1102 1103 return self._apply_func_to_corr(return_imag) 1104 1105 1106def _sort_vectors(vec_set, ts): 1107 """Helper function used to find a set of Eigenvectors consistent over all timeslices""" 1108 reference_sorting = np.array(vec_set[ts]) 1109 N = reference_sorting.shape[0] 1110 sorted_vec_set = [] 1111 for t in range(len(vec_set)): 1112 if vec_set[t] is None: 1113 sorted_vec_set.append(None) 1114 elif not t == ts: 1115 perms = [list(o) for o in permutations([i for i in range(N)], N)] 1116 best_score = 0 1117 for perm in perms: 1118 current_score = 1 1119 for k in range(N): 1120 new_sorting = reference_sorting.copy() 1121 new_sorting[perm[k], :] = vec_set[t][k] 1122 current_score *= abs(np.linalg.det(new_sorting)) 1123 if current_score > best_score: 1124 best_score = current_score 1125 best_perm = perm 1126 sorted_vec_set.append([vec_set[t][k] for k in best_perm]) 1127 else: 1128 sorted_vec_set.append(vec_set[t]) 1129 1130 return sorted_vec_set 1131 1132 1133def _GEVP_solver(Gt, G0): # Just so normalization an sorting does not need to be repeated. Here we could later put in some checks 1134 sp_val, sp_vecs = scipy.linalg.eigh(Gt, G0) 1135 sp_vecs = [sp_vecs[:, np.argsort(sp_val)[-i]] for i in range(1, sp_vecs.shape[0] + 1)] 1136 sp_vecs = [v / np.sqrt((v.T @ G0 @ v)) for v in sp_vecs] 1137 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 else: 278 279 all_vecs = [] 280 for t in range(self.T): 281 try: 282 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") 283 for i in range(self.N): 284 for j in range(self.N): 285 G0[i, j] = self.content[t0][i, j].value 286 Gt[i, j] = self.content[t][i, j].value 287 288 sp_vecs = _GEVP_solver(Gt, G0) 289 if sorted_list == "Eigenvalue": 290 sp_vec = sp_vecs[state] 291 all_vecs.append(sp_vec) 292 else: 293 all_vecs.append(sp_vecs) 294 except Exception: 295 all_vecs.append(None) 296 if sorted_list == "Eigenvector": 297 if (ts is None): 298 raise Exception("ts is required for the Eigenvector sorting method.") 299 all_vecs = _sort_vectors(all_vecs, ts) 300 all_vecs = [a[state] for a in all_vecs] 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): 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 """ 729 if self.N != 1: 730 raise Exception("Correlator must be projected before plotting") 731 732 if auto_gamma: 733 self.gamma_method() 734 735 if x_range is None: 736 x_range = [0, self.T - 1] 737 738 fig = plt.figure() 739 ax1 = fig.add_subplot(111) 740 741 x, y, y_err = self.plottable() 742 ax1.errorbar(x, y, y_err, label=self.tag) 743 if logscale: 744 ax1.set_yscale('log') 745 else: 746 if y_range is None: 747 try: 748 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)]) 749 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)]) 750 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 751 except Exception: 752 pass 753 else: 754 ax1.set_ylim(y_range) 755 if comp: 756 if isinstance(comp, (Corr, list)): 757 for corr in comp if isinstance(comp, list) else [comp]: 758 if auto_gamma: 759 corr.gamma_method() 760 x, y, y_err = corr.plottable() 761 plt.errorbar(x, y, y_err, label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 762 else: 763 raise Exception("'comp' must be a correlator or a list of correlators.") 764 765 if plateau: 766 if isinstance(plateau, Obs): 767 if auto_gamma: 768 plateau.gamma_method() 769 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 770 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 771 else: 772 raise Exception("'plateau' must be an Obs") 773 if self.prange: 774 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 775 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 776 777 if fit_res: 778 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 779 ax1.plot(x_samples, 780 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 781 ls='-', marker=',', lw=2) 782 783 ax1.set_xlabel(r'$x_0 / a$') 784 if ylabel: 785 ax1.set_ylabel(ylabel) 786 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 787 788 handles, labels = ax1.get_legend_handles_labels() 789 if labels: 790 ax1.legend() 791 plt.draw() 792 793 if save: 794 if isinstance(save, str): 795 fig.savefig(save) 796 else: 797 raise Exception("'save' has to be a string.") 798 799 def spaghetti_plot(self, logscale=True): 800 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 801 802 Parameters 803 ---------- 804 logscale : bool 805 Determines whether the scale of the y-axis is logarithmic or standard. 806 """ 807 if self.N != 1: 808 raise Exception("Correlator needs to be projected first.") 809 810 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])) 811 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 812 813 for name in mc_names: 814 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 815 816 fig = plt.figure() 817 ax = fig.add_subplot(111) 818 for dat in data: 819 ax.plot(x0_vals, dat, ls='-', marker='') 820 821 if logscale is True: 822 ax.set_yscale('log') 823 824 ax.set_xlabel(r'$x_0 / a$') 825 plt.title(name) 826 plt.draw() 827 828 def dump(self, filename, datatype="json.gz", **kwargs): 829 """Dumps the Corr into a file of chosen type 830 Parameters 831 ---------- 832 filename : str 833 Name of the file to be saved. 834 datatype : str 835 Format of the exported file. Supported formats include 836 "json.gz" and "pickle" 837 path : str 838 specifies a custom path for the file (default '.') 839 """ 840 if datatype == "json.gz": 841 from .input.json import dump_to_json 842 if 'path' in kwargs: 843 file_name = kwargs.get('path') + '/' + filename 844 else: 845 file_name = filename 846 dump_to_json(self, file_name) 847 elif datatype == "pickle": 848 dump_object(self, filename, **kwargs) 849 else: 850 raise Exception("Unknown datatype " + str(datatype)) 851 852 def print(self, range=[0, None]): 853 print(self.__repr__(range)) 854 855 def __repr__(self, range=[0, None]): 856 content_string = "" 857 858 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 859 860 if self.tag is not None: 861 content_string += "Description: " + self.tag + "\n" 862 if self.N != 1: 863 return content_string 864 865 if range[1]: 866 range[1] += 1 867 content_string += 'x0/a\tCorr(x0/a)\n------------------\n' 868 for i, sub_corr in enumerate(self.content[range[0]:range[1]]): 869 if sub_corr is None: 870 content_string += str(i + range[0]) + '\n' 871 else: 872 content_string += str(i + range[0]) 873 for element in sub_corr: 874 content_string += '\t' + ' ' * int(element >= 0) + str(element) 875 content_string += '\n' 876 return content_string 877 878 def __str__(self): 879 return self.__repr__() 880 881 # We define the basic operations, that can be performed with correlators. 882 # While */+- get defined here, they only work for Corr*Obs and not Obs*Corr. 883 # This is because Obs*Corr checks Obs.__mul__ first and does not catch an exception. 884 # One could try and tell Obs to check if the y in __mul__ is a Corr and 885 886 def __add__(self, y): 887 if isinstance(y, Corr): 888 if ((self.N != y.N) or (self.T != y.T)): 889 raise Exception("Addition of Corrs with different shape") 890 newcontent = [] 891 for t in range(self.T): 892 if (self.content[t] is None) or (y.content[t] is None): 893 newcontent.append(None) 894 else: 895 newcontent.append(self.content[t] + y.content[t]) 896 return Corr(newcontent) 897 898 elif isinstance(y, (Obs, int, float, CObs)): 899 newcontent = [] 900 for t in range(self.T): 901 if (self.content[t] is None): 902 newcontent.append(None) 903 else: 904 newcontent.append(self.content[t] + y) 905 return Corr(newcontent, prange=self.prange) 906 elif isinstance(y, np.ndarray): 907 if y.shape == (self.T,): 908 return Corr(list((np.array(self.content).T + y).T)) 909 else: 910 raise ValueError("operands could not be broadcast together") 911 else: 912 raise TypeError("Corr + wrong type") 913 914 def __mul__(self, y): 915 if isinstance(y, Corr): 916 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 917 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 918 newcontent = [] 919 for t in range(self.T): 920 if (self.content[t] is None) or (y.content[t] is None): 921 newcontent.append(None) 922 else: 923 newcontent.append(self.content[t] * y.content[t]) 924 return Corr(newcontent) 925 926 elif isinstance(y, (Obs, int, float, CObs)): 927 newcontent = [] 928 for t in range(self.T): 929 if (self.content[t] is None): 930 newcontent.append(None) 931 else: 932 newcontent.append(self.content[t] * y) 933 return Corr(newcontent, prange=self.prange) 934 elif isinstance(y, np.ndarray): 935 if y.shape == (self.T,): 936 return Corr(list((np.array(self.content).T * y).T)) 937 else: 938 raise ValueError("operands could not be broadcast together") 939 else: 940 raise TypeError("Corr * wrong type") 941 942 def __truediv__(self, y): 943 if isinstance(y, Corr): 944 if not((self.N == 1 or y.N == 1 or self.N == y.N) and self.T == y.T): 945 raise Exception("Multiplication of Corr object requires N=N or N=1 and T=T") 946 newcontent = [] 947 for t in range(self.T): 948 if (self.content[t] is None) or (y.content[t] is None): 949 newcontent.append(None) 950 else: 951 newcontent.append(self.content[t] / y.content[t]) 952 for t in range(self.T): 953 if newcontent[t] is None: 954 continue 955 if np.isnan(np.sum(newcontent[t]).value): 956 newcontent[t] = None 957 958 if all([item is None for item in newcontent]): 959 raise Exception("Division returns completely undefined correlator") 960 return Corr(newcontent) 961 962 elif isinstance(y, (Obs, CObs)): 963 if isinstance(y, Obs): 964 if y.value == 0: 965 raise Exception('Division by zero will return undefined correlator') 966 if isinstance(y, CObs): 967 if y.is_zero(): 968 raise Exception('Division by zero will return undefined correlator') 969 970 newcontent = [] 971 for t in range(self.T): 972 if (self.content[t] is None): 973 newcontent.append(None) 974 else: 975 newcontent.append(self.content[t] / y) 976 return Corr(newcontent, prange=self.prange) 977 978 elif isinstance(y, (int, float)): 979 if y == 0: 980 raise Exception('Division by zero will return undefined correlator') 981 newcontent = [] 982 for t in range(self.T): 983 if (self.content[t] is None): 984 newcontent.append(None) 985 else: 986 newcontent.append(self.content[t] / y) 987 return Corr(newcontent, prange=self.prange) 988 elif isinstance(y, np.ndarray): 989 if y.shape == (self.T,): 990 return Corr(list((np.array(self.content).T / y).T)) 991 else: 992 raise ValueError("operands could not be broadcast together") 993 else: 994 raise TypeError('Corr / wrong type') 995 996 def __neg__(self): 997 newcontent = [None if (item is None) else -1. * item for item in self.content] 998 return Corr(newcontent, prange=self.prange) 999 1000 def __sub__(self, y): 1001 return self + (-y) 1002 1003 def __pow__(self, y): 1004 if isinstance(y, (Obs, int, float, CObs)): 1005 newcontent = [None if (item is None) else item**y for item in self.content] 1006 return Corr(newcontent, prange=self.prange) 1007 else: 1008 raise TypeError('Type of exponent not supported') 1009 1010 def __abs__(self): 1011 newcontent = [None if (item is None) else np.abs(item) for item in self.content] 1012 return Corr(newcontent, prange=self.prange) 1013 1014 # The numpy functions: 1015 def sqrt(self): 1016 return self**0.5 1017 1018 def log(self): 1019 newcontent = [None if (item is None) else np.log(item) for item in self.content] 1020 return Corr(newcontent, prange=self.prange) 1021 1022 def exp(self): 1023 newcontent = [None if (item is None) else np.exp(item) for item in self.content] 1024 return Corr(newcontent, prange=self.prange) 1025 1026 def _apply_func_to_corr(self, func): 1027 newcontent = [None if (item is None) else func(item) for item in self.content] 1028 for t in range(self.T): 1029 if newcontent[t] is None: 1030 continue 1031 if np.isnan(np.sum(newcontent[t]).value): 1032 newcontent[t] = None 1033 if all([item is None for item in newcontent]): 1034 raise Exception('Operation returns undefined correlator') 1035 return Corr(newcontent) 1036 1037 def sin(self): 1038 return self._apply_func_to_corr(np.sin) 1039 1040 def cos(self): 1041 return self._apply_func_to_corr(np.cos) 1042 1043 def tan(self): 1044 return self._apply_func_to_corr(np.tan) 1045 1046 def sinh(self): 1047 return self._apply_func_to_corr(np.sinh) 1048 1049 def cosh(self): 1050 return self._apply_func_to_corr(np.cosh) 1051 1052 def tanh(self): 1053 return self._apply_func_to_corr(np.tanh) 1054 1055 def arcsin(self): 1056 return self._apply_func_to_corr(np.arcsin) 1057 1058 def arccos(self): 1059 return self._apply_func_to_corr(np.arccos) 1060 1061 def arctan(self): 1062 return self._apply_func_to_corr(np.arctan) 1063 1064 def arcsinh(self): 1065 return self._apply_func_to_corr(np.arcsinh) 1066 1067 def arccosh(self): 1068 return self._apply_func_to_corr(np.arccosh) 1069 1070 def arctanh(self): 1071 return self._apply_func_to_corr(np.arctanh) 1072 1073 # Right hand side operations (require tweak in main module to work) 1074 def __radd__(self, y): 1075 return self + y 1076 1077 def __rsub__(self, y): 1078 return -self + y 1079 1080 def __rmul__(self, y): 1081 return self * y 1082 1083 def __rtruediv__(self, y): 1084 return (self / y) ** (-1) 1085 1086 @property 1087 def real(self): 1088 def return_real(obs_OR_cobs): 1089 if isinstance(obs_OR_cobs, CObs): 1090 return obs_OR_cobs.real 1091 else: 1092 return obs_OR_cobs 1093 1094 return self._apply_func_to_corr(return_real) 1095 1096 @property 1097 def imag(self): 1098 def return_imag(obs_OR_cobs): 1099 if isinstance(obs_OR_cobs, CObs): 1100 return obs_OR_cobs.imag 1101 else: 1102 return obs_OR_cobs * 0 # So it stays the right type 1103 1104 return self._apply_func_to_corr(return_imag)
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 else: 278 279 all_vecs = [] 280 for t in range(self.T): 281 try: 282 G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double") 283 for i in range(self.N): 284 for j in range(self.N): 285 G0[i, j] = self.content[t0][i, j].value 286 Gt[i, j] = self.content[t][i, j].value 287 288 sp_vecs = _GEVP_solver(Gt, G0) 289 if sorted_list == "Eigenvalue": 290 sp_vec = sp_vecs[state] 291 all_vecs.append(sp_vec) 292 else: 293 all_vecs.append(sp_vecs) 294 except Exception: 295 all_vecs.append(None) 296 if sorted_list == "Eigenvector": 297 if (ts is None): 298 raise Exception("ts is required for the Eigenvector sorting method.") 299 all_vecs = _sort_vectors(all_vecs, ts) 300 all_vecs = [a[state] for a in all_vecs] 301 302 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
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)
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
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)
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
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)
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
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)
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
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)
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
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
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
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.")
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
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.")
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
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.')
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
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
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
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)
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
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
Sets the attribute prange of the Corr object.
View Source
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): 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 """ 729 if self.N != 1: 730 raise Exception("Correlator must be projected before plotting") 731 732 if auto_gamma: 733 self.gamma_method() 734 735 if x_range is None: 736 x_range = [0, self.T - 1] 737 738 fig = plt.figure() 739 ax1 = fig.add_subplot(111) 740 741 x, y, y_err = self.plottable() 742 ax1.errorbar(x, y, y_err, label=self.tag) 743 if logscale: 744 ax1.set_yscale('log') 745 else: 746 if y_range is None: 747 try: 748 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)]) 749 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)]) 750 ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)]) 751 except Exception: 752 pass 753 else: 754 ax1.set_ylim(y_range) 755 if comp: 756 if isinstance(comp, (Corr, list)): 757 for corr in comp if isinstance(comp, list) else [comp]: 758 if auto_gamma: 759 corr.gamma_method() 760 x, y, y_err = corr.plottable() 761 plt.errorbar(x, y, y_err, label=corr.tag, mfc=plt.rcParams['axes.facecolor']) 762 else: 763 raise Exception("'comp' must be a correlator or a list of correlators.") 764 765 if plateau: 766 if isinstance(plateau, Obs): 767 if auto_gamma: 768 plateau.gamma_method() 769 ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau)) 770 ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-') 771 else: 772 raise Exception("'plateau' must be an Obs") 773 if self.prange: 774 ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',') 775 ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',') 776 777 if fit_res: 778 x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05) 779 ax1.plot(x_samples, 780 fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples), 781 ls='-', marker=',', lw=2) 782 783 ax1.set_xlabel(r'$x_0 / a$') 784 if ylabel: 785 ax1.set_ylabel(ylabel) 786 ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5]) 787 788 handles, labels = ax1.get_legend_handles_labels() 789 if labels: 790 ax1.legend() 791 plt.draw() 792 793 if save: 794 if isinstance(save, str): 795 fig.savefig(save) 796 else: 797 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.
View Source
799 def spaghetti_plot(self, logscale=True): 800 """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations. 801 802 Parameters 803 ---------- 804 logscale : bool 805 Determines whether the scale of the y-axis is logarithmic or standard. 806 """ 807 if self.N != 1: 808 raise Exception("Correlator needs to be projected first.") 809 810 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])) 811 x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None] 812 813 for name in mc_names: 814 data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T 815 816 fig = plt.figure() 817 ax = fig.add_subplot(111) 818 for dat in data: 819 ax.plot(x0_vals, dat, ls='-', marker='') 820 821 if logscale is True: 822 ax.set_yscale('log') 823 824 ax.set_xlabel(r'$x_0 / a$') 825 plt.title(name) 826 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
828 def dump(self, filename, datatype="json.gz", **kwargs): 829 """Dumps the Corr into a file of chosen type 830 Parameters 831 ---------- 832 filename : str 833 Name of the file to be saved. 834 datatype : str 835 Format of the exported file. Supported formats include 836 "json.gz" and "pickle" 837 path : str 838 specifies a custom path for the file (default '.') 839 """ 840 if datatype == "json.gz": 841 from .input.json import dump_to_json 842 if 'path' in kwargs: 843 file_name = kwargs.get('path') + '/' + filename 844 else: 845 file_name = filename 846 dump_to_json(self, file_name) 847 elif datatype == "pickle": 848 dump_object(self, filename, **kwargs) 849 else: 850 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 '.')