From 344da7a3d9e290bcdececadcac3dcf077ebbebca Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 26 Feb 2022 08:45:33 +0000 Subject: [PATCH 1/5] feat: first version of strictly positive semi-definite covariance --- pyerrors/obs.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index aea75752..8089cf2e 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -1390,8 +1390,6 @@ def covariance(obs1, obs2, correlation=False, **kwargs): dvalue = 0 e_gamma = {} e_dvalue = {} - e_n_tauint = {} - e_rho = {} for e_name in obs1.mc_names: @@ -1435,14 +1433,10 @@ def covariance(obs1, obs2, correlation=False, **kwargs): e_N += np.sum(np.ones_like(idl_d[r_name])) e_gamma[e_name] /= gamma_div[:w_max] - e_rho[e_name] = e_gamma[e_name][:w_max] / e_gamma[e_name][0] - e_n_tauint[e_name] = np.cumsum(np.concatenate(([0.5], e_rho[e_name][1:]))) - # Make sure no entry of tauint is smaller than 0.5 - e_n_tauint[e_name][e_n_tauint[e_name] < 0.5] = 0.500000000001 + tau_int = min(obs1.e_tauint[e_name], obs2.e_tauint[e_name]) - window = min(obs1.e_windowsize[e_name], obs2.e_windowsize[e_name]) # Bias correction hep-lat/0306017 eq. (49) - e_dvalue[e_name] = 2 * (e_n_tauint[e_name][window] + obs1.tau_exp[e_name] * np.abs(e_rho[e_name][window + 1])) * (1 + (2 * window + 1) / e_N) * e_gamma[e_name][0] / e_N + e_dvalue[e_name] = 2 * (tau_int) * (1 + 1 / e_N) * e_gamma[e_name][0] / e_N dvalue += e_dvalue[e_name] From 2ba59f90c0619830b92cfb019b5b1aa581688633 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 28 Feb 2022 13:18:04 +0000 Subject: [PATCH 2/5] Version number bumped to 2.1.0+dev --- pyerrors/version.py | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyerrors/version.py b/pyerrors/version.py index c71c5a3d..6cb8ec9d 100644 --- a/pyerrors/version.py +++ b/pyerrors/version.py @@ -1 +1 @@ -__version__ = "2.0.0-rc.2+dev" +__version__ = "2.1.0+dev" diff --git a/setup.py b/setup.py index 0edc33a9..b578a727 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,7 @@ from setuptools import setup, find_packages setup(name='pyerrors', - version='2.0.0-rc.2+dev', + version='2.1.0+dev', description='Error analysis for lattice QCD', author='Fabian Joswig', author_email='fabian.joswig@ed.ac.uk', From 0f749fd107a4cddb665f765170f4740ad10e5bae Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 28 Feb 2022 13:25:09 +0000 Subject: [PATCH 3/5] refactor: instantiation of Obs in import_jackknife slightly optimized --- pyerrors/obs.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 8089cf2e..9ad3a0d4 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -1470,7 +1470,8 @@ def import_jackknife(jacks, name, idl=None): length = len(jacks) - 1 prj = (np.ones((length, length)) - (length - 1) * np.identity(length)) samples = jacks[1:] @ prj - new_obs = Obs([samples], [name], idl=idl) + mean = np.mean(samples) + new_obs = Obs([samples - mean], [name], idl=idl, means=[mean]) new_obs._value = jacks[0] return new_obs From 42df2542885d980720f92647b8dd93cd68b7ce3d Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 28 Feb 2022 13:26:48 +0000 Subject: [PATCH 4/5] refactor: else case for empty observables in Obs.__init__ simplified. --- pyerrors/obs.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 9ad3a0d4..cdf7bce8 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -90,6 +90,9 @@ class Obs: self.deltas = {} self._covobs = {} + self._value = 0 + self.N = 0 + self.is_merged = {} self.idl = {} if len(samples): if idl is not None: @@ -110,8 +113,6 @@ class Obs: for name, sample in sorted(zip(names, samples)): self.idl[name] = range(1, len(sample) + 1) - self._value = 0 - self.N = 0 if kwargs.get("means") is not None: for name, sample, mean in sorted(zip(names, samples, kwargs.get("means"))): self.shape[name] = len(self.idl[name]) @@ -129,13 +130,6 @@ class Obs: self._value += self.shape[name] * self.r_values[name] self._value /= self.N - self.is_merged = {} - - else: - self._value = 0 - self.is_merged = {} - self.N = 0 - self._dvalue = 0.0 self.ddvalue = 0.0 self.reweighted = False From 498a251072afaafd8f14c8ef4d7fa7e3b60268e1 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 28 Feb 2022 13:43:49 +0000 Subject: [PATCH 5/5] refactor!: if clause in Obs.__init__ eliminated, empty observables need to be initialized with means=[] from now on. --- pyerrors/input/json.py | 6 ++-- pyerrors/obs.py | 67 +++++++++++++++++++++--------------------- tests/obs_test.py | 4 +-- 3 files changed, 38 insertions(+), 39 deletions(-) diff --git a/pyerrors/input/json.py b/pyerrors/input/json.py index a6060d5f..c4c2aa4b 100644 --- a/pyerrors/input/json.py +++ b/pyerrors/input/json.py @@ -358,7 +358,7 @@ def _parse_json_dict(json_dict, verbose=True, full_output=False): ret = Obs([[ddi[0] + values[0] for ddi in di] for di in od['deltas']], od['names'], idl=od['idl']) ret.is_merged = od['is_merged'] else: - ret = Obs([], []) + ret = Obs([], [], means=[]) ret._value = values[0] for name in cd: co = cd[name][0] @@ -383,7 +383,7 @@ def _parse_json_dict(json_dict, verbose=True, full_output=False): ret.append(Obs([list(di[:, i] + values[i]) for di in od['deltas']], od['names'], idl=od['idl'])) ret[-1].is_merged = od['is_merged'] else: - ret.append(Obs([], [])) + ret.append(Obs([], [], means=[])) ret[-1]._value = values[i] print('Created Obs with means= ', values[i]) for name in cd: @@ -410,7 +410,7 @@ def _parse_json_dict(json_dict, verbose=True, full_output=False): ret.append(Obs([di[:, i] + values[i] for di in od['deltas']], od['names'], idl=od['idl'])) ret[-1].is_merged = od['is_merged'] else: - ret.append(Obs([], [])) + ret.append(Obs([], [], means=[])) ret[-1]._value = values[i] for name in cd: co = cd[name][i] diff --git a/pyerrors/obs.py b/pyerrors/obs.py index cdf7bce8..1d459247 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -94,41 +94,40 @@ class Obs: self.N = 0 self.is_merged = {} self.idl = {} - if len(samples): - if idl is not None: - for name, idx in sorted(zip(names, idl)): - if isinstance(idx, range): - self.idl[name] = idx - elif isinstance(idx, (list, np.ndarray)): - dc = np.unique(np.diff(idx)) - if np.any(dc < 0): - raise Exception("Unsorted idx for idl[%s]" % (name)) - if len(dc) == 1: - self.idl[name] = range(idx[0], idx[-1] + dc[0], dc[0]) - else: - self.idl[name] = list(idx) + if idl is not None: + for name, idx in sorted(zip(names, idl)): + if isinstance(idx, range): + self.idl[name] = idx + elif isinstance(idx, (list, np.ndarray)): + dc = np.unique(np.diff(idx)) + if np.any(dc < 0): + raise Exception("Unsorted idx for idl[%s]" % (name)) + if len(dc) == 1: + self.idl[name] = range(idx[0], idx[-1] + dc[0], dc[0]) else: - raise Exception('incompatible type for idl[%s].' % (name)) - else: - for name, sample in sorted(zip(names, samples)): - self.idl[name] = range(1, len(sample) + 1) + self.idl[name] = list(idx) + else: + raise Exception('incompatible type for idl[%s].' % (name)) + else: + for name, sample in sorted(zip(names, samples)): + self.idl[name] = range(1, len(sample) + 1) - if kwargs.get("means") is not None: - for name, sample, mean in sorted(zip(names, samples, kwargs.get("means"))): - self.shape[name] = len(self.idl[name]) - self.N += self.shape[name] - self.r_values[name] = mean - self.deltas[name] = sample - else: - for name, sample in sorted(zip(names, samples)): - self.shape[name] = len(self.idl[name]) - self.N += self.shape[name] - if len(sample) != self.shape[name]: - raise Exception('Incompatible samples and idx for %s: %d vs. %d' % (name, len(sample), self.shape[name])) - self.r_values[name] = np.mean(sample) - self.deltas[name] = sample - self.r_values[name] - self._value += self.shape[name] * self.r_values[name] - self._value /= self.N + if kwargs.get("means") is not None: + for name, sample, mean in sorted(zip(names, samples, kwargs.get("means"))): + self.shape[name] = len(self.idl[name]) + self.N += self.shape[name] + self.r_values[name] = mean + self.deltas[name] = sample + else: + for name, sample in sorted(zip(names, samples)): + self.shape[name] = len(self.idl[name]) + self.N += self.shape[name] + if len(sample) != self.shape[name]: + raise Exception('Incompatible samples and idx for %s: %d vs. %d' % (name, len(sample), self.shape[name])) + self.r_values[name] = np.mean(sample) + self.deltas[name] = sample - self.r_values[name] + self._value += self.shape[name] * self.r_values[name] + self._value /= self.N self._dvalue = 0.0 self.ddvalue = 0.0 @@ -1522,7 +1521,7 @@ def cov_Obs(means, cov, name, grad=None): co : Covobs Covobs to be embedded into the Obs """ - o = Obs([], []) + o = Obs([], [], means=[]) o._value = co.value o.names.append(co.name) o._covobs[co.name] = co diff --git a/tests/obs_test.py b/tests/obs_test.py index 2fed00af..14f18c81 100644 --- a/tests/obs_test.py +++ b/tests/obs_test.py @@ -653,7 +653,7 @@ def test_covariance_symmetry(): def test_empty_obs(): o = pe.Obs([np.random.rand(100)], ['test']) - q = o + pe.Obs([], []) + q = o + pe.Obs([], [], means=[]) assert q == o @@ -704,4 +704,4 @@ def test_merge_idx(): new_idx = pe.obs._merge_idx(idl) assert(new_idx[-1] > new_idx[0]) for i in range(1, len(new_idx)): - assert(new_idx[i - 1] < new_idx[i]) \ No newline at end of file + assert(new_idx[i - 1] < new_idx[i])