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

199    def symmetric(self):
 200        """ Symmetrize the correlator around x0=0."""
-201        if self.T % 2 != 0:
-202            raise Exception("Can not symmetrize odd T")
-203
-204        if np.argmax(np.abs(self.content)) != 0:
-205            warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
-206
-207        newcontent = [self.content[0]]
-208        for t in range(1, self.T):
-209            if (self.content[t] is None) or (self.content[self.T - t] is None):
-210                newcontent.append(None)
-211            else:
-212                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
-213        if(all([x is None for x in newcontent])):
-214            raise Exception("Corr could not be symmetrized: No redundant values")
-215        return Corr(newcontent, prange=self.prange)
+201        if self.N != 1:
+202            raise Exception('symmetric cannot be safely applied to multi-dimensional correlators.')
+203        if self.T % 2 != 0:
+204            raise Exception("Can not symmetrize odd T")
+205
+206        if np.argmax(np.abs(self.content)) != 0:
+207            warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
+208
+209        newcontent = [self.content[0]]
+210        for t in range(1, self.T):
+211            if (self.content[t] is None) or (self.content[self.T - t] is None):
+212                newcontent.append(None)
+213            else:
+214                newcontent.append(0.5 * (self.content[t] + self.content[self.T - t]))
+215        if(all([x is None for x in newcontent])):
+216            raise Exception("Corr could not be symmetrized: No redundant values")
+217        return Corr(newcontent, prange=self.prange)
@@ -2986,25 +3021,27 @@ timeslice and the error on each timeslice.

- 
217    def anti_symmetric(self):
-218        """Anti-symmetrize the correlator around x0=0."""
-219        if self.T % 2 != 0:
-220            raise Exception("Can not symmetrize odd T")
-221
-222        test = 1 * self
-223        test.gamma_method()
-224        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
-225            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning)
-226
-227        newcontent = [self.content[0]]
-228        for t in range(1, self.T):
-229            if (self.content[t] is None) or (self.content[self.T - t] is None):
-230                newcontent.append(None)
-231            else:
-232                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
-233        if(all([x is None for x in newcontent])):
-234            raise Exception("Corr could not be symmetrized: No redundant values")
-235        return Corr(newcontent, prange=self.prange)
+            
219    def anti_symmetric(self):
+220        """Anti-symmetrize the correlator around x0=0."""
+221        if self.N != 1:
+222            raise Exception('anti_symmetric cannot be safely applied to multi-dimensional correlators.')
+223        if self.T % 2 != 0:
+224            raise Exception("Can not symmetrize odd T")
+225
+226        test = 1 * self
+227        test.gamma_method()
+228        if not all([o.is_zero_within_error(3) for o in test.content[0]]):
+229            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning)
+230
+231        newcontent = [self.content[0]]
+232        for t in range(1, self.T):
+233            if (self.content[t] is None) or (self.content[self.T - t] is None):
+234                newcontent.append(None)
+235            else:
+236                newcontent.append(0.5 * (self.content[t] - self.content[self.T - t]))
+237        if(all([x is None for x in newcontent])):
+238            raise Exception("Corr could not be symmetrized: No redundant values")
+239        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -3024,13 +3061,13 @@ timeslice and the error on each timeslice.
- 
237    def matrix_symmetric(self):
-238        """Symmetrizes the correlator matrices on every timeslice."""
-239        if self.N > 1:
-240            transposed = [None if len(list(filter(None, np.asarray(G).flatten()))) < self.N ** 2 else G.T for G in self.content]
-241            return 0.5 * (Corr(transposed) + self)
-242        if self.N == 1:
-243            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
+            
241    def matrix_symmetric(self):
+242        """Symmetrizes the correlator matrices on every timeslice."""
+243        if self.N > 1:
+244            transposed = [None if _check_for_none(self, G) else G.T for G in self.content]
+245            return 0.5 * (Corr(transposed) + self)
+246        if self.N == 1:
+247            raise Exception("Trying to symmetrize a correlator matrix, that already has N=1.")
 

 
 
