lattice/pyerrors
1
0
Fork 0
mirror of https://github.com/fjosw/pyerrors.git synced 2025-05-14 19:43:41 +02:00
Code Issues Projects Releases Packages Wiki Activity

fix: bug in automatic window for irregular chains fixed

Browse source
This commit is contained in:
Simon Kuberski 2022-10-17 14:38:23 +02:00
parent 5e9cc3a807
commit fd4c866fdd
2 changed files with 23 additions and 6 deletions
Download patch file Download diff file Expand all files Collapse all files

19
pyerrors/obs.py
View file

@ -279,14 +279,20 @@ class Obs:
tmp = self.e_rho[e_name][i + 1:w_max] + np.concatenate([self.e_rho[e_name][i - 1::-1], self.e_rho[e_name][1:w_max - 2 * i]]) - 2 * self.e_rho[e_name][i] * self.e_rho[e_name][1:w_max - i]
self.e_drho[e_name][i] = np.sqrt(np.sum(tmp ** 2) / e_N)
_compute_drho(1)
# detect regular step size in irregular MC chain
gapsize = 1
for r_name in e_content[e_name][:1]:
if not isinstance(self.idl[r_name], range):
gapsize = np.min(np.diff(self.idl[r_name]))
_compute_drho(gapsize)
if self.tau_exp[e_name] > 0:
texp = self.tau_exp[e_name]
# Critical slowing down analysis
if w_max // 2 <= 1:
raise Exception("Need at least 8 samples for tau_exp error analysis")
for n in range(1, w_max // 2):
_compute_drho(n + 1)
for n in range(gapsize, w_max // 2, gapsize):
_compute_drho(n + gapsize)
if (self.e_rho[e_name][n] - self.N_sigma[e_name] * self.e_drho[e_name][n]) < 0 or n >= w_max // 2 - 2:
# Bias correction hep-lat/0306017 eq. (49) included
self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n + 1) / e_N) / (1 + 1 / e_N) + texp * np.abs(self.e_rho[e_name][n + 1]) # The absolute makes sure, that the tail contribution is always positive
@ -305,12 +311,13 @@ class Obs:
self.e_windowsize[e_name] = 0
else:
# Standard automatic windowing procedure
tau = self.S[e_name] / np.log((2 * self.e_n_tauint[e_name][1:] + 1) / (2 * self.e_n_tauint[e_name][1:] - 1))
g_w = np.exp(- np.arange(1, w_max) / tau) - tau / np.sqrt(np.arange(1, w_max) * e_N)
tau = self.S[e_name] / np.log((2 * self.e_n_tauint[e_name][gapsize::gapsize] + 1) / (2 * self.e_n_tauint[e_name][gapsize::gapsize] - 1))
g_w = np.exp(- np.arange(1, len(tau) + 1) / tau) - tau / np.sqrt(np.arange(1, len(tau) + 1) * e_N)
for n in range(1, w_max):
if n < w_max // 2 - 2:
_compute_drho(n + 1)
_compute_drho(gapsize * n + gapsize)
if g_w[n - 1] < 0 or n >= w_max - 1:
n *= gapsize
self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n + 1) / e_N) / (1 + 1 / e_N) # Bias correction hep-lat/0306017 eq. (49)
self.e_dtauint[e_name] = self.e_n_dtauint[e_name][n]
self.e_dvalue[e_name] = np.sqrt(2 * self.e_tauint[e_name] * e_gamma[e_name][0] * (1 + 1 / e_N) / e_N)