From 9dd672688f4c548f5f672c6c06504f593a570bf3 Mon Sep 17 00:00:00 2001 From: fjosw Date: Mon, 28 Nov 2022 15:32:37 +0000 Subject: [PATCH] Documentation updated --- docs/pyerrors/input/hadrons.html | 1046 +++++++++++++++--------------- 1 file changed, 525 insertions(+), 521 deletions(-) diff --git a/docs/pyerrors/input/hadrons.html b/docs/pyerrors/input/hadrons.html index 61e9aef4..e6ff3b26 100644 --- a/docs/pyerrors/input/hadrons.html +++ b/docs/pyerrors/input/hadrons.html @@ -244,319 +244,321 @@ 142 for explore_file in explore_path.iterdir(): 143 if explore_file.is_file(): 144 stem = explore_file.with_suffix("").with_suffix("").as_posix().split("/")[-1] -145 -146 file_list = [] -147 for dir in directories: -148 tmp_path = Path(path + "/" + dir) -149 file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5") -150 -151 corr_data = {} +145 else: +146 continue +147 +148 file_list = [] +149 for dir in directories: +150 tmp_path = Path(path + "/" + dir) +151 file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5") 152 -153 for diagram in diagrams: -154 corr_data[diagram] = [] -155 -156 for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))): -157 h5file = h5py.File(hd5_file) -158 -159 if n_file == 0: -160 if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...": -161 raise Exception("Routine is only implemented for files containing inversions on all timeslices.") -162 -163 Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0] +153 corr_data = {} +154 +155 for diagram in diagrams: +156 corr_data[diagram] = [] +157 +158 for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))): +159 h5file = h5py.File(hd5_file) +160 +161 if n_file == 0: +162 if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...": +163 raise Exception("Routine is only implemented for files containing inversions on all timeslices.") 164 -165 identifier = [] -166 for in_file in range(len(h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.keys()) - 1): -167 encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file)) -168 full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_") -169 my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3]) -170 identifier.append(my_tuple) -171 identifier = tuple(identifier) -172 # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case. -173 -174 for diagram in diagrams: -175 real_data = np.zeros(Nt) -176 for x0 in range(Nt): -177 raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:]["re"].astype(np.double) -178 real_data += np.roll(raw_data, -x0) -179 real_data /= Nt -180 -181 corr_data[diagram].append(real_data) -182 h5file.close() -183 -184 res_dict[str(identifier)] = {} +165 Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0] +166 +167 identifier = [] +168 for in_file in range(len(h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.keys()) - 1): +169 encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file)) +170 full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_") +171 my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3]) +172 identifier.append(my_tuple) +173 identifier = tuple(identifier) +174 # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case. +175 +176 for diagram in diagrams: +177 real_data = np.zeros(Nt) +178 for x0 in range(Nt): +179 raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:]["re"].astype(np.double) +180 real_data += np.roll(raw_data, -x0) +181 real_data /= Nt +182 +183 corr_data[diagram].append(real_data) +184 h5file.close() 185 -186 for diagram in diagrams: +186 res_dict[str(identifier)] = {} 187 -188 tmp_data = np.array(corr_data[diagram]) +188 for diagram in diagrams: 189 -190 l_obs = [] -191 for c in tmp_data.T: -192 l_obs.append(Obs([c], [ens_id], idl=[idx])) -193 -194 corr = Corr(l_obs) -195 corr.tag = str(identifier) -196 -197 res_dict[str(identifier)][diagram] = corr +190 tmp_data = np.array(corr_data[diagram]) +191 +192 l_obs = [] +193 for c in tmp_data.T: +194 l_obs.append(Obs([c], [ens_id], idl=[idx])) +195 +196 corr = Corr(l_obs) +197 corr.tag = str(identifier) 198 -199 return res_dict +199 res_dict[str(identifier)][diagram] = corr 200 -201 -202class Npr_matrix(np.ndarray): +201 return res_dict +202 203 -204 def __new__(cls, input_array, mom_in=None, mom_out=None): -205 obj = np.asarray(input_array).view(cls) -206 obj.mom_in = mom_in -207 obj.mom_out = mom_out -208 return obj -209 -210 @property -211 def g5H(self): -212 """Gamma_5 hermitean conjugate -213 -214 Uses the fact that the propagator is gamma5 hermitean, so just the -215 in and out momenta of the propagator are exchanged. -216 """ -217 return Npr_matrix(self, -218 mom_in=self.mom_out, -219 mom_out=self.mom_in) -220 -221 def _propagate_mom(self, other, name): -222 s_mom = getattr(self, name, None) -223 o_mom = getattr(other, name, None) -224 if s_mom is not None and o_mom is not None: -225 if not np.allclose(s_mom, o_mom): -226 raise Exception(name + ' does not match.') -227 return o_mom if o_mom is not None else s_mom -228 -229 def __matmul__(self, other): -230 return self.__new__(Npr_matrix, -231 super().__matmul__(other), -232 self._propagate_mom(other, 'mom_in'), -233 self._propagate_mom(other, 'mom_out')) -234 -235 def __array_finalize__(self, obj): -236 if obj is None: -237 return -238 self.mom_in = getattr(obj, 'mom_in', None) -239 self.mom_out = getattr(obj, 'mom_out', None) -240 -241 -242def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None): -243 """Read hadrons ExternalLeg hdf5 file and output an array of CObs -244 -245 Parameters -246 ---------- -247 path : str -248 path to the files to read -249 filestem : str -250 namestem of the files to read -251 ens_id : str -252 name of the ensemble, required for internal bookkeeping -253 idl : range -254 If specified only configurations in the given range are read in. -255 """ -256 -257 files, idx = _get_files(path, filestem, idl) +204class Npr_matrix(np.ndarray): +205 +206 def __new__(cls, input_array, mom_in=None, mom_out=None): +207 obj = np.asarray(input_array).view(cls) +208 obj.mom_in = mom_in +209 obj.mom_out = mom_out +210 return obj +211 +212 @property +213 def g5H(self): +214 """Gamma_5 hermitean conjugate +215 +216 Uses the fact that the propagator is gamma5 hermitean, so just the +217 in and out momenta of the propagator are exchanged. +218 """ +219 return Npr_matrix(self, +220 mom_in=self.mom_out, +221 mom_out=self.mom_in) +222 +223 def _propagate_mom(self, other, name): +224 s_mom = getattr(self, name, None) +225 o_mom = getattr(other, name, None) +226 if s_mom is not None and o_mom is not None: +227 if not np.allclose(s_mom, o_mom): +228 raise Exception(name + ' does not match.') +229 return o_mom if o_mom is not None else s_mom +230 +231 def __matmul__(self, other): +232 return self.__new__(Npr_matrix, +233 super().__matmul__(other), +234 self._propagate_mom(other, 'mom_in'), +235 self._propagate_mom(other, 'mom_out')) +236 +237 def __array_finalize__(self, obj): +238 if obj is None: +239 return +240 self.mom_in = getattr(obj, 'mom_in', None) +241 self.mom_out = getattr(obj, 'mom_out', None) +242 +243 +244def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None): +245 """Read hadrons ExternalLeg hdf5 file and output an array of CObs +246 +247 Parameters +248 ---------- +249 path : str +250 path to the files to read +251 filestem : str +252 namestem of the files to read +253 ens_id : str +254 name of the ensemble, required for internal bookkeeping +255 idl : range +256 If specified only configurations in the given range are read in. +257 """ 258 -259 mom = None +259 files, idx = _get_files(path, filestem, idl) 260 -261 corr_data = [] -262 for hd5_file in files: -263 file = h5py.File(path + '/' + hd5_file, "r") -264 raw_data = file['ExternalLeg/corr'][0][0].view('complex') -265 corr_data.append(raw_data) -266 if mom is None: -267 mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) -268 file.close() -269 corr_data = np.array(corr_data) -270 -271 rolled_array = np.rollaxis(corr_data, 0, 5) +261 mom = None +262 +263 corr_data = [] +264 for hd5_file in files: +265 file = h5py.File(path + '/' + hd5_file, "r") +266 raw_data = file['ExternalLeg/corr'][0][0].view('complex') +267 corr_data.append(raw_data) +268 if mom is None: +269 mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) +270 file.close() +271 corr_data = np.array(corr_data) 272 -273 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) -274 for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]): -275 real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx]) -276 imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx]) -277 matrix[si, sj, ci, cj] = CObs(real, imag) -278 -279 return Npr_matrix(matrix, mom_in=mom) +273 rolled_array = np.rollaxis(corr_data, 0, 5) +274 +275 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) +276 for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]): +277 real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx]) +278 imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx]) +279 matrix[si, sj, ci, cj] = CObs(real, imag) 280 -281 -282def read_Bilinear_hd5(path, filestem, ens_id, idl=None): -283 """Read hadrons Bilinear hdf5 file and output an array of CObs -284 -285 Parameters -286 ---------- -287 path : str -288 path to the files to read -289 filestem : str -290 namestem of the files to read -291 ens_id : str -292 name of the ensemble, required for internal bookkeeping -293 idl : range -294 If specified only configurations in the given range are read in. -295 """ -296 -297 files, idx = _get_files(path, filestem, idl) +281 return Npr_matrix(matrix, mom_in=mom) +282 +283 +284def read_Bilinear_hd5(path, filestem, ens_id, idl=None): +285 """Read hadrons Bilinear hdf5 file and output an array of CObs +286 +287 Parameters +288 ---------- +289 path : str +290 path to the files to read +291 filestem : str +292 namestem of the files to read +293 ens_id : str +294 name of the ensemble, required for internal bookkeeping +295 idl : range +296 If specified only configurations in the given range are read in. +297 """ 298 -299 mom_in = None -300 mom_out = None -301 -302 corr_data = {} -303 for hd5_file in files: -304 file = h5py.File(path + '/' + hd5_file, "r") -305 for i in range(16): -306 name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8') -307 if name not in corr_data: -308 corr_data[name] = [] -309 raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex') -310 corr_data[name].append(raw_data) -311 if mom_in is None: -312 mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) -313 if mom_out is None: -314 mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float) -315 -316 file.close() +299 files, idx = _get_files(path, filestem, idl) +300 +301 mom_in = None +302 mom_out = None +303 +304 corr_data = {} +305 for hd5_file in files: +306 file = h5py.File(path + '/' + hd5_file, "r") +307 for i in range(16): +308 name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8') +309 if name not in corr_data: +310 corr_data[name] = [] +311 raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex') +312 corr_data[name].append(raw_data) +313 if mom_in is None: +314 mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) +315 if mom_out is None: +316 mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float) 317 -318 result_dict = {} +318 file.close() 319 -320 for key, data in corr_data.items(): -321 local_data = np.array(data) -322 -323 rolled_array = np.rollaxis(local_data, 0, 5) +320 result_dict = {} +321 +322 for key, data in corr_data.items(): +323 local_data = np.array(data) 324 -325 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) -326 for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]): -327 real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx]) -328 imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx]) -329 matrix[si, sj, ci, cj] = CObs(real, imag) -330 -331 result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out) +325 rolled_array = np.rollaxis(local_data, 0, 5) +326 +327 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) +328 for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]): +329 real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx]) +330 imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx]) +331 matrix[si, sj, ci, cj] = CObs(real, imag) 332 -333 return result_dict +333 result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out) 334 -335 -336def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]): -337 """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs -338 -339 Parameters -340 ---------- -341 path : str -342 path to the files to read -343 filestem : str -344 namestem of the files to read -345 ens_id : str -346 name of the ensemble, required for internal bookkeeping -347 idl : range -348 If specified only configurations in the given range are read in. -349 vertices : list -350 Vertex functions to be extracted. -351 """ -352 -353 files, idx = _get_files(path, filestem, idl) +335 return result_dict +336 +337 +338def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]): +339 """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs +340 +341 Parameters +342 ---------- +343 path : str +344 path to the files to read +345 filestem : str +346 namestem of the files to read +347 ens_id : str +348 name of the ensemble, required for internal bookkeeping +349 idl : range +350 If specified only configurations in the given range are read in. +351 vertices : list +352 Vertex functions to be extracted. +353 """ 354 -355 mom_in = None -356 mom_out = None -357 -358 vertex_names = [] -359 for vertex in vertices: -360 vertex_names += _get_lorentz_names(vertex) -361 -362 corr_data = {} +355 files, idx = _get_files(path, filestem, idl) +356 +357 mom_in = None +358 mom_out = None +359 +360 vertex_names = [] +361 for vertex in vertices: +362 vertex_names += _get_lorentz_names(vertex) 363 -364 tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_' +364 corr_data = {} 365 -366 for hd5_file in files: -367 file = h5py.File(path + '/' + hd5_file, "r") -368 -369 for i in range(32): -370 name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8')) -371 if name in vertex_names: -372 if name not in corr_data: -373 corr_data[name] = [] -374 raw_data = file[tree + str(i) + '/corr'][0][0].view('complex') -375 corr_data[name].append(raw_data) -376 if mom_in is None: -377 mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) -378 if mom_out is None: -379 mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float) -380 -381 file.close() +366 tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_' +367 +368 for hd5_file in files: +369 file = h5py.File(path + '/' + hd5_file, "r") +370 +371 for i in range(32): +372 name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8')) +373 if name in vertex_names: +374 if name not in corr_data: +375 corr_data[name] = [] +376 raw_data = file[tree + str(i) + '/corr'][0][0].view('complex') +377 corr_data[name].append(raw_data) +378 if mom_in is None: +379 mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) +380 if mom_out is None: +381 mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float) 382 -383 intermediate_dict = {} +383 file.close() 384 -385 for vertex in vertices: -386 lorentz_names = _get_lorentz_names(vertex) -387 for v_name in lorentz_names: -388 if v_name in [('SigmaXY', 'SigmaZT'), -389 ('SigmaXT', 'SigmaYZ'), -390 ('SigmaYZ', 'SigmaXT'), -391 ('SigmaZT', 'SigmaXY')]: -392 sign = -1 -393 else: -394 sign = 1 -395 if vertex not in intermediate_dict: -396 intermediate_dict[vertex] = sign * np.array(corr_data[v_name]) -397 else: -398 intermediate_dict[vertex] += sign * np.array(corr_data[v_name]) -399 -400 result_dict = {} +385 intermediate_dict = {} +386 +387 for vertex in vertices: +388 lorentz_names = _get_lorentz_names(vertex) +389 for v_name in lorentz_names: +390 if v_name in [('SigmaXY', 'SigmaZT'), +391 ('SigmaXT', 'SigmaYZ'), +392 ('SigmaYZ', 'SigmaXT'), +393 ('SigmaZT', 'SigmaXY')]: +394 sign = -1 +395 else: +396 sign = 1 +397 if vertex not in intermediate_dict: +398 intermediate_dict[vertex] = sign * np.array(corr_data[v_name]) +399 else: +400 intermediate_dict[vertex] += sign * np.array(corr_data[v_name]) 401 -402 for key, data in intermediate_dict.items(): +402 result_dict = {} 403 -404 rolled_array = np.moveaxis(data, 0, 8) +404 for key, data in intermediate_dict.items(): 405 -406 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) -407 for index in np.ndindex(rolled_array.shape[:-1]): -408 real = Obs([rolled_array[index].real], [ens_id], idl=[idx]) -409 imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx]) -410 matrix[index] = CObs(real, imag) -411 -412 result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out) +406 rolled_array = np.moveaxis(data, 0, 8) +407 +408 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) +409 for index in np.ndindex(rolled_array.shape[:-1]): +410 real = Obs([rolled_array[index].real], [ens_id], idl=[idx]) +411 imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx]) +412 matrix[index] = CObs(real, imag) 413 -414 return result_dict +414 result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out) 415 -416 -417def _get_lorentz_names(name): -418 lorentz_index = ['X', 'Y', 'Z', 'T'] -419 -420 res = [] +416 return result_dict +417 +418 +419def _get_lorentz_names(name): +420 lorentz_index = ['X', 'Y', 'Z', 'T'] 421 -422 if name == "TT": -423 for i in range(4): -424 for j in range(i + 1, 4): -425 res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[i] + lorentz_index[j])) -426 return res -427 -428 if name == "TTtilde": -429 for i in range(4): -430 for j in range(i + 1, 4): -431 for k in range(4): -432 for o in range(k + 1, 4): -433 fac = epsilon_tensor_rank4(i, j, k, o) -434 if not np.isclose(fac, 0.0): -435 res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[k] + lorentz_index[o])) -436 return res -437 -438 assert len(name) == 2 +422 res = [] +423 +424 if name == "TT": +425 for i in range(4): +426 for j in range(i + 1, 4): +427 res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[i] + lorentz_index[j])) +428 return res +429 +430 if name == "TTtilde": +431 for i in range(4): +432 for j in range(i + 1, 4): +433 for k in range(4): +434 for o in range(k + 1, 4): +435 fac = epsilon_tensor_rank4(i, j, k, o) +436 if not np.isclose(fac, 0.0): +437 res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[k] + lorentz_index[o])) +438 return res 439 -440 if 'S' in name or 'P' in name: -441 if not set(name) <= set(['S', 'P']): -442 raise Exception("'" + name + "' is not a Lorentz scalar") -443 -444 g_names = {'S': 'Identity', -445 'P': 'Gamma5'} -446 -447 res.append((g_names[name[0]], g_names[name[1]])) +440 assert len(name) == 2 +441 +442 if 'S' in name or 'P' in name: +443 if not set(name) <= set(['S', 'P']): +444 raise Exception("'" + name + "' is not a Lorentz scalar") +445 +446 g_names = {'S': 'Identity', +447 'P': 'Gamma5'} 448 -449 else: -450 if not set(name) <= set(['V', 'A']): -451 raise Exception("'" + name + "' is not a Lorentz scalar") -452 -453 for ind in lorentz_index: -454 res.append(('Gamma' + ind + (name[0] == 'A') * 'Gamma5', -455 'Gamma' + ind + (name[1] == 'A') * 'Gamma5')) -456 -457 return res +449 res.append((g_names[name[0]], g_names[name[1]])) +450 +451 else: +452 if not set(name) <= set(['V', 'A']): +453 raise Exception("'" + name + "' is not a Lorentz scalar") +454 +455 for ind in lorentz_index: +456 res.append(('Gamma' + ind + (name[0] == 'A') * 'Gamma5', +457 'Gamma' + ind + (name[1] == 'A') * 'Gamma5')) +458 +459 return res @@ -697,61 +699,63 @@ If specified only configurations in the given range are read in. 143 for explore_file in explore_path.iterdir(): 144 if explore_file.is_file(): 145 stem = explore_file.with_suffix("").with_suffix("").as_posix().split("/")[-1] -146 -147 file_list = [] -148 for dir in directories: -149 tmp_path = Path(path + "/" + dir) -150 file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5") -151 -152 corr_data = {} +146 else: +147 continue +148 +149 file_list = [] +150 for dir in directories: +151 tmp_path = Path(path + "/" + dir) +152 file_list.append((tmp_path / stem).as_posix() + tmp_path.suffix + ".h5") 153 -154 for diagram in diagrams: -155 corr_data[diagram] = [] -156 -157 for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))): -158 h5file = h5py.File(hd5_file) -159 -160 if n_file == 0: -161 if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...": -162 raise Exception("Routine is only implemented for files containing inversions on all timeslices.") -163 -164 Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0] +154 corr_data = {} +155 +156 for diagram in diagrams: +157 corr_data[diagram] = [] +158 +159 for n_file, (hd5_file, n_traj) in enumerate(zip(file_list, list(idx))): +160 h5file = h5py.File(hd5_file) +161 +162 if n_file == 0: +163 if h5file["DistillationContraction/Metadata"].attrs.get("TimeSources")[0].decode() != "0...": +164 raise Exception("Routine is only implemented for files containing inversions on all timeslices.") 165 -166 identifier = [] -167 for in_file in range(len(h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.keys()) - 1): -168 encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file)) -169 full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_") -170 my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3]) -171 identifier.append(my_tuple) -172 identifier = tuple(identifier) -173 # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case. -174 -175 for diagram in diagrams: -176 real_data = np.zeros(Nt) -177 for x0 in range(Nt): -178 raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:]["re"].astype(np.double) -179 real_data += np.roll(raw_data, -x0) -180 real_data /= Nt -181 -182 corr_data[diagram].append(real_data) -183 h5file.close() -184 -185 res_dict[str(identifier)] = {} +166 Nt = h5file["DistillationContraction/Metadata"].attrs.get("Nt")[0] +167 +168 identifier = [] +169 for in_file in range(len(h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.keys()) - 1): +170 encoded_info = h5file["DistillationContraction/Metadata/DmfInputFiles"].attrs.get("DmfInputFiles_" + str(in_file)) +171 full_info = encoded_info[0].decode().split("/")[-1].replace(".h5", "").split("_") +172 my_tuple = (full_info[0], full_info[1][1:], full_info[2], full_info[3]) +173 identifier.append(my_tuple) +174 identifier = tuple(identifier) +175 # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case. +176 +177 for diagram in diagrams: +178 real_data = np.zeros(Nt) +179 for x0 in range(Nt): +180 raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:]["re"].astype(np.double) +181 real_data += np.roll(raw_data, -x0) +182 real_data /= Nt +183 +184 corr_data[diagram].append(real_data) +185 h5file.close() 186 -187 for diagram in diagrams: +187 res_dict[str(identifier)] = {} 188 -189 tmp_data = np.array(corr_data[diagram]) +189 for diagram in diagrams: 190 -191 l_obs = [] -192 for c in tmp_data.T: -193 l_obs.append(Obs([c], [ens_id], idl=[idx])) -194 -195 corr = Corr(l_obs) -196 corr.tag = str(identifier) -197 -198 res_dict[str(identifier)][diagram] = corr +191 tmp_data = np.array(corr_data[diagram]) +192 +193 l_obs = [] +194 for c in tmp_data.T: +195 l_obs.append(Obs([c], [ens_id], idl=[idx])) +196 +197 corr = Corr(l_obs) +198 corr.tag = str(identifier) 199 -200 return res_dict +200 res_dict[str(identifier)][diagram] = corr +201 +202 return res_dict @@ -784,44 +788,44 @@ If specified only configurations in the given range are read in. - 
203class Npr_matrix(np.ndarray):
-204
-205    def __new__(cls, input_array, mom_in=None, mom_out=None):
-206        obj = np.asarray(input_array).view(cls)
-207        obj.mom_in = mom_in
-208        obj.mom_out = mom_out
-209        return obj
-210
-211    @property
-212    def g5H(self):
-213        """Gamma_5 hermitean conjugate
-214
-215        Uses the fact that the propagator is gamma5 hermitean, so just the
-216        in and out momenta of the propagator are exchanged.
-217        """
-218        return Npr_matrix(self,
-219                          mom_in=self.mom_out,
-220                          mom_out=self.mom_in)
-221
-222    def _propagate_mom(self, other, name):
-223        s_mom = getattr(self, name, None)
-224        o_mom = getattr(other, name, None)
-225        if s_mom is not None and o_mom is not None:
-226            if not np.allclose(s_mom, o_mom):
-227                raise Exception(name + ' does not match.')
-228        return o_mom if o_mom is not None else s_mom
-229
-230    def __matmul__(self, other):
-231        return self.__new__(Npr_matrix,
-232                            super().__matmul__(other),
-233                            self._propagate_mom(other, 'mom_in'),
-234                            self._propagate_mom(other, 'mom_out'))
-235
-236    def __array_finalize__(self, obj):
-237        if obj is None:
-238            return
-239        self.mom_in = getattr(obj, 'mom_in', None)
-240        self.mom_out = getattr(obj, 'mom_out', None)
+            
205class Npr_matrix(np.ndarray):
+206
+207    def __new__(cls, input_array, mom_in=None, mom_out=None):
+208        obj = np.asarray(input_array).view(cls)
+209        obj.mom_in = mom_in
+210        obj.mom_out = mom_out
+211        return obj
+212
+213    @property
+214    def g5H(self):
+215        """Gamma_5 hermitean conjugate
+216
+217        Uses the fact that the propagator is gamma5 hermitean, so just the
+218        in and out momenta of the propagator are exchanged.
+219        """
+220        return Npr_matrix(self,
+221                          mom_in=self.mom_out,
+222                          mom_out=self.mom_in)
+223
+224    def _propagate_mom(self, other, name):
+225        s_mom = getattr(self, name, None)
+226        o_mom = getattr(other, name, None)
+227        if s_mom is not None and o_mom is not None:
+228            if not np.allclose(s_mom, o_mom):
+229                raise Exception(name + ' does not match.')
+230        return o_mom if o_mom is not None else s_mom
+231
+232    def __matmul__(self, other):
+233        return self.__new__(Npr_matrix,
+234                            super().__matmul__(other),
+235                            self._propagate_mom(other, 'mom_in'),
+236                            self._propagate_mom(other, 'mom_out'))
+237
+238    def __array_finalize__(self, obj):
+239        if obj is None:
+240            return
+241        self.mom_in = getattr(obj, 'mom_in', None)
+242        self.mom_out = getattr(obj, 'mom_out', None)
 

 
 
