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

 
 
@@ -1139,49 +1145,49 @@ in and out momenta of the propagator are exchanged.
- 
259def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
-260    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
-261
-262    Parameters
-263    ----------
-264    path : str
-265        path to the files to read
-266    filestem : str
-267        namestem of the files to read
-268    ens_id : str
-269        name of the ensemble, required for internal bookkeeping
-270    idl : range
-271        If specified only configurations in the given range are read in.
-272
-273    Returns
-274    -------
-275    result : Npr_matrix
-276        read Cobs-matrix
-277    """
-278
-279    files, idx = _get_files(path, filestem, idl)
-280
-281    mom = None
-282
-283    corr_data = []
-284    for hd5_file in files:
-285        file = h5py.File(path + '/' + hd5_file, "r")
-286        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
-287        corr_data.append(raw_data)
-288        if mom is None:
-289            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-290        file.close()
-291    corr_data = np.array(corr_data)
-292
-293    rolled_array = np.rollaxis(corr_data, 0, 5)
-294
-295    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-296    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
-297        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
-298        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
-299        matrix[si, sj, ci, cj] = CObs(real, imag)
-300
-301    return Npr_matrix(matrix, mom_in=mom)
+            
262def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
+263    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
+264
+265    Parameters
+266    ----------
+267    path : str
+268        path to the files to read
+269    filestem : str
+270        namestem of the files to read
+271    ens_id : str
+272        name of the ensemble, required for internal bookkeeping
+273    idl : range
+274        If specified only configurations in the given range are read in.
+275
+276    Returns
+277    -------
+278    result : Npr_matrix
+279        read Cobs-matrix
+280    """
+281
+282    files, idx = _get_files(path, filestem, idl)
+283
+284    mom = None
+285
+286    corr_data = []
+287    for hd5_file in files:
+288        file = h5py.File(path + '/' + hd5_file, "r")
+289        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
+290        corr_data.append(raw_data)
+291        if mom is None:
+292            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+293        file.close()
+294    corr_data = np.array(corr_data)
+295
+296    rolled_array = np.rollaxis(corr_data, 0, 5)
+297
+298    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+299    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
+300        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
+301        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
+302        matrix[si, sj, ci, cj] = CObs(real, imag)
+303
+304    return Npr_matrix(matrix, mom_in=mom)
 

 
 
@@ -1221,63 +1227,63 @@ read Cobs-matrix
- 
304def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
-305    """Read hadrons Bilinear hdf5 file and output an array of CObs
-306
-307    Parameters
-308    ----------
-309    path : str
-310        path to the files to read
-311    filestem : str
-312        namestem of the files to read
-313    ens_id : str
-314        name of the ensemble, required for internal bookkeeping
-315    idl : range
-316        If specified only configurations in the given range are read in.
-317
-318    Returns
-319    -------
-320    result_dict: dict[Npr_matrix]
-321        extracted Bilinears
-322    """
-323
-324    files, idx = _get_files(path, filestem, idl)
-325
-326    mom_in = None
-327    mom_out = None
+            
307def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
+308    """Read hadrons Bilinear hdf5 file and output an array of CObs
+309
+310    Parameters
+311    ----------
+312    path : str
+313        path to the files to read
+314    filestem : str
+315        namestem of the files to read
+316    ens_id : str
+317        name of the ensemble, required for internal bookkeeping
+318    idl : range
+319        If specified only configurations in the given range are read in.
+320
+321    Returns
+322    -------
+323    result_dict: dict[Npr_matrix]
+324        extracted Bilinears
+325    """
+326
+327    files, idx = _get_files(path, filestem, idl)
 328
-329    corr_data = {}
-330    for hd5_file in files:
-331        file = h5py.File(path + '/' + hd5_file, "r")
-332        for i in range(16):
-333            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
-334            if name not in corr_data:
-335                corr_data[name] = []
-336            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
-337            corr_data[name].append(raw_data)
-338            if mom_in is None:
-339                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-340            if mom_out is None:
-341                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
-342
-343        file.close()
-344
-345    result_dict = {}
-346
-347    for key, data in corr_data.items():
-348        local_data = np.array(data)
+329    mom_in = None
+330    mom_out = None
+331
+332    corr_data = {}
+333    for hd5_file in files:
+334        file = h5py.File(path + '/' + hd5_file, "r")
+335        for i in range(16):
+336            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
+337            if name not in corr_data:
+338                corr_data[name] = []
+339            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
+340            corr_data[name].append(raw_data)
+341            if mom_in is None:
+342                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+343            if mom_out is None:
+344                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
+345
+346        file.close()
+347
+348    result_dict = {}
 349
