feat: linalg.jack_matmul now also works with irregular monte carlo

chains
This commit is contained in:
Fabian Joswig 2021-11-18 11:17:20 +00:00
parent 0954ebee6e
commit 28bf0f1701
2 changed files with 15 additions and 9 deletions

View file

@ -187,40 +187,46 @@ def jack_matmul(*operands):
""" """
if any(isinstance(o[0, 0], CObs) for o in 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): def _exp_to_jack(matrix):
base_matrix = np.empty_like(matrix) base_matrix = np.empty_like(matrix)
for (n, m), entry in np.ndenumerate(matrix): for (n, m), entry in np.ndenumerate(matrix):
base_matrix[n, m] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife() base_matrix[n, m] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife()
return base_matrix return base_matrix
def _imp_from_jack(matrix, name): def _imp_from_jack(matrix):
base_matrix = np.empty_like(matrix) base_matrix = np.empty_like(matrix)
for (n, m), entry in np.ndenumerate(matrix): for (n, m), entry in np.ndenumerate(matrix):
base_matrix[n, m] = CObs(import_jackknife(entry.real, name), base_matrix[n, m] = CObs(import_jackknife(entry.real, name, [idl]),
import_jackknife(entry.imag, name)) import_jackknife(entry.imag, name, [idl]))
return base_matrix return base_matrix
r = _exp_to_jack(operands[0]) r = _exp_to_jack(operands[0])
for op in operands[1:]: for op in operands[1:]:
r = r @ _exp_to_jack(op) r = r @ _exp_to_jack(op)
return _imp_from_jack(r, op.ravel()[0].real.names[0]) return _imp_from_jack(r)
else: else:
name = operands[0][0, 0].names[0]
idl = operands[0][0, 0].idl[name]
def _exp_to_jack(matrix): def _exp_to_jack(matrix):
base_matrix = np.empty_like(matrix) base_matrix = np.empty_like(matrix)
for (n, m), entry in np.ndenumerate(matrix): for (n, m), entry in np.ndenumerate(matrix):
base_matrix[n, m] = entry.export_jackknife() base_matrix[n, m] = entry.export_jackknife()
return base_matrix return base_matrix
def _imp_from_jack(matrix, name): def _imp_from_jack(matrix):
base_matrix = np.empty_like(matrix) base_matrix = np.empty_like(matrix)
for (n, m), entry in np.ndenumerate(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 return base_matrix
r = _exp_to_jack(operands[0]) r = _exp_to_jack(operands[0])
for op in operands[1:]: for op in operands[1:]:
r = r @ _exp_to_jack(op) r = r @ _exp_to_jack(op)
return _imp_from_jack(r, op.ravel()[0].names[0]) return _imp_from_jack(r)
def inv(x): def inv(x):

View file

@ -1559,7 +1559,7 @@ def load_object(path):
return pickle.load(file) return pickle.load(file)
def import_jackknife(jacks, name): def import_jackknife(jacks, name, idl=None):
"""Imports jackknife samples and returns an Obs """Imports jackknife samples and returns an Obs
Parameters Parameters
@ -1573,7 +1573,7 @@ def import_jackknife(jacks, name):
length = len(jacks) - 1 length = len(jacks) - 1
prj = (np.ones((length, length)) - (length - 1) * np.identity(length)) prj = (np.ones((length, length)) - (length - 1) * np.identity(length))
samples = jacks[1:] @ prj samples = jacks[1:] @ prj
new_obs = Obs([samples], [name]) new_obs = Obs([samples], [name], idl=idl)
new_obs._value = jacks[0] new_obs._value = jacks[0]
return new_obs return new_obs