@@ -1076,44 +1080,44 @@ in and out momenta of the propagator are exchanged.
- 
243def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
-244    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
-245
-246    Parameters
-247    ----------
-248    path : str
-249        path to the files to read
-250    filestem : str
-251        namestem of the files to read
-252    ens_id : str
-253        name of the ensemble, required for internal bookkeeping
-254    idl : range
-255        If specified only configurations in the given range are read in.
-256    """
-257
-258    files, idx = _get_files(path, filestem, idl)
+            
245def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
+246    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
+247
+248    Parameters
+249    ----------
+250    path : str
+251        path to the files to read
+252    filestem : str
+253        namestem of the files to read
+254    ens_id : str
+255        name of the ensemble, required for internal bookkeeping
+256    idl : range
+257        If specified only configurations in the given range are read in.
+258    """
 259
-260    mom = None
+260    files, idx = _get_files(path, filestem, idl)
 261
-262    corr_data = []
-263    for hd5_file in files:
-264        file = h5py.File(path + '/' + hd5_file, "r")
-265        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
-266        corr_data.append(raw_data)
-267        if mom is None:
-268            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-269        file.close()
-270    corr_data = np.array(corr_data)
-271
-272    rolled_array = np.rollaxis(corr_data, 0, 5)
+262    mom = None
+263
+264    corr_data = []
+265    for hd5_file in files:
+266        file = h5py.File(path + '/' + hd5_file, "r")
+267        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
+268        corr_data.append(raw_data)
+269        if mom is None:
+270            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+271        file.close()
+272    corr_data = np.array(corr_data)
 273
