mirror of
https://github.com/fjosw/pyerrors.git
synced 2025-03-15 06:40:24 +01:00
Merge pull request #124 from s-kuberski/fix/irregular
fix: bug in automatic window for irregular chains fixed
This commit is contained in:
commit
760916692d
2 changed files with 52 additions and 10 deletions
|
@ -279,21 +279,33 @@ 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)
|
||||
gaps = []
|
||||
for r_name in e_content[e_name]:
|
||||
if isinstance(self.idl[r_name], range):
|
||||
gaps.append(1)
|
||||
else:
|
||||
gaps.append(np.min(np.diff(self.idl[r_name])))
|
||||
|
||||
if not np.all([gi == gaps[0] for gi in gaps]):
|
||||
raise Exception(f"Replica for ensemble {e_name} are not equally spaced.", gaps)
|
||||
else:
|
||||
gapsize = gaps[0]
|
||||
|
||||
_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
|
||||
self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n / gapsize + 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
|
||||
self.e_dtauint[e_name] = np.sqrt(self.e_n_dtauint[e_name][n] ** 2 + texp ** 2 * self.e_drho[e_name][n + 1] ** 2)
|
||||
# Error of tau_exp neglected so far, missing term: self.e_rho[e_name][n + 1] ** 2 * d_tau_exp ** 2
|
||||
self.e_dvalue[e_name] = np.sqrt(2 * self.e_tauint[e_name] * e_gamma[e_name][0] * (1 + 1 / e_N) / e_N)
|
||||
self.e_ddvalue[e_name] = self.e_dvalue[e_name] * np.sqrt((n + 0.5) / e_N)
|
||||
self.e_ddvalue[e_name] = self.e_dvalue[e_name] * np.sqrt((n / gapsize + 0.5) / e_N)
|
||||
self.e_windowsize[e_name] = n
|
||||
break
|
||||
else:
|
||||
|
@ -305,16 +317,17 @@ 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:
|
||||
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)
|
||||
n *= gapsize
|
||||
self.e_tauint[e_name] = self.e_n_tauint[e_name][n] * (1 + (2 * n / gapsize + 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)
|
||||
self.e_ddvalue[e_name] = self.e_dvalue[e_name] * np.sqrt((n + 0.5) / e_N)
|
||||
self.e_ddvalue[e_name] = self.e_dvalue[e_name] * np.sqrt((n / gapsize + 0.5) / e_N)
|
||||
self.e_windowsize[e_name] = n
|
||||
break
|
||||
|
||||
|
|
|
@ -628,9 +628,19 @@ def test_gamma_method_irregular():
|
|||
ao = pe.Obs([[carr[i] for i in range(len(carr)) if i % 2 == 1]], ['a'], idl=[[i for i in range(len(carr)) if i % 2 == 1]])
|
||||
ao.gamma_method()
|
||||
|
||||
arrt = [carr[i] for i in range(len(carr)) if i % 2 == 1]
|
||||
idlt = [i for i in range(len(carr)) if i % 2 == 1]
|
||||
for el in [int(e) for e in N * np.random.uniform(size=10)]:
|
||||
arrt = arrt[:el] + arrt[el + 1:]
|
||||
idlt = idlt[:el] + idlt[el + 1:]
|
||||
ai = pe.Obs([arrt], ['a'], idl=[idlt])
|
||||
ai.gamma_method()
|
||||
|
||||
assert(ae.e_tauint['a'] < a.e_tauint['a'])
|
||||
assert((ae.e_tauint['a'] - 4 * ae.e_dtauint['a'] < ao.e_tauint['a']))
|
||||
assert((ae.e_tauint['a'] + 4 * ae.e_dtauint['a'] > ao.e_tauint['a']))
|
||||
assert((ai.e_tauint['a'] - 4 * ai.e_dtauint['a'] < ao.e_tauint['a']))
|
||||
assert((ai.e_tauint['a'] + 4 * ai.e_dtauint['a'] > ao.e_tauint['a']))
|
||||
|
||||
a = pe.pseudo_Obs(1, .1, 'a', samples=10)
|
||||
a.idl['a'] = range(4, 15)
|
||||
|
@ -639,6 +649,25 @@ def test_gamma_method_irregular():
|
|||
ol = [a, b]
|
||||
o = (ol[0] - ol[1]) / (ol[1])
|
||||
|
||||
N = 1000
|
||||
dat = gen_autocorrelated_array(np.random.normal(1, .2, size=N), .8)
|
||||
|
||||
idl_a = list(range(0, 1001, 1))
|
||||
idl_a.remove(101)
|
||||
|
||||
oa = pe.Obs([dat], ["ens1"], idl=[idl_a])
|
||||
oa.gamma_method()
|
||||
tau_a = oa.e_tauint["ens1"]
|
||||
|
||||
idl_b = list(range(0, 10010, 10))
|
||||
idl_b.remove(1010)
|
||||
|
||||
ob = pe.Obs([dat], ["ens1"], idl=[idl_b])
|
||||
ob.gamma_method()
|
||||
tau_b = ob.e_tauint["ens1"]
|
||||
|
||||
assert np.isclose(tau_a, tau_b)
|
||||
|
||||
|
||||
def test_covariance_is_variance():
|
||||
value = np.random.normal(5, 10)
|
||||
|
|
Loading…
Add table
Reference in a new issue