From ec20ee38a66b7be8c5dc544f8944615329089076 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 13 Dec 2021 17:06:03 +0000 Subject: [PATCH] feat!: covariance replaced by covariance2, window altered to minimum of the window of the two observables. Tests adjusted. --- pyerrors/obs.py | 72 +------------------------------------------- tests/covobs_test.py | 2 +- tests/fits_test.py | 3 +- tests/obs_test.py | 12 ++++---- 4 files changed, 10 insertions(+), 79 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 34bbf46d..14138b92 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -1334,76 +1334,6 @@ def covariance(obs1, obs2, correlation=False, **kwargs): is constrained to the maximum value in order to make sure that covariance matrices are positive semidefinite. - Parameters - ---------- - obs1 : Obs - First Obs - obs2 : Obs - Second Obs - correlation : bool - if true the correlation instead of the covariance is - returned (default False) - """ - if set(obs1.names).isdisjoint(set(obs2.names)): - return 0. - - for name in sorted(set(obs1.names + obs2.names)): - if (obs1.shape.get(name) != obs2.shape.get(name)) and (obs1.shape.get(name) is not None) and (obs2.shape.get(name) is not None): - raise Exception('Shapes of ensemble', name, 'do not fit') - if (1 != len(set([len(idx) for idx in [obs1.idl[name], obs2.idl[name], _merge_idx([obs1.idl[name], obs2.idl[name]])]]))): - raise Exception('Shapes of ensemble', name, 'do not fit') - - if not hasattr(obs1, 'e_dvalue') or not hasattr(obs2, 'e_dvalue'): - raise Exception('The gamma method has to be applied to both Obs first.') - - dvalue = 0 - - for e_name in obs1.mc_names: - - if e_name not in obs2.e_names: - continue - - gamma = 0 - r_length = [] - for r_name in obs1.e_content[e_name]: - if r_name not in obs2.e_content[e_name]: - continue - - r_length.append(len(obs1.deltas[r_name])) - - gamma += np.sum(obs1.deltas[r_name] * obs2.deltas[r_name]) - - e_N = np.sum(r_length) - - tau_combined = (obs1.e_tauint[e_name] + obs2.e_tauint[e_name]) / 2 - dvalue += gamma / e_N * (1 + 1 / e_N) / e_N * 2 * tau_combined - - for e_name in obs1.cov_names: - - if e_name not in obs2.cov_names: - continue - - dvalue += float(np.dot(np.transpose(obs1.covobs[e_name].grad), np.dot(obs1.covobs[e_name].cov, obs2.covobs[e_name].grad))) - - if np.abs(dvalue / obs1.dvalue / obs2.dvalue) > 1.0: - dvalue = np.sign(dvalue) * obs1.dvalue * obs2.dvalue - - if correlation: - dvalue = dvalue / obs1.dvalue / obs2.dvalue - - return dvalue - - -def covariance2(obs1, obs2, correlation=False, **kwargs): - """Alternative implementation of the covariance of two observables. - - covariance(obs, obs) is equal to obs.dvalue ** 2 - The gamma method has to be applied first to both observables. - - If abs(covariance(obs1, obs2)) > obs1.dvalue * obs2.dvalue, the covariance - is constrained to the maximum value in order to make sure that covariance - matrices are positive semidefinite. - Keyword arguments ----------------- correlation -- if true the correlation instead of the covariance is @@ -1503,7 +1433,7 @@ def covariance2(obs1, obs2, correlation=False, **kwargs): # 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 - window = max(obs1.e_windowsize[e_name], obs2.e_windowsize[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 diff --git a/tests/covobs_test.py b/tests/covobs_test.py index 513fb351..215c75d3 100644 --- a/tests/covobs_test.py +++ b/tests/covobs_test.py @@ -40,7 +40,7 @@ def test_covobs(): [o.gamma_method() for o in cl] assert(pe.covariance(cl[0], cl[1]) == cov[0][1]) - assert(pe.covariance2(cl[0], cl[1]) == cov[1][0]) + assert(pe.covariance(cl[0], cl[1]) == cov[1][0]) do = cl[0] * cl[1] assert(np.array_equal(do.covobs['rAP'].grad, np.transpose([pi[1], pi[0]]).reshape(2, 1))) diff --git a/tests/fits_test.py b/tests/fits_test.py index b37abc25..8a9e0843 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -83,6 +83,8 @@ def test_least_squares(): assert math.isclose(pcov[i, i], betac[i].dvalue ** 2, abs_tol=1e-3) assert math.isclose(pe.covariance(betac[0], betac[1]), pcov[0, 1], abs_tol=1e-3) + +def test_correlated_fit(): num_samples = 400 N = 10 @@ -101,7 +103,6 @@ def test_least_squares(): c = cholesky(r, lower=True) y = np.dot(c, x) - x = np.arange(N) for linear in [True, False]: data = [] diff --git a/tests/obs_test.py b/tests/obs_test.py index fe8826f8..989a0064 100644 --- a/tests/obs_test.py +++ b/tests/obs_test.py @@ -555,7 +555,7 @@ def test_gamma_method_irregular(): assert((ae.e_tauint['a'] + 4 * ae.e_dtauint['a'] > ao.e_tauint['a'])) -def test_covariance2_symmetry(): +def test_covariance_symmetry(): value1 = np.random.normal(5, 10) dvalue1 = np.abs(np.random.normal(0, 1)) test_obs1 = pe.pseudo_Obs(value1, dvalue1, 't') @@ -564,8 +564,8 @@ def test_covariance2_symmetry(): dvalue2 = np.abs(np.random.normal(0, 1)) test_obs2 = pe.pseudo_Obs(value2, dvalue2, 't') test_obs2.gamma_method() - cov_ab = pe.covariance2(test_obs1, test_obs2) - cov_ba = pe.covariance2(test_obs2, test_obs1) + cov_ab = pe.covariance(test_obs1, test_obs2) + cov_ba = pe.covariance(test_obs2, test_obs1) assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float64).eps assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps) @@ -578,10 +578,10 @@ def test_covariance2_symmetry(): idx = [i + 1 for i in range(len(configs)) if configs[i] == 1] a = pe.Obs([zero_arr], ['t'], idl=[idx]) a.gamma_method() - assert np.isclose(a.dvalue**2, pe.covariance2(a, a), atol=100, rtol=1e-4) + assert np.isclose(a.dvalue**2, pe.covariance(a, a), atol=100, rtol=1e-4) - cov_ab = pe.covariance2(test_obs1, a) - cov_ba = pe.covariance2(a, test_obs1) + cov_ab = pe.covariance(test_obs1, a) + cov_ba = pe.covariance(a, test_obs1) assert np.abs(cov_ab - cov_ba) <= 10 * np.finfo(np.float64).eps assert np.abs(cov_ab) < test_obs1.dvalue * test_obs2.dvalue * (1 + 10 * np.finfo(np.float64).eps)