-274    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-275    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
-276        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
-277        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
-278        matrix[si, sj, ci, cj] = CObs(real, imag)
-279
-280    return Npr_matrix(matrix, mom_in=mom)
+274    rolled_array = np.rollaxis(corr_data, 0, 5)
+275
+276    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+277    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
+278        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
+279        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
+280        matrix[si, sj, ci, cj] = CObs(real, imag)
+281
+282    return Npr_matrix(matrix, mom_in=mom)
 

 
 
@@ -1146,58 +1150,58 @@ If specified only configurations in the given range are read in.
- 
283def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
-284    """Read hadrons Bilinear hdf5 file and output an array of CObs
-285
-286    Parameters
-287    ----------
-288    path : str
-289        path to the files to read
-290    filestem : str
-291        namestem of the files to read
-292    ens_id : str
-293        name of the ensemble, required for internal bookkeeping
-294    idl : range
-295        If specified only configurations in the given range are read in.
-296    """
-297
-298    files, idx = _get_files(path, filestem, idl)
+            
285def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
+286    """Read hadrons Bilinear hdf5 file and output an array of CObs
+287
+288    Parameters
+289    ----------
+290    path : str
+291        path to the files to read
+292    filestem : str
+293        namestem of the files to read
+294    ens_id : str
+295        name of the ensemble, required for internal bookkeeping
+296    idl : range
+297        If specified only configurations in the given range are read in.
+298    """
 299