@@ -3050,87 +3087,87 @@ timeslice and the error on each timeslice.
- 
245    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
-246        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
-247
-248        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
-249        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
-250        ```python
-251        C.GEVP(t0=2)[0]  # Ground state vector(s)
-252        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
-253        ```
-254
-255        Parameters
-256        ----------
-257        t0 : int
-258            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
-259        ts : int
-260            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
-261            If sort="Eigenvector" it gives a reference point for the sorting method.
-262        sort : string
-263            If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned.
-264            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
-265            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
-266              The reference state is identified by its eigenvalue at $t=t_s$.
-267
-268        Other Parameters
-269        ----------------
-270        state : int
-271           Returns only the vector(s) for a specified state. The lowest state is zero.
-272        '''
-273
-274        if self.N == 1:
-275            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
-276        if ts is not None:
-277            if (ts <= t0):
-278                raise Exception("ts has to be larger than t0.")
-279
-280        if "sorted_list" in kwargs:
-281            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
-282            sort = kwargs.get("sorted_list")
+            
249    def GEVP(self, t0, ts=None, sort="Eigenvalue", **kwargs):
+250        r'''Solve the generalized eigenvalue problem on the correlator matrix and returns the corresponding eigenvectors.
+251
+252        The eigenvectors are sorted according to the descending eigenvalues, the zeroth eigenvector(s) correspond to the
+253        largest eigenvalue(s). The eigenvector(s) for the individual states can be accessed via slicing
+254        ```python
+255        C.GEVP(t0=2)[0]  # Ground state vector(s)
+256        C.GEVP(t0=2)[:3]  # Vectors for the lowest three states
+257        ```
+258
+259        Parameters
+260        ----------
+261        t0 : int
+262            The time t0 for the right hand side of the GEVP according to $G(t)v_i=\lambda_i G(t_0)v_i$
+263        ts : int
+264            fixed time $G(t_s)v_i=\lambda_i G(t_0)v_i$ if sort=None.
+265            If sort="Eigenvector" it gives a reference point for the sorting method.
+266        sort : string
+267            If this argument is set, a list of self.T vectors per state is returned. If it is set to None, only one vector is returned.
+268            - "Eigenvalue": The eigenvector is chosen according to which eigenvalue it belongs individually on every timeslice.
+269            - "Eigenvector": Use the method described in arXiv:2004.10472 to find the set of v(t) belonging to the state.
+270              The reference state is identified by its eigenvalue at $t=t_s$.
+271
+272        Other Parameters
+273        ----------------
+274        state : int
+275           Returns only the vector(s) for a specified state. The lowest state is zero.
+276        '''
+277
+278        if self.N == 1:
+279            raise Exception("GEVP methods only works on correlator matrices and not single correlators.")
+280        if ts is not None:
+281            if (ts <= t0):
+282                raise Exception("ts has to be larger than t0.")
 283
-284        symmetric_corr = self.matrix_symmetric()
-285        if sort is None:
-286            if (ts is None):
-287                raise Exception("ts is required if sort=None.")
-288            if (self.content[t0] is None) or (self.content[ts] is None):
-289                raise Exception("Corr not defined at t0/ts.")
-290            G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
-291            for i in range(self.N):
-292                for j in range(self.N):
-293                    G0[i, j] = symmetric_corr[t0][i, j].value
-294                    Gt[i, j] = symmetric_corr[ts][i, j].value
-295
-296            reordered_vecs = _GEVP_solver(Gt, G0)
-297
-298        elif sort in ["Eigenvalue", "Eigenvector"]:
-299            if sort == "Eigenvalue" and ts is not None:
-300                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
-301            all_vecs = [None] * (t0 + 1)
-302            for t in range(t0 + 1, self.T):
-303                try:
-304                    G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
-305                    for i in range(self.N):
-306                        for j in range(self.N):
-307                            G0[i, j] = symmetric_corr[t0][i, j].value
-308                            Gt[i, j] = symmetric_corr[t][i, j].value
-309
-310                    all_vecs.append(_GEVP_solver(Gt, G0))
-311                except Exception:
-312                    all_vecs.append(None)
-313            if sort == "Eigenvector":
-314                if (ts is None):
-315                    raise Exception("ts is required for the Eigenvector sorting method.")
-316                all_vecs = _sort_vectors(all_vecs, ts)
-317
-318            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
-319        else:
-320            raise Exception("Unkown value for 'sort'.")
+284        if "sorted_list" in kwargs:
+285            warnings.warn("Argument 'sorted_list' is deprecated, use 'sort' instead.", DeprecationWarning)
+286            sort = kwargs.get("sorted_list")
+287
+288        symmetric_corr = self.matrix_symmetric()
+289        if sort is None:
+290            if (ts is None):
+291                raise Exception("ts is required if sort=None.")
+292            if (self.content[t0] is None) or (self.content[ts] is None):
+293                raise Exception("Corr not defined at t0/ts.")
+294            G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
+295            for i in range(self.N):
+296                for j in range(self.N):
+297                    G0[i, j] = symmetric_corr[t0][i, j].value
+298                    Gt[i, j] = symmetric_corr[ts][i, j].value
+299
+300            reordered_vecs = _GEVP_solver(Gt, G0)
+301
+302        elif sort in ["Eigenvalue", "Eigenvector"]:
+303            if sort == "Eigenvalue" and ts is not None:
+304                warnings.warn("ts has no effect when sorting by eigenvalue is chosen.", RuntimeWarning)
+305            all_vecs = [None] * (t0 + 1)
+306            for t in range(t0 + 1, self.T):
+307                try:
+308                    G0, Gt = np.empty([self.N, self.N], dtype="double"), np.empty([self.N, self.N], dtype="double")
+309                    for i in range(self.N):
+310                        for j in range(self.N):
+311                            G0[i, j] = symmetric_corr[t0][i, j].value
+312                            Gt[i, j] = symmetric_corr[t][i, j].value
+313
+314                    all_vecs.append(_GEVP_solver(Gt, G0))
+315                except Exception:
+316                    all_vecs.append(None)
+317            if sort == "Eigenvector":
+318                if (ts is None):
+319                    raise Exception("ts is required for the Eigenvector sorting method.")
+320                all_vecs = _sort_vectors(all_vecs, ts)
 321
-322        if "state" in kwargs:
-323            return reordered_vecs[kwargs.get("state")]
-324        else:
-325            return reordered_vecs
+322            reordered_vecs = [[v[s] if v is not None else None for v in all_vecs] for s in range(self.N)]
+323        else:
+324            raise Exception("Unkown value for 'sort'.")
+325
+326        if "state" in kwargs:
+327            return reordered_vecs[kwargs.get("state")]
+328        else:
+329            return reordered_vecs
 

 
 
@@ -3181,18 +3218,18 @@ Returns only the vector(s) for a specified state. The lowest state is zero.
- 
327    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
-328        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
-329
-330        Parameters
-331        ----------
-332        state : int
-333            The state one is interested in ordered by energy. The lowest state is zero.
-334
-335        All other parameters are identical to the ones of Corr.GEVP.
-336        """
-337        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
-338        return self.projected(vec)
+            
331    def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"):
+332        """Determines the eigenvalue of the GEVP by solving and projecting the correlator
+333
+334        Parameters
+335        ----------
+336        state : int
+337            The state one is interested in ordered by energy. The lowest state is zero.
+338
+339        All other parameters are identical to the ones of Corr.GEVP.
+340        """
+341        vec = self.GEVP(t0, ts=ts, sort=sort)[state]
+342        return self.projected(vec)
 

 
 
@@ -3220,46 +3257,46 @@ The state one is interested in ordered by energy. The lowest state is zero.
- 
340    def Hankel(self, N, periodic=False):
-341        """Constructs an NxN Hankel matrix
-342
-343        C(t) c(t+1) ... c(t+n-1)
-344        C(t+1) c(t+2) ... c(t+n)
-345        .................
-346        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
-347
-348        Parameters
-349        ----------
-350        N : int
-351            Dimension of the Hankel matrix
-352        periodic : bool, optional
-353            determines whether the matrix is extended periodically
-354        """
-355
-356        if self.N != 1:
-357            raise Exception("Multi-operator Prony not implemented!")
-358
-359        array = np.empty([N, N], dtype="object")
-360        new_content = []
-361        for t in range(self.T):
-362            new_content.append(array.copy())
-363
-364        def wrap(i):
-365            while i >= self.T:
-366                i -= self.T
-367            return i
-368
-369        for t in range(self.T):
-370            for i in range(N):
-371                for j in range(N):
-372                    if periodic:
-373                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
-374                    elif (t + i + j) >= self.T:
-375                        new_content[t] = None
-376                    else:
-377                        new_content[t][i, j] = self.content[t + i + j][0]
-378
-379        return Corr(new_content)
+            
344    def Hankel(self, N, periodic=False):
+345        """Constructs an NxN Hankel matrix
+346
+347        C(t) c(t+1) ... c(t+n-1)
+348        C(t+1) c(t+2) ... c(t+n)
+349        .................
+350        C(t+(n-1)) c(t+n) ... c(t+2(n-1))
+351
+352        Parameters
+353        ----------
+354        N : int
+355            Dimension of the Hankel matrix
+356        periodic : bool, optional
+357            determines whether the matrix is extended periodically
+358        """
+359
+360        if self.N != 1:
+361            raise Exception("Multi-operator Prony not implemented!")
+362
+363        array = np.empty([N, N], dtype="object")
+364        new_content = []
+365        for t in range(self.T):
+366            new_content.append(array.copy())
+367
+368        def wrap(i):
+369            while i >= self.T:
+370                i -= self.T
+371            return i
+372
+373        for t in range(self.T):
+374            for i in range(N):
+375                for j in range(N):
+376                    if periodic:
+377                        new_content[t][i, j] = self.content[wrap(t + i + j)][0]
+378                    elif (t + i + j) >= self.T:
+379                        new_content[t] = None
+380                    else:
+381                        new_content[t][i, j] = self.content[t + i + j][0]
+382
+383        return Corr(new_content)
 

 
 
@@ -3293,15 +3330,15 @@ determines whether the matrix is extended periodically
- 
381    def roll(self, dt):
-382        """Periodically shift the correlator by dt timeslices
-383
-384        Parameters
-385        ----------
-386        dt : int
-387            number of timeslices
-388        """
-389        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
+            
385    def roll(self, dt):
+386        """Periodically shift the correlator by dt timeslices
+387
+388        Parameters
+389        ----------
+390        dt : int
+391            number of timeslices
+392        """
+393        return Corr(list(np.roll(np.array(self.content, dtype=object), dt)))
 

 
 
@@ -3328,9 +3365,9 @@ number of timeslices
- 
391    def reverse(self):
-392        """Reverse the time ordering of the Corr"""
-393        return Corr(self.content[:: -1])
+            
395    def reverse(self):
+396        """Reverse the time ordering of the Corr"""
+397        return Corr(self.content[:: -1])
 

 
 
@@ -3350,23 +3387,23 @@ number of timeslices
- 
395    def thin(self, spacing=2, offset=0):
-396        """Thin out a correlator to suppress correlations
-397
-398        Parameters
-399        ----------
-400        spacing : int
-401            Keep only every 'spacing'th entry of the correlator
-402        offset : int
-403            Offset the equal spacing
-404        """
-405        new_content = []
-406        for t in range(self.T):
-407            if (offset + t) % spacing != 0:
-408                new_content.append(None)
-409            else:
-410                new_content.append(self.content[t])
-411        return Corr(new_content)
+            
399    def thin(self, spacing=2, offset=0):
+400        """Thin out a correlator to suppress correlations
+401
+402        Parameters
+403        ----------
+404        spacing : int
+405            Keep only every 'spacing'th entry of the correlator
+406        offset : int
+407            Offset the equal spacing
+408        """
+409        new_content = []
+410        for t in range(self.T):
+411            if (offset + t) % spacing != 0:
+412                new_content.append(None)
+413            else:
+414                new_content.append(self.content[t])
+415        return Corr(new_content)
 

 
 
@@ -3395,32 +3432,34 @@ Offset the equal spacing
- 
413    def correlate(self, partner):
-414        """Correlate the correlator with another correlator or Obs
-415
-416        Parameters
-417        ----------
-418        partner : Obs or Corr
-419            partner to correlate the correlator with.
-420            Can either be an Obs which is correlated with all entries of the
-421            correlator or a Corr of same length.
-422        """
-423        new_content = []
-424        for x0, t_slice in enumerate(self.content):
-425            if t_slice is None:
-426                new_content.append(None)
-427            else:
-428                if isinstance(partner, Corr):
-429                    if partner.content[x0] is None:
-430                        new_content.append(None)
-431                    else:
-432                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
-433                elif isinstance(partner, Obs):  # Should this include CObs?
-434                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
-435                else:
-436                    raise Exception("Can only correlate with an Obs or a Corr.")
-437
-438        return Corr(new_content)
+            
417    def correlate(self, partner):
+418        """Correlate the correlator with another correlator or Obs
+419
+420        Parameters
+421        ----------
+422        partner : Obs or Corr
+423            partner to correlate the correlator with.
+424            Can either be an Obs which is correlated with all entries of the
+425            correlator or a Corr of same length.
+426        """
+427        if self.N != 1:
+428            raise Exception("Only one-dimensional correlators can be safely correlated.")
+429        new_content = []
+430        for x0, t_slice in enumerate(self.content):
+431            if _check_for_none(self, t_slice):
+432                new_content.append(None)
+433            else:
+434                if isinstance(partner, Corr):
+435                    if _check_for_none(partner, partner.content[x0]):
+436                        new_content.append(None)
+437                    else:
+438                        new_content.append(np.array([correlate(o, partner.content[x0][0]) for o in t_slice]))
+439                elif isinstance(partner, Obs):  # Should this include CObs?
+440                    new_content.append(np.array([correlate(o, partner) for o in t_slice]))
+441                else:
+442                    raise Exception("Can only correlate with an Obs or a Corr.")
+443
+444        return Corr(new_content)
 

 
 
@@ -3449,26 +3488,28 @@ correlator or a Corr of same length.
- 
440    def reweight(self, weight, **kwargs):
-441        """Reweight the correlator.
-442
-443        Parameters
-444        ----------
-445        weight : Obs
-446            Reweighting factor. An Observable that has to be defined on a superset of the
-447            configurations in obs[i].idl for all i.
-448        all_configs : bool
-449            if True, the reweighted observables are normalized by the average of
-450            the reweighting factor on all configurations in weight.idl and not
-451            on the configurations in obs[i].idl.
-452        """
-453        new_content = []
-454        for t_slice in self.content:
-455            if t_slice is None:
-456                new_content.append(None)
-457            else:
-458                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
-459        return Corr(new_content)
+            
446    def reweight(self, weight, **kwargs):
+447        """Reweight the correlator.
+448
+449        Parameters
+450        ----------
+451        weight : Obs
+452            Reweighting factor. An Observable that has to be defined on a superset of the
+453            configurations in obs[i].idl for all i.
+454        all_configs : bool
+455            if True, the reweighted observables are normalized by the average of
+456            the reweighting factor on all configurations in weight.idl and not
+457            on the configurations in obs[i].idl.
+458        """
+459        if self.N != 1:
+460            raise Exception("Reweighting only implemented for one-dimensional correlators.")
+461        new_content = []
+462        for t_slice in self.content:
+463            if _check_for_none(self, t_slice):
+464                new_content.append(None)
+465            else:
+466                new_content.append(np.array(reweight(weight, t_slice, **kwargs)))
+467        return Corr(new_content)
 

 
 
@@ -3500,33 +3541,35 @@ on the configurations in obs[i].idl.
- 
461    def T_symmetry(self, partner, parity=+1):
-462        """Return the time symmetry average of the correlator and its partner
-463
-464        Parameters
-465        ----------
-466        partner : Corr
-467            Time symmetry partner of the Corr
-468        partity : int
-469            Parity quantum number of the correlator, can be +1 or -1
-470        """
-471        if not isinstance(partner, Corr):
-472            raise Exception("T partner has to be a Corr object.")
-473        if parity not in [+1, -1]:
-474            raise Exception("Parity has to be +1 or -1.")
-475        T_partner = parity * partner.reverse()
-476
-477        t_slices = []
-478        test = (self - T_partner)
-479        test.gamma_method()
-480        for x0, t_slice in enumerate(test.content):
-481            if t_slice is not None:
-482                if not t_slice[0].is_zero_within_error(5):
-483                    t_slices.append(x0)
-484        if t_slices:
-485            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
+            
469    def T_symmetry(self, partner, parity=+1):
+470        """Return the time symmetry average of the correlator and its partner
+471
+472        Parameters
+473        ----------
+474        partner : Corr
+475            Time symmetry partner of the Corr
+476        partity : int
+477            Parity quantum number of the correlator, can be +1 or -1
+478        """
+479        if self.N != 1:
+480            raise Exception("T_symmetry only implemented for one-dimensional correlators.")
+481        if not isinstance(partner, Corr):
+482            raise Exception("T partner has to be a Corr object.")
+483        if parity not in [+1, -1]:
+484            raise Exception("Parity has to be +1 or -1.")
+485        T_partner = parity * partner.reverse()
 486
-487        return (self + T_partner) / 2
+487        t_slices = []
+488        test = (self - T_partner)
+489        test.gamma_method()
+490        for x0, t_slice in enumerate(test.content):
+491            if t_slice is not None:
+492                if not t_slice[0].is_zero_within_error(5):
+493                    t_slices.append(x0)
+494        if t_slices:
+495            warnings.warn("T symmetry partners do not agree within 5 sigma on time slices " + str(t_slices) + ".", RuntimeWarning)
+496
+497        return (self + T_partner) / 2
 

 
 
@@ -3555,57 +3598,59 @@ Parity quantum number of the correlator, can be +1 or -1
- 
489    def deriv(self, variant="symmetric"):
-490        """Return the first derivative of the correlator with respect to x0.
-491
-492        Parameters
-493        ----------
-494        variant : str
-495            decides which definition of the finite differences derivative is used.
-496            Available choice: symmetric, forward, backward, improved, default: symmetric
-497        """
-498        if variant == "symmetric":
-499            newcontent = []
-500            for t in range(1, self.T - 1):
-501                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-502                    newcontent.append(None)
-503                else:
-504                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
-505            if(all([x is None for x in newcontent])):
-506                raise Exception('Derivative is undefined at all timeslices')
-507            return Corr(newcontent, padding=[1, 1])
-508        elif variant == "forward":
-509            newcontent = []
-510            for t in range(self.T - 1):
-511                if (self.content[t] is None) or (self.content[t + 1] is None):
-512                    newcontent.append(None)
-513                else:
-514                    newcontent.append(self.content[t + 1] - self.content[t])
-515            if(all([x is None for x in newcontent])):
-516                raise Exception("Derivative is undefined at all timeslices")
-517            return Corr(newcontent, padding=[0, 1])
-518        elif variant == "backward":
-519            newcontent = []
-520            for t in range(1, self.T):
-521                if (self.content[t - 1] is None) or (self.content[t] is None):
-522                    newcontent.append(None)
-523                else:
-524                    newcontent.append(self.content[t] - self.content[t - 1])
-525            if(all([x is None for x in newcontent])):
-526                raise Exception("Derivative is undefined at all timeslices")
-527            return Corr(newcontent, padding=[1, 0])
-528        elif variant == "improved":
-529            newcontent = []
-530            for t in range(2, self.T - 2):
-531                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):
-532                    newcontent.append(None)
-533                else:
-534                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
-535            if(all([x is None for x in newcontent])):
-536                raise Exception('Derivative is undefined at all timeslices')
-537            return Corr(newcontent, padding=[2, 2])
-538        else:
-539            raise Exception("Unknown variant.")
+            
499    def deriv(self, variant="symmetric"):
+500        """Return the first derivative of the correlator with respect to x0.
+501
+502        Parameters
+503        ----------
+504        variant : str
+505            decides which definition of the finite differences derivative is used.
+506            Available choice: symmetric, forward, backward, improved, default: symmetric
+507        """
+508        if self.N != 1:
+509            raise Exception("deriv only implemented for one-dimensional correlators.")
+510        if variant == "symmetric":
+511            newcontent = []
+512            for t in range(1, self.T - 1):
+513                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+514                    newcontent.append(None)
+515                else:
+516                    newcontent.append(0.5 * (self.content[t + 1] - self.content[t - 1]))
+517            if(all([x is None for x in newcontent])):
+518                raise Exception('Derivative is undefined at all timeslices')
+519            return Corr(newcontent, padding=[1, 1])
+520        elif variant == "forward":
+521            newcontent = []
+522            for t in range(self.T - 1):
+523                if (self.content[t] is None) or (self.content[t + 1] is None):
+524                    newcontent.append(None)
+525                else:
+526                    newcontent.append(self.content[t + 1] - self.content[t])
+527            if(all([x is None for x in newcontent])):
+528                raise Exception("Derivative is undefined at all timeslices")
+529            return Corr(newcontent, padding=[0, 1])
+530        elif variant == "backward":
+531            newcontent = []
+532            for t in range(1, self.T):
+533                if (self.content[t - 1] is None) or (self.content[t] is None):
+534                    newcontent.append(None)
+535                else:
+536                    newcontent.append(self.content[t] - self.content[t - 1])
+537            if(all([x is None for x in newcontent])):
+538                raise Exception("Derivative is undefined at all timeslices")
+539            return Corr(newcontent, padding=[1, 0])
+540        elif variant == "improved":
+541            newcontent = []
+542            for t in range(2, self.T - 2):
+543                if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
+544                    newcontent.append(None)
+545                else:
+546                    newcontent.append((1 / 12) * (self.content[t - 2] - 8 * self.content[t - 1] + 8 * self.content[t + 1] - self.content[t + 2]))
+547            if(all([x is None for x in newcontent])):
+548                raise Exception('Derivative is undefined at all timeslices')
+549            return Corr(newcontent, padding=[2, 2])
+550        else:
+551            raise Exception("Unknown variant.")
 

 
 
@@ -3633,37 +3678,39 @@ Available choice: symmetric, forward, backward, improved, default: symmetric
     
-            
541    def second_deriv(self, variant="symmetric"):
-542        """Return the second derivative of the correlator with respect to x0.
-543
-544        Parameters
-545        ----------
-546        variant : str
-547            decides which definition of the finite differences derivative is used.
-548            Available choice: symmetric, improved, default: symmetric
-549        """
-550        if variant == "symmetric":
-551            newcontent = []
-552            for t in range(1, self.T - 1):
-553                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
-554                    newcontent.append(None)
-555                else:
-556                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
-557            if(all([x is None for x in newcontent])):
-558                raise Exception("Derivative is undefined at all timeslices")
-559            return Corr(newcontent, padding=[1, 1])
-560        elif variant == "improved":
-561            newcontent = []
-562            for t in range(2, self.T - 2):
-563                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):
-564                    newcontent.append(None)
-565                else:
-566                    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]))
-567            if(all([x is None for x in newcontent])):
-568                raise Exception("Derivative is undefined at all timeslices")
-569            return Corr(newcontent, padding=[2, 2])
-570        else:
-571            raise Exception("Unknown variant.")
+            
553    def second_deriv(self, variant="symmetric"):
+554        """Return the second derivative of the correlator with respect to x0.
+555
+556        Parameters
+557        ----------
+558        variant : str
+559            decides which definition of the finite differences derivative is used.
+560            Available choice: symmetric, improved, default: symmetric
+561        """
+562        if self.N != 1:
+563            raise Exception("second_deriv only implemented for one-dimensional correlators.")
+564        if variant == "symmetric":
+565            newcontent = []
+566            for t in range(1, self.T - 1):
+567                if (self.content[t - 1] is None) or (self.content[t + 1] is None):
+568                    newcontent.append(None)
+569                else:
+570                    newcontent.append((self.content[t + 1] - 2 * self.content[t] + self.content[t - 1]))
+571            if(all([x is None for x in newcontent])):
+572                raise Exception("Derivative is undefined at all timeslices")
+573            return Corr(newcontent, padding=[1, 1])
+574        elif variant == "improved":
+575            newcontent = []
+576            for t in range(2, self.T - 2):
+577                if (self.content[t - 2] is None) or (self.content[t - 1] is None) or (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t + 2] is None):
+578                    newcontent.append(None)
+579                else:
+580                    newcontent.append((1 / 12) * (-self.content[t + 2] + 16 * self.content[t + 1] - 30 * self.content[t] + 16 * self.content[t - 1] - self.content[t - 2]))
+581            if(all([x is None for x in newcontent])):
+582                raise Exception("Derivative is undefined at all timeslices")
+583            return Corr(newcontent, padding=[2, 2])
+584        else:
+585            raise Exception("Unknown variant.")
 

 
 
@@ -3691,70 +3738,70 @@ Available choice: symmetric, improved, default: symmetric
 
     

     
-            
573    def m_eff(self, variant='log', guess=1.0):
-574        """Returns the effective mass of the correlator as correlator object
-575
-576        Parameters
-577        ----------
-578        variant : str
-579            log : uses the standard effective mass log(C(t) / C(t+1))
-580            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.
-581            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.
-582            See, e.g., arXiv:1205.5380
-583            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
-584        guess : float
-585            guess for the root finder, only relevant for the root variant
-586        """
-587        if self.N != 1:
-588            raise Exception('Correlator must be projected before getting m_eff')
-589        if variant == 'log':
-590            newcontent = []
-591            for t in range(self.T - 1):
-592                if (self.content[t] is None) or (self.content[t + 1] is None):
-593                    newcontent.append(None)
-594                else:
-595                    newcontent.append(self.content[t] / self.content[t + 1])
-596            if(all([x is None for x in newcontent])):
-597                raise Exception('m_eff is undefined at all timeslices')
-598
-599            return np.log(Corr(newcontent, padding=[0, 1]))
-600
-601        elif variant in ['periodic', 'cosh', 'sinh']:
-602            if variant in ['periodic', 'cosh']:
-603                func = anp.cosh
-604            else:
-605                func = anp.sinh
-606
-607            def root_function(x, d):
-608                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
-609
-610            newcontent = []
-611            for t in range(self.T - 1):
-612                if (self.content[t] is None) or (self.content[t + 1] is None):
-613                    newcontent.append(None)
-614                # Fill the two timeslices in the middle of the lattice with their predecessors
-615                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
-616                    newcontent.append(newcontent[-1])
-617                else:
-618                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
-619            if(all([x is None for x in newcontent])):
-620                raise Exception('m_eff is undefined at all timeslices')
-621
-622            return Corr(newcontent, padding=[0, 1])
+            
587    def m_eff(self, variant='log', guess=1.0):
+588        """Returns the effective mass of the correlator as correlator object
+589
+590        Parameters
+591        ----------
+592        variant : str
+593            log : uses the standard effective mass log(C(t) / C(t+1))
+594            cosh, periodic : Use periodicitiy of the correlator by solving C(t) / C(t+1) = cosh(m * (t - T/2)) / cosh(m * (t + 1 - T/2)) for m.
+595            sinh : Use anti-periodicitiy of the correlator by solving C(t) / C(t+1) = sinh(m * (t - T/2)) / sinh(m * (t + 1 - T/2)) for m.
+596            See, e.g., arXiv:1205.5380
+597            arccosh : Uses the explicit form of the symmetrized correlator (not recommended)
+598        guess : float
+599            guess for the root finder, only relevant for the root variant
+600        """
+601        if self.N != 1:
+602            raise Exception('Correlator must be projected before getting m_eff')
+603        if variant == 'log':
+604            newcontent = []
+605            for t in range(self.T - 1):
+606                if (self.content[t] is None) or (self.content[t + 1] is None):
+607                    newcontent.append(None)
+608                else:
+609                    newcontent.append(self.content[t] / self.content[t + 1])
+610            if(all([x is None for x in newcontent])):
+611                raise Exception('m_eff is undefined at all timeslices')
+612
+613            return np.log(Corr(newcontent, padding=[0, 1]))
+614
+615        elif variant in ['periodic', 'cosh', 'sinh']:
+616            if variant in ['periodic', 'cosh']:
+617                func = anp.cosh
+618            else:
+619                func = anp.sinh
+620
+621            def root_function(x, d):
+622                return func(x * (t - self.T / 2)) / func(x * (t + 1 - self.T / 2)) - d
 623
-624        elif variant == 'arccosh':
-625            newcontent = []
-626            for t in range(1, self.T - 1):
-627                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
-628                    newcontent.append(None)
-629                else:
-630                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
-631            if(all([x is None for x in newcontent])):
-632                raise Exception("m_eff is undefined at all timeslices")
-633            return np.arccosh(Corr(newcontent, padding=[1, 1]))
-634
-635        else:
-636            raise Exception('Unknown variant.')
+624            newcontent = []
+625            for t in range(self.T - 1):
+626                if (self.content[t] is None) or (self.content[t + 1] is None):
+627                    newcontent.append(None)
+628                # Fill the two timeslices in the middle of the lattice with their predecessors
+629                elif variant == 'sinh' and t in [self.T / 2, self.T / 2 - 1]:
+630                    newcontent.append(newcontent[-1])
+631                else:
+632                    newcontent.append(np.abs(find_root(self.content[t][0] / self.content[t + 1][0], root_function, guess=guess)))
+633            if(all([x is None for x in newcontent])):
+634                raise Exception('m_eff is undefined at all timeslices')
+635
+636            return Corr(newcontent, padding=[0, 1])
+637
+638        elif variant == 'arccosh':
+639            newcontent = []
+640            for t in range(1, self.T - 1):
+641                if (self.content[t] is None) or (self.content[t + 1] is None) or (self.content[t - 1] is None):
+642                    newcontent.append(None)
+643                else:
+644                    newcontent.append((self.content[t + 1] + self.content[t - 1]) / (2 * self.content[t]))
+645            if(all([x is None for x in newcontent])):
+646                raise Exception("m_eff is undefined at all timeslices")
+647            return np.arccosh(Corr(newcontent, padding=[1, 1]))
+648
+649        else:
+650            raise Exception('Unknown variant.')
 

 
 
@@ -3787,39 +3834,39 @@ guess for the root finder, only relevant for the root variant
 
     

     
-            
638    def fit(self, function, fitrange=None, silent=False, **kwargs):
-639        r'''Fits function to the data
-640
-641        Parameters
-642        ----------
-643        function : obj
-644            function to fit to the data. See fits.least_squares for details.
-645        fitrange : list
-646            Two element list containing the timeslices on which the fit is supposed to start and stop.
-647            Caution: This range is inclusive as opposed to standard python indexing.
-648            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
-649            If not specified, self.prange or all timeslices are used.
-650        silent : bool
-651            Decides whether output is printed to the standard output.
-652        '''
-653        if self.N != 1:
-654            raise Exception("Correlator must be projected before fitting")
-655
-656        if fitrange is None:
-657            if self.prange:
-658                fitrange = self.prange
-659            else:
-660                fitrange = [0, self.T - 1]
-661        else:
-662            if not isinstance(fitrange, list):
-663                raise Exception("fitrange has to be a list with two elements")
-664            if len(fitrange) != 2:
-665                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
-666
-667        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-668        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
-669        result = least_squares(xs, ys, function, silent=silent, **kwargs)
-670        return result
+            
652    def fit(self, function, fitrange=None, silent=False, **kwargs):
+653        r'''Fits function to the data
+654
+655        Parameters
+656        ----------
+657        function : obj
+658            function to fit to the data. See fits.least_squares for details.
+659        fitrange : list
+660            Two element list containing the timeslices on which the fit is supposed to start and stop.
+661            Caution: This range is inclusive as opposed to standard python indexing.
+662            `fitrange=[4, 6]` corresponds to the three entries 4, 5 and 6.
+663            If not specified, self.prange or all timeslices are used.
+664        silent : bool
+665            Decides whether output is printed to the standard output.
+666        '''
+667        if self.N != 1:
+668            raise Exception("Correlator must be projected before fitting")
+669
+670        if fitrange is None:
+671            if self.prange:
+672                fitrange = self.prange
+673            else:
+674                fitrange = [0, self.T - 1]
+675        else:
+676            if not isinstance(fitrange, list):
+677                raise Exception("fitrange has to be a list with two elements")
+678            if len(fitrange) != 2:
+679                raise Exception("fitrange has to have exactly two elements [fit_start, fit_stop]")
+680
+681        xs = [x for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+682        ys = [self.content[x][0] for x in range(fitrange[0], fitrange[1] + 1) if not self.content[x] is None]
+683        result = least_squares(xs, ys, function, silent=silent, **kwargs)
+684        return result
 

 
 
@@ -3853,42 +3900,42 @@ Decides whether output is printed to the standard output.
 
     

     
-            
672    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
-673        """ Extract a plateau value from a Corr object
-674
-675        Parameters
-676        ----------
-677        plateau_range : list
-678            list with two entries, indicating the first and the last timeslice
-679            of the plateau region.
-680        method : str
-681            method to extract the plateau.
-682                'fit' fits a constant to the plateau region
-683                'avg', 'average' or 'mean' just average over the given timeslices.
-684        auto_gamma : bool
-685            apply gamma_method with default parameters to the Corr. Defaults to None
-686        """
-687        if not plateau_range:
-688            if self.prange:
-689                plateau_range = self.prange
-690            else:
-691                raise Exception("no plateau range provided")
-692        if self.N != 1:
-693            raise Exception("Correlator must be projected before getting a plateau.")
-694        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
-695            raise Exception("plateau is undefined at all timeslices in plateaurange.")
-696        if auto_gamma:
-697            self.gamma_method()
-698        if method == "fit":
-699            def const_func(a, t):
-700                return a[0]
-701            return self.fit(const_func, plateau_range)[0]
-702        elif method in ["avg", "average", "mean"]:
-703            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
-704            return returnvalue
-705
-706        else:
-707            raise Exception("Unsupported plateau method: " + method)
+            
686    def plateau(self, plateau_range=None, method="fit", auto_gamma=False):
+687        """ Extract a plateau value from a Corr object
+688
+689        Parameters
+690        ----------
+691        plateau_range : list
+692            list with two entries, indicating the first and the last timeslice
+693            of the plateau region.
+694        method : str
+695            method to extract the plateau.
+696                'fit' fits a constant to the plateau region
+697                'avg', 'average' or 'mean' just average over the given timeslices.
+698        auto_gamma : bool
+699            apply gamma_method with default parameters to the Corr. Defaults to None
+700        """
+701        if not plateau_range:
+702            if self.prange:
+703                plateau_range = self.prange
+704            else:
+705                raise Exception("no plateau range provided")
+706        if self.N != 1:
+707            raise Exception("Correlator must be projected before getting a plateau.")
+708        if(all([self.content[t] is None for t in range(plateau_range[0], plateau_range[1] + 1)])):
+709            raise Exception("plateau is undefined at all timeslices in plateaurange.")
+710        if auto_gamma:
+711            self.gamma_method()
+712        if method == "fit":
+713            def const_func(a, t):
+714                return a[0]
+715            return self.fit(const_func, plateau_range)[0]
+716        elif method in ["avg", "average", "mean"]:
+717            returnvalue = np.mean([item[0] for item in self.content[plateau_range[0]:plateau_range[1] + 1] if item is not None])
+718            return returnvalue
+719
+720        else:
+721            raise Exception("Unsupported plateau method: " + method)
 

 
 
@@ -3922,17 +3969,17 @@ apply gamma_method with default parameters to the Corr. Defaults to None
 
     

     
-            
709    def set_prange(self, prange):
-710        """Sets the attribute prange of the Corr object."""
-711        if not len(prange) == 2:
-712            raise Exception("prange must be a list or array with two values")
-713        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
-714            raise Exception("Start and end point must be integers")
-715        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
-716            raise Exception("Start and end point must define a range in the interval 0,T")
-717
-718        self.prange = prange
-719        return
+            
723    def set_prange(self, prange):
+724        """Sets the attribute prange of the Corr object."""
+725        if not len(prange) == 2:
+726            raise Exception("prange must be a list or array with two values")
+727        if not ((isinstance(prange[0], int)) and (isinstance(prange[1], int))):
+728            raise Exception("Start and end point must be integers")
+729        if not (0 <= prange[0] <= self.T and 0 <= prange[1] <= self.T and prange[0] < prange[1]):
+730            raise Exception("Start and end point must define a range in the interval 0,T")
+731
+732        self.prange = prange
+733        return
 

 
 
@@ -3952,118 +3999,118 @@ apply gamma_method with default parameters to the Corr. Defaults to None
 
     

     
-            
721    def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None):
-722        """Plots the correlator using the tag of the correlator as label if available.
-723
-724        Parameters
-725        ----------
-726        x_range : list
-727            list of two values, determining the range of the x-axis e.g. [4, 8]
-728        comp : Corr or list of Corr
-729            Correlator or list of correlators which are plotted for comparison.
-730            The tags of these correlators are used as labels if available.
-731        logscale : bool
-732            Sets y-axis to logscale
-733        plateau : Obs
-734            Plateau value to be visualized in the figure
-735        fit_res : Fit_result
-736            Fit_result object to be visualized
-737        ylabel : str
-738            Label for the y-axis
-739        save : str
-740            path to file in which the figure should be saved
-741        auto_gamma : bool
-742            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
-743        hide_sigma : float
-744            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
-745        references : list
-746            List of floating point values that are displayed as horizontal lines for reference.
-747        """
-748        if self.N != 1:
-749            raise Exception("Correlator must be projected before plotting")
-750
-751        if auto_gamma:
-752            self.gamma_method()
-753
-754        if x_range is None:
-755            x_range = [0, self.T - 1]
-756
-757        fig = plt.figure()
-758        ax1 = fig.add_subplot(111)
-759
-760        x, y, y_err = self.plottable()
-761        if hide_sigma:
-762            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
-763        else:
-764            hide_from = None
-765        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
-766        if logscale:
-767            ax1.set_yscale('log')
-768        else:
-769            if y_range is None:
-770                try:
-771                    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)])
-772                    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)])
-773                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
-774                except Exception:
-775                    pass
-776            else:
-777                ax1.set_ylim(y_range)
-778        if comp:
-779            if isinstance(comp, (Corr, list)):
-780                for corr in comp if isinstance(comp, list) else [comp]:
-781                    if auto_gamma:
-782                        corr.gamma_method()
-783                    x, y, y_err = corr.plottable()
-784                    if hide_sigma:
-785                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
-786                    else:
-787                        hide_from = None
-788                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
-789            else:
-790                raise Exception("'comp' must be a correlator or a list of correlators.")
-791
-792        if plateau:
-793            if isinstance(plateau, Obs):
-794                if auto_gamma:
-795                    plateau.gamma_method()
-796                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
-797                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
-798            else:
-799                raise Exception("'plateau' must be an Obs")
-800
-801        if references:
-802            if isinstance(references, list):
-803                for ref in references:
-804                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
-805            else:
-806                raise Exception("'references' must be a list of floating pint values.")
-807
-808        if self.prange:
-809            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
-810            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
-811
-812        if fit_res:
-813            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
-814            ax1.plot(x_samples,
-815                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
-816                     ls='-', marker=',', lw=2)
-817
-818        ax1.set_xlabel(r'$x_0 / a$')
-819        if ylabel:
-820            ax1.set_ylabel(ylabel)
-821        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
-822
-823        handles, labels = ax1.get_legend_handles_labels()
-824        if labels:
-825            ax1.legend()
-826        plt.draw()
-827
-828        if save:
-829            if isinstance(save, str):
-830                fig.savefig(save)
-831            else:
-832                raise Exception("'save' has to be a string.")
+            
735    def show(self, x_range=None, comp=None, y_range=None, logscale=False, plateau=None, fit_res=None, ylabel=None, save=None, auto_gamma=False, hide_sigma=None, references=None):
+736        """Plots the correlator using the tag of the correlator as label if available.
+737
+738        Parameters
+739        ----------
+740        x_range : list
+741            list of two values, determining the range of the x-axis e.g. [4, 8]
+742        comp : Corr or list of Corr
+743            Correlator or list of correlators which are plotted for comparison.
+744            The tags of these correlators are used as labels if available.
+745        logscale : bool
+746            Sets y-axis to logscale
+747        plateau : Obs
+748            Plateau value to be visualized in the figure
+749        fit_res : Fit_result
+750            Fit_result object to be visualized
+751        ylabel : str
+752            Label for the y-axis
+753        save : str
+754            path to file in which the figure should be saved
+755        auto_gamma : bool
+756            Apply the gamma method with standard parameters to all correlators and plateau values before plotting.
+757        hide_sigma : float
+758            Hides data points from the first value on which is consistent with zero within 'hide_sigma' standard errors.
+759        references : list
+760            List of floating point values that are displayed as horizontal lines for reference.
+761        """
+762        if self.N != 1:
+763            raise Exception("Correlator must be projected before plotting")
+764
+765        if auto_gamma:
+766            self.gamma_method()
+767
+768        if x_range is None:
+769            x_range = [0, self.T - 1]
+770
+771        fig = plt.figure()
+772        ax1 = fig.add_subplot(111)
+773
+774        x, y, y_err = self.plottable()
+775        if hide_sigma:
+776            hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+777        else:
+778            hide_from = None
+779        ax1.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=self.tag)
+780        if logscale:
+781            ax1.set_yscale('log')
+782        else:
+783            if y_range is None:
+784                try:
+785                    y_min = min([(x[0].value - x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+786                    y_max = max([(x[0].value + x[0].dvalue) for x in self.content[x_range[0]: x_range[1] + 1] if (x is not None) and x[0].dvalue < 2 * np.abs(x[0].value)])
+787                    ax1.set_ylim([y_min - 0.1 * (y_max - y_min), y_max + 0.1 * (y_max - y_min)])
+788                except Exception:
+789                    pass
+790            else:
+791                ax1.set_ylim(y_range)
+792        if comp:
+793            if isinstance(comp, (Corr, list)):
+794                for corr in comp if isinstance(comp, list) else [comp]:
+795                    if auto_gamma:
+796                        corr.gamma_method()
+797                    x, y, y_err = corr.plottable()
+798                    if hide_sigma:
+799                        hide_from = np.argmax((hide_sigma * np.array(y_err[1:])) > np.abs(y[1:])) - 1
+800                    else:
+801                        hide_from = None
+802                    plt.errorbar(x[:hide_from], y[:hide_from], y_err[:hide_from], label=corr.tag, mfc=plt.rcParams['axes.facecolor'])
+803            else:
+804                raise Exception("'comp' must be a correlator or a list of correlators.")
+805
+806        if plateau:
+807            if isinstance(plateau, Obs):
+808                if auto_gamma:
+809                    plateau.gamma_method()
+810                ax1.axhline(y=plateau.value, linewidth=2, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--', label=str(plateau))
+811                ax1.axhspan(plateau.value - plateau.dvalue, plateau.value + plateau.dvalue, alpha=0.25, color=plt.rcParams['text.color'], ls='-')
+812            else:
+813                raise Exception("'plateau' must be an Obs")
+814
+815        if references:
+816            if isinstance(references, list):
+817                for ref in references:
+818                    ax1.axhline(y=ref, linewidth=1, color=plt.rcParams['text.color'], alpha=0.6, marker=',', ls='--')
+819            else:
+820                raise Exception("'references' must be a list of floating pint values.")
+821
+822        if self.prange:
+823            ax1.axvline(self.prange[0], 0, 1, ls='-', marker=',')
+824            ax1.axvline(self.prange[1], 0, 1, ls='-', marker=',')
+825
+826        if fit_res:
+827            x_samples = np.arange(x_range[0], x_range[1] + 1, 0.05)
+828            ax1.plot(x_samples,
+829                     fit_res.fit_function([o.value for o in fit_res.fit_parameters], x_samples),
+830                     ls='-', marker=',', lw=2)
+831
+832        ax1.set_xlabel(r'$x_0 / a$')
+833        if ylabel:
+834            ax1.set_ylabel(ylabel)
+835        ax1.set_xlim([x_range[0] - 0.5, x_range[1] + 0.5])
+836
+837        handles, labels = ax1.get_legend_handles_labels()
+838        if labels:
+839            ax1.legend()
+840        plt.draw()
+841
+842        if save:
+843            if isinstance(save, str):
+844                fig.savefig(save)
+845            else:
+846                raise Exception("'save' has to be a string.")
 

 
 
@@ -4109,34 +4156,34 @@ List of floating point values that are displayed as horizontal lines for referen
 
     

     
-            
834    def spaghetti_plot(self, logscale=True):
-835        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
-836
-837        Parameters
-838        ----------
-839        logscale : bool
-840            Determines whether the scale of the y-axis is logarithmic or standard.
-841        """
-842        if self.N != 1:
-843            raise Exception("Correlator needs to be projected first.")
-844
-845        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]))
-846        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
-847
-848        for name in mc_names:
-849            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+            
848    def spaghetti_plot(self, logscale=True):
+849        """Produces a spaghetti plot of the correlator suited to monitor exceptional configurations.
 850
-851            fig = plt.figure()
-852            ax = fig.add_subplot(111)
-853            for dat in data:
-854                ax.plot(x0_vals, dat, ls='-', marker='')
-855
-856            if logscale is True:
-857                ax.set_yscale('log')
+851        Parameters
+852        ----------
+853        logscale : bool
+854            Determines whether the scale of the y-axis is logarithmic or standard.
+855        """
+856        if self.N != 1:
+857            raise Exception("Correlator needs to be projected first.")
 858
-859            ax.set_xlabel(r'$x_0 / a$')
-860            plt.title(name)
-861            plt.draw()
+859        mc_names = list(set([item for sublist in [o[0].mc_names for o in self.content if o is not None] for item in sublist]))
+860        x0_vals = [n for (n, o) in zip(np.arange(self.T), self.content) if o is not None]
+861
+862        for name in mc_names:
+863            data = np.array([o[0].deltas[name] + o[0].r_values[name] for o in self.content if o is not None]).T
+864
+865            fig = plt.figure()
+866            ax = fig.add_subplot(111)
+867            for dat in data:
+868                ax.plot(x0_vals, dat, ls='-', marker='')
+869
+870            if logscale is True:
+871                ax.set_yscale('log')
+872
+873            ax.set_xlabel(r'$x_0 / a$')
+874            plt.title(name)
+875            plt.draw()
 

 
 
@@ -4163,29 +4210,29 @@ Determines whether the scale of the y-axis is logarithmic or standard.
 
     

     
-            
863    def dump(self, filename, datatype="json.gz", **kwargs):
-864        """Dumps the Corr into a file of chosen type
-865        Parameters
-866        ----------
-867        filename : str
-868            Name of the file to be saved.
-869        datatype : str
-870            Format of the exported file. Supported formats include
-871            "json.gz" and "pickle"
-872        path : str
-873            specifies a custom path for the file (default '.')
-874        """
-875        if datatype == "json.gz":
-876            from .input.json import dump_to_json
-877            if 'path' in kwargs:
-878                file_name = kwargs.get('path') + '/' + filename
-879            else:
-880                file_name = filename
-881            dump_to_json(self, file_name)
-882        elif datatype == "pickle":
-883            dump_object(self, filename, **kwargs)
-884        else:
-885            raise Exception("Unknown datatype " + str(datatype))
+            
877    def dump(self, filename, datatype="json.gz", **kwargs):
+878        """Dumps the Corr into a file of chosen type
+879        Parameters
+880        ----------
+881        filename : str
+882            Name of the file to be saved.
+883        datatype : str
+884            Format of the exported file. Supported formats include
+885            "json.gz" and "pickle"
+886        path : str
+887            specifies a custom path for the file (default '.')
+888        """
+889        if datatype == "json.gz":
+890            from .input.json import dump_to_json
+891            if 'path' in kwargs:
+892                file_name = kwargs.get('path') + '/' + filename
+893            else:
+894                file_name = filename
+895            dump_to_json(self, file_name)
+896        elif datatype == "pickle":
+897            dump_object(self, filename, **kwargs)
+898        else:
+899            raise Exception("Unknown datatype " + str(datatype))
 

 
 
@@ -4217,8 +4264,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
887    def print(self, range=[0, None]):
-888        print(self.__repr__(range))
+            
901    def print(self, range=[0, None]):
+902        print(self.__repr__(range))
 

 
 
@@ -4236,8 +4283,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1050    def sqrt(self):
-1051        return self**0.5
+            
1064    def sqrt(self):
+1065        return self ** 0.5
 

 
 
@@ -4255,9 +4302,9 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1053    def log(self):
-1054        newcontent = [None if (item is None) else np.log(item) for item in self.content]
-1055        return Corr(newcontent, prange=self.prange)
+            
1067    def log(self):
+1068        newcontent = [None if _check_for_none(self, item) else np.log(item) for item in self.content]
+1069        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4275,9 +4322,9 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1057    def exp(self):
-1058        newcontent = [None if (item is None) else np.exp(item) for item in self.content]
-1059        return Corr(newcontent, prange=self.prange)
+            
1071    def exp(self):
+1072        newcontent = [None if _check_for_none(self, item) else np.exp(item) for item in self.content]
+1073        return Corr(newcontent, prange=self.prange)
 

 
 
@@ -4295,8 +4342,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1072    def sin(self):
-1073        return self._apply_func_to_corr(np.sin)
+            
1086    def sin(self):
+1087        return self._apply_func_to_corr(np.sin)
 

 
 
@@ -4314,8 +4361,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1075    def cos(self):
-1076        return self._apply_func_to_corr(np.cos)
+            
1089    def cos(self):
+1090        return self._apply_func_to_corr(np.cos)
 

 
 
@@ -4333,8 +4380,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1078    def tan(self):
-1079        return self._apply_func_to_corr(np.tan)
+            
1092    def tan(self):
+1093        return self._apply_func_to_corr(np.tan)
 

 
 
@@ -4352,8 +4399,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1081    def sinh(self):
-1082        return self._apply_func_to_corr(np.sinh)
+            
1095    def sinh(self):
+1096        return self._apply_func_to_corr(np.sinh)
 

 
 
@@ -4371,8 +4418,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1084    def cosh(self):
-1085        return self._apply_func_to_corr(np.cosh)
+            
1098    def cosh(self):
+1099        return self._apply_func_to_corr(np.cosh)
 

 
 
@@ -4390,8 +4437,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1087    def tanh(self):
-1088        return self._apply_func_to_corr(np.tanh)
+            
1101    def tanh(self):
+1102        return self._apply_func_to_corr(np.tanh)
 

 
 
@@ -4409,8 +4456,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1090    def arcsin(self):
-1091        return self._apply_func_to_corr(np.arcsin)
+            
1104    def arcsin(self):
+1105        return self._apply_func_to_corr(np.arcsin)
 

 
 
@@ -4428,8 +4475,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1093    def arccos(self):
-1094        return self._apply_func_to_corr(np.arccos)
+            
1107    def arccos(self):
+1108        return self._apply_func_to_corr(np.arccos)
 

 
 
@@ -4447,8 +4494,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1096    def arctan(self):
-1097        return self._apply_func_to_corr(np.arctan)
+            
1110    def arctan(self):
+1111        return self._apply_func_to_corr(np.arctan)
 

 
 
@@ -4466,8 +4513,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1099    def arcsinh(self):
-1100        return self._apply_func_to_corr(np.arcsinh)
+            
1113    def arcsinh(self):
+1114        return self._apply_func_to_corr(np.arcsinh)
 

 
 
@@ -4485,8 +4532,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1102    def arccosh(self):
-1103        return self._apply_func_to_corr(np.arccosh)
+            
1116    def arccosh(self):
+1117        return self._apply_func_to_corr(np.arccosh)
 

 
 
@@ -4504,8 +4551,8 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1105    def arctanh(self):
-1106        return self._apply_func_to_corr(np.arctanh)
+            
1119    def arctanh(self):
+1120        return self._apply_func_to_corr(np.arctanh)
 

 
 
@@ -4545,62 +4592,62 @@ specifies a custom path for the file (default '.')
 
     

     
-            
1141    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
-1142        r''' Project large correlation matrix to lowest states
-1143
-1144        This method can be used to reduce the size of an (N x N) correlation matrix
-1145        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
-1146        is still small.
-1147
-1148        Parameters
-1149        ----------
-1150        Ntrunc: int
-1151            Rank of the target matrix.
-1152        tproj: int
-1153            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
-1154            The default value is 3.
-1155        t0proj: int
-1156            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
-1157            discouraged for O(a) improved theories, since the correctness of the procedure
-1158            cannot be granted in this case. The default value is 2.
-1159        basematrix : Corr
-1160            Correlation matrix that is used to determine the eigenvectors of the
-1161            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
-1162            is is not specified.
-1163
-1164        Notes
-1165        -----
-1166        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
-1167        the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
-1168        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
-1169        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
-1170        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
-1171        correlation matrix and to remove some noise that is added by irrelevant operators.
-1172        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
-1173        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
-1174        '''
-1175
-1176        if self.N == 1:
-1177            raise Exception('Method cannot be applied to one-dimensional correlators.')
-1178        if basematrix is None:
-1179            basematrix = self
-1180        if Ntrunc >= basematrix.N:
-1181            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
-1182        if basematrix.N != self.N:
-1183            raise Exception('basematrix and targetmatrix have to be of the same size.')
-1184
-1185        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
-1186
-1187        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
-1188        rmat = []
-1189        for t in range(basematrix.T):
-1190            for i in range(Ntrunc):
-1191                for j in range(Ntrunc):
-1192                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
-1193            rmat.append(np.copy(tmpmat))
-1194
-1195        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
-1196        return Corr(newcontent)
+            
1155    def prune(self, Ntrunc, tproj=3, t0proj=2, basematrix=None):
+1156        r''' Project large correlation matrix to lowest states
+1157
+1158        This method can be used to reduce the size of an (N x N) correlation matrix
+1159        to (Ntrunc x Ntrunc) by solving a GEVP at very early times where the noise
+1160        is still small.
+1161
+1162        Parameters
+1163        ----------
+1164        Ntrunc: int
+1165            Rank of the target matrix.
+1166        tproj: int
+1167            Time where the eigenvectors are evaluated, corresponds to ts in the GEVP method.
+1168            The default value is 3.
+1169        t0proj: int
+1170            Time where the correlation matrix is inverted. Choosing t0proj=1 is strongly
+1171            discouraged for O(a) improved theories, since the correctness of the procedure
+1172            cannot be granted in this case. The default value is 2.
+1173        basematrix : Corr
+1174            Correlation matrix that is used to determine the eigenvectors of the
+1175            lowest states based on a GEVP. basematrix is taken to be the Corr itself if
+1176            is is not specified.
+1177
+1178        Notes
+1179        -----
+1180        We have the basematrix $C(t)$ and the target matrix $G(t)$. We start by solving
+1181        the GEVP $$C(t) v_n(t, t_0) = \lambda_n(t, t_0) C(t_0) v_n(t, t_0)$$ where $t \equiv t_\mathrm{proj}$
+1182        and $t_0 \equiv t_{0, \mathrm{proj}}$. The target matrix is projected onto the subspace of the
+1183        resulting eigenvectors $v_n, n=1,\dots,N_\mathrm{trunc}$ via
+1184        $$G^\prime_{i, j}(t) = (v_i, G(t) v_j)$$. This allows to reduce the size of a large
+1185        correlation matrix and to remove some noise that is added by irrelevant operators.
+1186        This may allow to use the GEVP on $G(t)$ at late times such that the theoretically motivated
+1187        bound $t_0 \leq t/2$ holds, since the condition number of $G(t)$ is decreased, compared to $C(t)$.
+1188        '''
+1189
+1190        if self.N == 1:
+1191            raise Exception('Method cannot be applied to one-dimensional correlators.')
+1192        if basematrix is None:
+1193            basematrix = self
+1194        if Ntrunc >= basematrix.N:
+1195            raise Exception('Cannot truncate using Ntrunc <= %d' % (basematrix.N))
+1196        if basematrix.N != self.N:
+1197            raise Exception('basematrix and targetmatrix have to be of the same size.')
+1198
+1199        evecs = basematrix.GEVP(t0proj, tproj, sort=None)[:Ntrunc]
+1200
+1201        tmpmat = np.empty((Ntrunc, Ntrunc), dtype=object)
+1202        rmat = []
+1203        for t in range(basematrix.T):
+1204            for i in range(Ntrunc):
+1205                for j in range(Ntrunc):
+1206                    tmpmat[i][j] = evecs[i].T @ self[t] @ evecs[j]
+1207            rmat.append(np.copy(tmpmat))
+1208
+1209        newcontent = [None if (self.content[t] is None) else rmat[t] for t in range(self.T)]
+1210        return Corr(newcontent)