From ce4e73ec1a5f6c8c0edd1200a6a9524f95dfbb2d Mon Sep 17 00:00:00 2001 From: Simon Kuberski Date: Fri, 3 Dec 2021 14:04:27 +0100 Subject: [PATCH 1/9] Reweighting is now possible if the observable is defined only on a subset of replica of those of the RWF --- pyerrors/obs.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index b8bb371e..5570f509 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -1238,22 +1238,22 @@ def reweight(weight, obs, **kwargs): for i in range(len(obs)): if len(obs[i].cov_names): raise Exception('Error: Not possible to reweight an Obs that contains covobs!') - if sorted(weight.names) != sorted(obs[i].names): + if not set(obs[i].names).issubset(weight.names): raise Exception('Error: Ensembles do not fit') - for name in weight.names: + for name in obs[i].names: if not set(obs[i].idl[name]).issubset(weight.idl[name]): raise Exception('obs[%d] has to be defined on a subset of the configs in weight.idl[%s]!' % (i, name)) new_samples = [] w_deltas = {} - for name in sorted(weight.names): + for name in sorted(obs[i].names): w_deltas[name] = _reduce_deltas(weight.deltas[name], weight.idl[name], obs[i].idl[name]) new_samples.append((w_deltas[name] + weight.r_values[name]) * (obs[i].deltas[name] + obs[i].r_values[name])) - tmp_obs = Obs(new_samples, sorted(weight.names), idl=[obs[i].idl[name] for name in sorted(weight.names)]) + tmp_obs = Obs(new_samples, sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)]) if kwargs.get('all_configs'): new_weight = weight else: - new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(weight.names)], sorted(weight.names), idl=[obs[i].idl[name] for name in sorted(weight.names)]) + new_weight = Obs([w_deltas[name] + weight.r_values[name] for name in sorted(obs[i].names)], sorted(obs[i].names), idl=[obs[i].idl[name] for name in sorted(obs[i].names)]) result.append(derived_observable(lambda x, **kwargs: x[0] / x[1], [tmp_obs, new_weight], **kwargs)) result[-1].reweighted = True From 06dc6cd808fdc89820f434181951b03d99ae2ead Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 4 Dec 2021 11:55:11 +0000 Subject: [PATCH 2/9] fix: r_values in fits are now properly propagated --- pyerrors/fits.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 677f6eba..85e28557 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -301,8 +301,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs): result = [] for i in range(n_parms): - result.append(derived_observable(lambda x, **kwargs: x[0], list(x.ravel()) + list(y), man_grad=list(deriv_x[i]) + list(deriv_y[i]))) - result[-1]._value = out.beta[i] + result.append(derived_observable(lambda my_var, **kwargs: my_var[0] / x.ravel()[0].value * out.beta[i], list(x.ravel()) + list(y), man_grad=list(deriv_x[i]) + list(deriv_y[i]))) output.fit_parameters = result + const_par @@ -419,8 +418,7 @@ def _prior_fit(x, y, func, priors, silent=False, **kwargs): result = [] for i in range(n_parms): - result.append(derived_observable(lambda x, **kwargs: x[0], list(y) + list(loc_priors), man_grad=list(deriv[i]))) - result[-1]._value = params[i] + result.append(derived_observable(lambda x, **kwargs: x[0] / y[0].value * params[i], list(y) + list(loc_priors), man_grad=list(deriv[i]))) output.fit_parameters = result output.chisquare = chisqfunc(np.asarray(params)) @@ -614,8 +612,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs): result = [] for i in range(n_parms): - result.append(derived_observable(lambda x, **kwargs: x[0], list(y), man_grad=list(deriv[i]))) - result[-1]._value = fit_result.x[i] + result.append(derived_observable(lambda x, **kwargs: x[0] / y[0].value * fit_result.x[i], list(y), man_grad=list(deriv[i]))) output.fit_parameters = result + const_par From 28f1372cfd569633979abd497076c55704d73773 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 4 Dec 2021 13:02:45 +0000 Subject: [PATCH 3/9] fix: r_value propagation also adjusted in root module --- pyerrors/roots.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyerrors/roots.py b/pyerrors/roots.py index 99434a1c..1cb7b46f 100644 --- a/pyerrors/roots.py +++ b/pyerrors/roots.py @@ -33,6 +33,5 @@ def find_root(d, func, guess=1.0, **kwargs): da = jacobian(lambda u, v: func(v, u))(d.value, root[0]) deriv = - da / dx - res = derived_observable(lambda x, **kwargs: x[0], [d], man_grad=[deriv]) - res._value = root[0] + res = derived_observable(lambda x, **kwargs: x[0] / d.value * root[0], [d], man_grad=[deriv]) return res From 602827434221f4a36c6302a31497137e8eb359f0 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 4 Dec 2021 13:05:18 +0000 Subject: [PATCH 4/9] test: test for r_value propagation in fits added --- tests/fits_test.py | 44 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 41 insertions(+), 3 deletions(-) diff --git a/tests/fits_test.py b/tests/fits_test.py index 45e98e35..5d7f1de3 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -52,7 +52,7 @@ def test_least_squares(): outc = pe.least_squares(x, oyc, func) betac = outc.fit_parameters - + for i in range(2): betac[i].gamma_method(S=1.0) assert math.isclose(betac[i].value, popt[i], abs_tol=1e-5) @@ -97,7 +97,7 @@ def test_least_squares(): return p[1] * np.exp(-p[0] * x) fitp = pe.least_squares(x, data, fitf, expected_chisquare=True) - + fitpc = pe.least_squares(x, data, fitf, correlated_fit=True) for i in range(2): diff = fitp[i] - fitpc[i] @@ -170,7 +170,7 @@ def test_total_least_squares(): diffc = outc.fit_parameters[0] - betac[0] assert(diffc / betac[0] < 1e-3 * betac[0].dvalue) assert((outc.fit_parameters[1] - betac[1]).is_zero()) - + outc = pe.total_least_squares(oxc, oy, func) betac = outc.fit_parameters @@ -208,3 +208,41 @@ def test_odr_derivatives(): tfit = pe.fits.fit_general(x, y, func, base_step=0.1, step_ratio=1.1, num_steps=20) assert np.abs(np.max(np.array(list(fit1[1].deltas.values())) - np.array(list(tfit[1].deltas.values())))) < 10e-8 + + +def test_r_value_persistence(): + def f(a, x): + return a[0] + a[1] * x + + a = pe.pseudo_Obs(1.1, .1, 'a') + assert np.isclose(a.value, a.r_values['a']) + + a_2 = a ** 2 + assert np.isclose(a_2.value, a_2.r_values['a']) + + b = pe.pseudo_Obs(2.1, .2, 'b') + + y = [a, b] + [o.gamma_method() for o in y] + + fitp = pe.fits.least_squares([1, 2], y, f) + + assert np.isclose(fitp[0].value, fitp[0].r_values['a']) + assert np.isclose(fitp[0].value, fitp[0].r_values['b']) + assert np.isclose(fitp[1].value, fitp[1].r_values['a']) + assert np.isclose(fitp[1].value, fitp[1].r_values['b']) + + fitp = pe.fits.total_least_squares(y, y, f) + + assert np.isclose(fitp[0].value, fitp[0].r_values['a']) + assert np.isclose(fitp[0].value, fitp[0].r_values['b']) + assert np.isclose(fitp[1].value, fitp[1].r_values['a']) + assert np.isclose(fitp[1].value, fitp[1].r_values['b']) + + fitp = pe.fits.least_squares([1, 2], y, f, priors=y) + + assert np.isclose(fitp[0].value, fitp[0].r_values['a']) + assert np.isclose(fitp[0].value, fitp[0].r_values['b']) + assert np.isclose(fitp[1].value, fitp[1].r_values['a']) + assert np.isclose(fitp[1].value, fitp[1].r_values['b']) + From 2d72d730a56b5a392411d0c13f036f35ac166832 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 4 Dec 2021 13:09:30 +0000 Subject: [PATCH 5/9] test: r_value test for merge_obs added --- tests/obs_test.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/obs_test.py b/tests/obs_test.py index 2a66ce7c..ce9d7f41 100644 --- a/tests/obs_test.py +++ b/tests/obs_test.py @@ -381,6 +381,16 @@ def test_merge_obs(): assert diff == -(my_obs1.value + my_obs2.value) / 2 +def test_merge_obs_r_values(): + a1 = pe.pseudo_Obs(1.1, .1, 'a|1') + a2 = pe.pseudo_Obs(1.2, .1, 'a|2') + a = pe.merge_obs([a1, a2]) + + assert np.isclose(a.r_values['a|1'], a1.value) + assert np.isclose(a.r_values['a|2'], a2.value) + assert np.isclose(a.value, np.mean([a1.value, a2.value])) + + def test_correlate(): my_obs1 = pe.Obs([np.random.rand(100)], ['t']) my_obs2 = pe.Obs([np.random.rand(100)], ['t']) From 8163629efff1804dbd535eebd4ca21b0cb2250f4 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Sat, 4 Dec 2021 13:11:19 +0000 Subject: [PATCH 6/9] test: root test extended --- tests/roots_test.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/roots_test.py b/tests/roots_test.py index 8a4720d4..dbb27bbe 100644 --- a/tests/roots_test.py +++ b/tests/roots_test.py @@ -15,6 +15,7 @@ def test_root_linear(): my_root = pe.roots.find_root(my_obs, root_function) assert np.isclose(my_root.value, value) + assert np.isclose(my_root.value, my_root.r_values['t']) difference = my_obs - my_root assert difference.is_zero() From 12a93eafb02418c74907c57516c4f7a9c8e0c982 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 6 Dec 2021 10:54:06 +0000 Subject: [PATCH 7/9] feat: performance of export to jackknife improved --- pyerrors/obs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 5570f509..c1fcc8ca 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -623,9 +623,9 @@ class Obs: name = self.names[0] full_data = self.deltas[name] + self.r_values[name] n = full_data.size - mean = np.mean(full_data) + mean = self.value tmp_jacks = np.zeros(n + 1) - tmp_jacks[0] = self.value + tmp_jacks[0] = mean tmp_jacks[1:] = (n * mean - full_data) / (n - 1) return tmp_jacks From 7937635ca25e8d2a82ddc5e89bfbe0ae3ba0e429 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 6 Dec 2021 15:02:15 +0000 Subject: [PATCH 8/9] feat: check for name doublers in Obs.__init__ optimized --- pyerrors/obs.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index c1fcc8ca..38faa276 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -74,8 +74,10 @@ class Obs: if idl is not None: if len(idl) != len(names): raise Exception('Length of idl incompatible with samples and names.') - if len(names) != len(set(names)): - raise Exception('names are not unique.') + name_length = len(names) + if name_length > 1: + if name_length != len(set(names)): + raise Exception('names are not unique.') if not all(isinstance(x, str) for x in names): raise TypeError('All names have to be strings.') if min(len(x) for x in samples) <= 4: From 1f91175e50be85b3701b3909fcc64a853ff8c066 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 6 Dec 2021 15:20:36 +0000 Subject: [PATCH 9/9] feat: check for non string names in Obs.__init__ optimized --- pyerrors/obs.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 38faa276..c6c97197 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -78,8 +78,11 @@ class Obs: if name_length > 1: if name_length != len(set(names)): raise Exception('names are not unique.') - if not all(isinstance(x, str) for x in names): - raise TypeError('All names have to be strings.') + if not all(isinstance(x, str) for x in names): + raise TypeError('All names have to be strings.') + else: + if not isinstance(names[0], str): + raise TypeError('All names have to be strings.') if min(len(x) for x in samples) <= 4: raise Exception('Samples have to have at least 5 entries.')