-300    mom_in = None
-301    mom_out = None
-302
-303    corr_data = {}
-304    for hd5_file in files:
-305        file = h5py.File(path + '/' + hd5_file, "r")
-306        for i in range(16):
-307            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
-308            if name not in corr_data:
-309                corr_data[name] = []
-310            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
-311            corr_data[name].append(raw_data)
-312            if mom_in is None:
-313                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-314            if mom_out is None:
-315                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
-316
-317        file.close()
+300    files, idx = _get_files(path, filestem, idl)
+301
+302    mom_in = None
+303    mom_out = None
+304
+305    corr_data = {}
+306    for hd5_file in files:
+307        file = h5py.File(path + '/' + hd5_file, "r")
+308        for i in range(16):
+309            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
+310            if name not in corr_data:
+311                corr_data[name] = []
+312            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
+313            corr_data[name].append(raw_data)
+314            if mom_in is None:
+315                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+316            if mom_out is None:
+317                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
 318
-319    result_dict = {}
+319        file.close()
 320
-321    for key, data in corr_data.items():
-322        local_data = np.array(data)
-323
-324        rolled_array = np.rollaxis(local_data, 0, 5)
+321    result_dict = {}
+322
+323    for key, data in corr_data.items():
+324        local_data = np.array(data)
 325
