From fd4c866fddaeff295edb67a8a62f3a34fb759e07 Mon Sep 17 00:00:00 2001 From: Simon Kuberski Date: Mon, 17 Oct 2022 14:38:23 +0200 Subject: [PATCH 1/3] fix: bug in automatic window for irregular chains fixed --- pyerrors/obs.py | 19 +++++++++++++------ tests/obs_test.py | 10 ++++++++++ 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 24ef3720..481b71a8 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -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) diff --git a/tests/obs_test.py b/tests/obs_test.py index d58933c9..d996e065 100644 --- a/tests/obs_test.py +++ b/tests/obs_test.py @@ -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) From 89b0c37e7604a640e3373252a092e8a2c0a52b64 Mon Sep 17 00:00:00 2001 From: Simon Kuberski Date: Mon, 17 Oct 2022 17:40:27 +0200 Subject: [PATCH 2/3] fix: corrected bias correction of tau_int for irregular chains --- pyerrors/obs.py | 20 +++++++++++++------- tests/obs_test.py | 19 +++++++++++++++++++ 2 files changed, 32 insertions(+), 7 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 481b71a8..03486f14 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -279,11 +279,17 @@ 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) - # 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])) + 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: @@ -295,7 +301,7 @@ class Obs: _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) @@ -318,7 +324,7 @@ class Obs: _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_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) diff --git a/tests/obs_test.py b/tests/obs_test.py index d996e065..29ee67ba 100644 --- a/tests/obs_test.py +++ b/tests/obs_test.py @@ -649,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) From 150a13136c947ea3c8a2bf52b695c56254350776 Mon Sep 17 00:00:00 2001 From: Simon Kuberski Date: Mon, 17 Oct 2022 17:50:24 +0200 Subject: [PATCH 3/3] fix: corrected bias correction of ddvalue for irregular chains --- pyerrors/obs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 03486f14..c8373d5e 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -305,7 +305,7 @@ class Obs: 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: @@ -327,7 +327,7 @@ class Obs: 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