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)
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) +@@ -3024,13 +3061,13 @@ timeslice and the error on each timeslice.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)
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.") +@@ -3050,87 +3087,87 @@ timeslice and the error on each timeslice.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.")
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") +@@ -3181,18 +3218,18 @@ Returns only the vector(s) for a specified state. The lowest state is zero.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
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) +@@ -3220,46 +3257,46 @@ The state one is interested in ordered by energy. The lowest state is zero.331 def Eigenvalue(self, t0, ts=None, state=0, sort="Eigenvalue"): +332 """Determines the eigenvalue of the GEVP by solving and projecting the correlator +333 +334 Parameters +335 ---------- +336 state : int +337 The state one is interested in ordered by energy. The lowest state is zero. +338 +339 All other parameters are identical to the ones of Corr.GEVP. +340 """ +341 vec = self.GEVP(t0, ts=ts, sort=sort)[state] +342 return self.projected(vec)
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) +@@ -3293,15 +3330,15 @@ determines whether the matrix is extended periodically344 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)
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))) + @@ -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]) + @@ -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) +@@ -3395,32 +3432,34 @@ Offset the equal spacing399 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)
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) +@@ -3449,26 +3488,28 @@ correlator or a Corr of same length.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)
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) +@@ -3500,33 +3541,35 @@ on the configurations in obs[i].idl.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)
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) +@@ -3555,57 +3598,59 @@ Parity quantum number of the correlator, can be +1 or -1469 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
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.") +@@ -3633,37 +3678,39 @@ Available choice: symmetric, forward, backward, improved, default: symmetric -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.")-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.") +@@ -3691,70 +3738,70 @@ Available choice: symmetric, improved, default: symmetric553 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.")-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]) +@@ -3787,39 +3834,39 @@ guess for the root finder, only relevant for the root variant587 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.')-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 +@@ -3853,42 +3900,42 @@ Decides whether output is printed to the standard output.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-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) +@@ -3922,17 +3969,17 @@ apply gamma_method with default parameters to the Corr. Defaults to None686 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)-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 +@@ -3952,118 +3999,118 @@ apply gamma_method with default parameters to the Corr. Defaults to None723 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-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.") +@@ -4109,34 +4156,34 @@ List of floating point values that are displayed as horizontal lines for referen735 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.")-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 +@@ -4163,29 +4210,29 @@ Determines whether the scale of the y-axis is logarithmic or standard.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()-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)) +@@ -4217,8 +4264,8 @@ specifies a custom path for the file (default '.')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))-887 def print(self, range=[0, None]): -888 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 + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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) + @@ -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)