-326        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-327        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
-328            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
-329            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
-330            matrix[si, sj, ci, cj] = CObs(real, imag)
-331
-332        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+326        rolled_array = np.rollaxis(local_data, 0, 5)
+327
+328        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+329        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
+330            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
+331            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
+332            matrix[si, sj, ci, cj] = CObs(real, imag)
 333
-334    return result_dict
+334        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+335
+336    return result_dict
 

 
 
@@ -1230,85 +1234,85 @@ If specified only configurations in the given range are read in.
- 
337def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
-338    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
-339
-340    Parameters
-341    ----------
-342    path : str
-343        path to the files to read
-344    filestem : str
-345        namestem of the files to read
-346    ens_id : str
-347        name of the ensemble, required for internal bookkeeping
-348    idl : range
-349        If specified only configurations in the given range are read in.
-350    vertices : list
-351        Vertex functions to be extracted.
-352    """
-353
-354    files, idx = _get_files(path, filestem, idl)
+            
339def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
+340    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
+341
+342    Parameters
+343    ----------
+344    path : str
+345        path to the files to read
+346    filestem : str
+347        namestem of the files to read
+348    ens_id : str
+349        name of the ensemble, required for internal bookkeeping
+350    idl : range
+351        If specified only configurations in the given range are read in.
+352    vertices : list
+353        Vertex functions to be extracted.
+354    """
 355
