From cf1c38462d31511a9d71ca6a94ec902dd421ad34 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 18 Nov 2021 11:02:31 +0000 Subject: [PATCH 1/3] test: test for jack_matmul with multiple operands added --- tests/linalg_test.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/linalg_test.py b/tests/linalg_test.py index 4cfa9107..58c876ee 100644 --- a/tests/linalg_test.py +++ b/tests/linalg_test.py @@ -92,6 +92,17 @@ def test_multi_dot(): assert e.is_zero(), t +def test_jack_multi_dot(): + for dim in [2, 4, 8]: + my_array = get_real_matrix(dim) + + tt = pe.linalg.jack_matmul(my_array, my_array, my_array) - pe.linalg.matmul(my_array, my_array, my_array) + + for t, e in np.ndenumerate(tt): + assert e.is_zero(atol=1e-1), t + assert np.isclose(e.value, 0.0) + + def test_matmul_irregular_histories(): dim = 2 length = 500 From 0954ebee6eb58c5e79cd7e589f28a731115f592c Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 18 Nov 2021 11:03:55 +0000 Subject: [PATCH 2/3] test: test for jack_matmul with mutliple operands extended --- tests/linalg_test.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/linalg_test.py b/tests/linalg_test.py index 58c876ee..5af84b52 100644 --- a/tests/linalg_test.py +++ b/tests/linalg_test.py @@ -99,6 +99,8 @@ def test_jack_multi_dot(): tt = pe.linalg.jack_matmul(my_array, my_array, my_array) - pe.linalg.matmul(my_array, my_array, my_array) for t, e in np.ndenumerate(tt): + e.gamma_method() + assert e.is_zero_within_error(0.01) assert e.is_zero(atol=1e-1), t assert np.isclose(e.value, 0.0) From 28bf0f1701a1740a4e1a9ccc99c94d17c71f31f6 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 18 Nov 2021 11:17:20 +0000 Subject: [PATCH 3/3] feat: linalg.jack_matmul now also works with irregular monte carlo chains --- pyerrors/linalg.py | 20 +++++++++++++------- pyerrors/obs.py | 4 ++-- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/pyerrors/linalg.py b/pyerrors/linalg.py index e9f5dbd5..075b9605 100644 --- a/pyerrors/linalg.py +++ b/pyerrors/linalg.py @@ -187,40 +187,46 @@ def jack_matmul(*operands): """ if any(isinstance(o[0, 0], CObs) for o in operands): + name = operands[0][0, 0].real.names[0] + idl = operands[0][0, 0].real.idl[name] + def _exp_to_jack(matrix): base_matrix = np.empty_like(matrix) for (n, m), entry in np.ndenumerate(matrix): base_matrix[n, m] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife() return base_matrix - def _imp_from_jack(matrix, name): + def _imp_from_jack(matrix): base_matrix = np.empty_like(matrix) for (n, m), entry in np.ndenumerate(matrix): - base_matrix[n, m] = CObs(import_jackknife(entry.real, name), - import_jackknife(entry.imag, name)) + base_matrix[n, m] = CObs(import_jackknife(entry.real, name, [idl]), + import_jackknife(entry.imag, name, [idl])) return base_matrix r = _exp_to_jack(operands[0]) for op in operands[1:]: r = r @ _exp_to_jack(op) - return _imp_from_jack(r, op.ravel()[0].real.names[0]) + return _imp_from_jack(r) else: + name = operands[0][0, 0].names[0] + idl = operands[0][0, 0].idl[name] + def _exp_to_jack(matrix): base_matrix = np.empty_like(matrix) for (n, m), entry in np.ndenumerate(matrix): base_matrix[n, m] = entry.export_jackknife() return base_matrix - def _imp_from_jack(matrix, name): + def _imp_from_jack(matrix): base_matrix = np.empty_like(matrix) for (n, m), entry in np.ndenumerate(matrix): - base_matrix[n, m] = import_jackknife(entry, name) + base_matrix[n, m] = import_jackknife(entry, name, [idl]) return base_matrix r = _exp_to_jack(operands[0]) for op in operands[1:]: r = r @ _exp_to_jack(op) - return _imp_from_jack(r, op.ravel()[0].names[0]) + return _imp_from_jack(r) def inv(x): diff --git a/pyerrors/obs.py b/pyerrors/obs.py index 1e594ea4..56a67e57 100644 --- a/pyerrors/obs.py +++ b/pyerrors/obs.py @@ -1559,7 +1559,7 @@ def load_object(path): return pickle.load(file) -def import_jackknife(jacks, name): +def import_jackknife(jacks, name, idl=None): """Imports jackknife samples and returns an Obs Parameters @@ -1573,7 +1573,7 @@ def import_jackknife(jacks, name): length = len(jacks) - 1 prj = (np.ones((length, length)) - (length - 1) * np.identity(length)) samples = jacks[1:] @ prj - new_obs = Obs([samples], [name]) + new_obs = Obs([samples], [name], idl=idl) new_obs._value = jacks[0] return new_obs