Merge branch 'develop' into feature/covobs

This commit is contained in:
Simon Kuberski 2021-11-30 13:33:47 +01:00
commit 0e75cf48bc
2 changed files with 92 additions and 27 deletions

View file

@ -63,7 +63,7 @@ def read_meson_hd5(path, filestem, ens_id, meson='meson_0', tree='meson', idl=No
for outputs of the Meson module. Can be altered to read input
from other modules with similar structures.
idl : range
If specified only conifgurations in the given range are read in.
If specified only configurations in the given range are read in.
"""
files, idx = _get_files(path, filestem, idl)
@ -134,12 +134,15 @@ def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
"""Read hadrons ExternalLeg hdf5 file and output an array of CObs
Parameters
-----------------
path -- path to the files to read
filestem -- namestem of the files to read
ens_id -- name of the ensemble, required for internal bookkeeping
----------
path : str
path to the files to read
filestem : str
namestem of the files to read
ens_id : str
name of the ensemble, required for internal bookkeeping
idl : range
If specified only conifgurations in the given range are read in.
If specified only configurations in the given range are read in.
"""
files, idx = _get_files(path, filestem, idl)
@ -171,12 +174,15 @@ def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
"""Read hadrons Bilinear hdf5 file and output an array of CObs
Parameters
-----------------
path -- path to the files to read
filestem -- namestem of the files to read
ens_id -- name of the ensemble, required for internal bookkeeping
----------
path : str
path to the files to read
filestem : str
namestem of the files to read
ens_id : str
name of the ensemble, required for internal bookkeeping
idl : range
If specified only conifgurations in the given range are read in.
If specified only configurations in the given range are read in.
"""
files, idx = _get_files(path, filestem, idl)
@ -216,3 +222,62 @@ def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
result_dict[key] = Npr_matrix(matrix.swapaxes(1, 2).reshape((12, 12), order='F'), mom_in=mom_in, mom_out=mom_out)
return result_dict
def read_Fourquark_hd5(path, filestem, ens_id, idl=None):
"""Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
Parameters
----------
path : str
path to the files to read
filestem : str
namestem of the files to read
ens_id : str
name of the ensemble, required for internal bookkeeping
idl : range
If specified only configurations in the given range are read in.
"""
files, idx = _get_files(path, filestem, idl)
mom_in = None
mom_out = None
corr_data = {}
tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
for hd5_file in files:
file = h5py.File(path + '/' + hd5_file, "r")
for i in range(1):
name = file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8') + '_' + file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8')
if name not in corr_data:
corr_data[name] = []
raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
corr_data[name].append(raw_data)
if mom_in is None:
mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(' '), dtype=int)
if mom_out is None:
mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(' '), dtype=int)
file.close()
result_dict = {}
for key, data in corr_data.items():
local_data = np.array(data)
rolled_array = np.moveaxis(local_data, 0, 8)
matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
for index in np.ndindex(rolled_array.shape[:-1]):
real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
matrix[index] = CObs(real, imag)
result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
# result_dict[key] = Npr_matrix(matrix.swapaxes(1, 2).reshape((12, 12), order='F'), mom_in=mom_in, mom_out=mom_out)
return result_dict

View file

@ -186,49 +186,49 @@ def jack_matmul(*operands):
For large matrices this is considerably faster compared to matmul.
"""
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]
if any(isinstance(o.flat[0], CObs) for o in operands):
name = operands[0].flat[0].real.names[0]
idl = operands[0].flat[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()
for index, entry in np.ndenumerate(matrix):
base_matrix[index] = entry.real.export_jackknife() + 1j * entry.imag.export_jackknife()
return base_matrix
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, [idl]),
import_jackknife(entry.imag, name, [idl]))
for index, entry in np.ndenumerate(matrix):
base_matrix[index] = 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:]:
if isinstance(op[0, 0], CObs):
if isinstance(op.flat[0], CObs):
r = r @ _exp_to_jack(op)
else:
r = r @ op
return _imp_from_jack(r)
else:
name = operands[0][0, 0].names[0]
idl = operands[0][0, 0].idl[name]
name = operands[0].flat[0].names[0]
idl = operands[0].flat[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()
for index, entry in np.ndenumerate(matrix):
base_matrix[index] = entry.export_jackknife()
return base_matrix
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, [idl])
for index, entry in np.ndenumerate(matrix):
base_matrix[index] = import_jackknife(entry, name, [idl])
return base_matrix
r = _exp_to_jack(operands[0])
for op in operands[1:]:
if isinstance(op[0, 0], Obs):
if isinstance(op.flat[0], Obs):
r = r @ _exp_to_jack(op)
else:
r = r @ op