-356    mom_in = None
-357    mom_out = None
-358
-359    vertex_names = []
-360    for vertex in vertices:
-361        vertex_names += _get_lorentz_names(vertex)
-362
-363    corr_data = {}
+356    files, idx = _get_files(path, filestem, idl)
+357
+358    mom_in = None
+359    mom_out = None
+360
+361    vertex_names = []
+362    for vertex in vertices:
+363        vertex_names += _get_lorentz_names(vertex)
 364
-365    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
+365    corr_data = {}
 366
-367    for hd5_file in files:
-368        file = h5py.File(path + '/' + hd5_file, "r")
-369
-370        for i in range(32):
-371            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
-372            if name in vertex_names:
-373                if name not in corr_data:
-374                    corr_data[name] = []
-375                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
-376                corr_data[name].append(raw_data)
-377                if mom_in is None:
-378                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-379                if mom_out is None:
-380                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
-381
-382        file.close()
+367    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
+368
+369    for hd5_file in files:
+370        file = h5py.File(path + '/' + hd5_file, "r")
+371
+372        for i in range(32):
+373            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
+374            if name in vertex_names:
+375                if name not in corr_data:
+376                    corr_data[name] = []
+377                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
+378                corr_data[name].append(raw_data)
+379                if mom_in is None:
+380                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+381                if mom_out is None:
+382                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
 383