-350        rolled_array = np.rollaxis(local_data, 0, 5)
-351
-352        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-353        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
-354            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
-355            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
-356            matrix[si, sj, ci, cj] = CObs(real, imag)
-357
-358        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
-359
-360    return result_dict
+350    for key, data in corr_data.items():
+351        local_data = np.array(data)
+352
+353        rolled_array = np.rollaxis(local_data, 0, 5)
+354
+355        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+356        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
+357            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
+358            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
+359            matrix[si, sj, ci, cj] = CObs(real, imag)
+360
+361        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+362
+363    return result_dict
 

 
 
@@ -1317,90 +1323,90 @@ extracted Bilinears
- 
363def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
-364    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
-365
-366    Parameters
-367    ----------
-368    path : str
-369        path to the files to read
-370    filestem : str
-371        namestem of the files to read
-372    ens_id : str
-373        name of the ensemble, required for internal bookkeeping
-374    idl : range
-375        If specified only configurations in the given range are read in.
-376    vertices : list
-377        Vertex functions to be extracted.
-378
-379    Returns
-380    -------
-381    result_dict : dict
-382        extracted fourquark matrizes
-383    """
-384
-385    files, idx = _get_files(path, filestem, idl)
-386
-387    mom_in = None
-388    mom_out = None
+            
366def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
+367    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
+368
+369    Parameters
+370    ----------
+371    path : str
+372        path to the files to read
+373    filestem : str
+374        namestem of the files to read
+375    ens_id : str
+376        name of the ensemble, required for internal bookkeeping
+377    idl : range
+378        If specified only configurations in the given range are read in.
+379    vertices : list
+380        Vertex functions to be extracted.
+381
+382    Returns
+383    -------
+384    result_dict : dict
+385        extracted fourquark matrizes
+386    """
+387
+388    files, idx = _get_files(path, filestem, idl)
 389
-390    vertex_names = []
-391    for vertex in vertices:
-392        vertex_names += _get_lorentz_names(vertex)
-393
-394    corr_data = {}
-395
-396    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
-397
-398    for hd5_file in files:
-399        file = h5py.File(path + '/' + hd5_file, "r")
+390    mom_in = None
+391    mom_out = None
+392
+393    vertex_names = []
+394    for vertex in vertices:
+395        vertex_names += _get_lorentz_names(vertex)
+396
+397    corr_data = {}
+398
+399    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
 400
-401        for i in range(32):
-402            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
-403            if name in vertex_names:
-404                if name not in corr_data:
-405                    corr_data[name] = []
-406                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
-407                corr_data[name].append(raw_data)
-408                if mom_in is None:
-409                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-410                if mom_out is None:
-411                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
-412
-413        file.close()
-414
-415    intermediate_dict = {}
-416
-417    for vertex in vertices:
-418        lorentz_names = _get_lorentz_names(vertex)
-419        for v_name in lorentz_names:
-420            if v_name in [('SigmaXY', 'SigmaZT'),
-421                          ('SigmaXT', 'SigmaYZ'),
-422                          ('SigmaYZ', 'SigmaXT'),
-423                          ('SigmaZT', 'SigmaXY')]:
-424                sign = -1
-425            else:
-426                sign = 1
-427            if vertex not in intermediate_dict:
-428                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
-429            else:
-430                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
-431
-432    result_dict = {}
-433
-434    for key, data in intermediate_dict.items():
-435
-436        rolled_array = np.moveaxis(data, 0, 8)
-437
-438        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-439        for index in np.ndindex(rolled_array.shape[:-1]):
-440            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
-441            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
-442            matrix[index] = CObs(real, imag)
-443
-444        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
-445
-446    return result_dict
+401    for hd5_file in files:
+402        file = h5py.File(path + '/' + hd5_file, "r")
+403
+404        for i in range(32):
+405            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
+406            if name in vertex_names:
+407                if name not in corr_data:
+408                    corr_data[name] = []
+409                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
+410                corr_data[name].append(raw_data)
+411                if mom_in is None:
+412                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+413                if mom_out is None:
+414                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
+415
+416        file.close()
+417
+418    intermediate_dict = {}
+419
+420    for vertex in vertices:
+421        lorentz_names = _get_lorentz_names(vertex)
+422        for v_name in lorentz_names:
+423            if v_name in [('SigmaXY', 'SigmaZT'),
+424                          ('SigmaXT', 'SigmaYZ'),
+425                          ('SigmaYZ', 'SigmaXT'),
+426                          ('SigmaZT', 'SigmaXY')]:
+427                sign = -1
+428            else:
+429                sign = 1
+430            if vertex not in intermediate_dict:
+431                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
+432            else:
+433                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
+434
+435    result_dict = {}
+436
+437    for key, data in intermediate_dict.items():
+438
+439        rolled_array = np.moveaxis(data, 0, 8)
+440
+441        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+442        for index in np.ndindex(rolled_array.shape[:-1]):
+443            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
+444            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
+445            matrix[index] = CObs(real, imag)
+446
+447        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+448
+449    return result_dict