-384    intermediate_dict = {}
+384        file.close()
 385
-386    for vertex in vertices:
-387        lorentz_names = _get_lorentz_names(vertex)
-388        for v_name in lorentz_names:
-389            if v_name in [('SigmaXY', 'SigmaZT'),
-390                          ('SigmaXT', 'SigmaYZ'),
-391                          ('SigmaYZ', 'SigmaXT'),
-392                          ('SigmaZT', 'SigmaXY')]:
-393                sign = -1
-394            else:
-395                sign = 1
-396            if vertex not in intermediate_dict:
-397                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
-398            else:
-399                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
-400
-401    result_dict = {}
+386    intermediate_dict = {}
+387
+388    for vertex in vertices:
+389        lorentz_names = _get_lorentz_names(vertex)
+390        for v_name in lorentz_names:
+391            if v_name in [('SigmaXY', 'SigmaZT'),
+392                          ('SigmaXT', 'SigmaYZ'),
+393                          ('SigmaYZ', 'SigmaXT'),
+394                          ('SigmaZT', 'SigmaXY')]:
+395                sign = -1
+396            else:
+397                sign = 1
+398            if vertex not in intermediate_dict:
+399                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
+400            else:
+401                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
 402
-403    for key, data in intermediate_dict.items():
+403    result_dict = {}
 404
-405        rolled_array = np.moveaxis(data, 0, 8)
+405    for key, data in intermediate_dict.items():
 406
-407        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-408        for index in np.ndindex(rolled_array.shape[:-1]):
-409            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
-410            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
-411            matrix[index] = CObs(real, imag)
-412
-413        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+407        rolled_array = np.moveaxis(data, 0, 8)
+408
+409        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+410        for index in np.ndindex(rolled_array.shape[:-1]):
+411            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
+412            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
+413            matrix[index] = CObs(real, imag)
 414
-415    return result_dict
+415        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+416
+417    return result_dict