diff --git a/docs/pyerrors/input/openQCD.html b/docs/pyerrors/input/openQCD.html
index b8974f9d..af63fa04 100644
--- a/docs/pyerrors/input/openQCD.html
+++ b/docs/pyerrors/input/openQCD.html
@@ -108,1159 +108,1181 @@
12from ..correlators import Corr
13
14
- 15def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
- 16 """Read rwms format from given folder structure. Returns a list of length nrw
- 17
- 18 Parameters
- 19 ----------
- 20 path : str
- 21 path that contains the data files
- 22 prefix : str
- 23 all files in path that start with prefix are considered as input files.
- 24 May be used together postfix to consider only special file endings.
- 25 Prefix is ignored, if the keyword 'files' is used.
- 26 version : str
- 27 version of openQCD, default 2.0
- 28 names : list
- 29 list of names that is assigned to the data according according
- 30 to the order in the file list. Use careful, if you do not provide file names!
- 31 r_start : list
- 32 list which contains the first config to be read for each replicum
- 33 r_stop : list
- 34 list which contains the last config to be read for each replicum
- 35 r_step : int
- 36 integer that defines a fixed step size between two measurements (in units of configs)
- 37 If not given, r_step=1 is assumed.
- 38 postfix : str
- 39 postfix of the file to read, e.g. '.ms1' for openQCD-files
- 40 files : list
- 41 list which contains the filenames to be read. No automatic detection of
- 42 files performed if given.
- 43 print_err : bool
- 44 Print additional information that is useful for debugging.
- 45
- 46 Returns
- 47 -------
- 48 rwms : Obs
- 49 Reweighting factors read
- 50 """
- 51 known_oqcd_versions = ['1.4', '1.6', '2.0']
- 52 if not (version in known_oqcd_versions):
- 53 raise Exception('Unknown openQCD version defined!')
- 54 print("Working with openQCD version " + version)
- 55 if 'postfix' in kwargs:
- 56 postfix = kwargs.get('postfix')
- 57 else:
- 58 postfix = ''
- 59 ls = []
- 60 for (dirpath, dirnames, filenames) in os.walk(path):
- 61 ls.extend(filenames)
- 62 break
- 63
- 64 if not ls:
- 65 raise Exception(f"Error, directory '{path}' not found")
- 66 if 'files' in kwargs:
- 67 ls = kwargs.get('files')
- 68 else:
- 69 for exc in ls:
- 70 if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'):
- 71 ls = list(set(ls) - set([exc]))
- 72 if len(ls) > 1:
- 73 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
- 74 replica = len(ls)
- 75
- 76 if 'r_start' in kwargs:
- 77 r_start = kwargs.get('r_start')
- 78 if len(r_start) != replica:
- 79 raise Exception('r_start does not match number of replicas')
- 80 r_start = [o if o else None for o in r_start]
- 81 else:
- 82 r_start = [None] * replica
- 83
- 84 if 'r_stop' in kwargs:
- 85 r_stop = kwargs.get('r_stop')
- 86 if len(r_stop) != replica:
- 87 raise Exception('r_stop does not match number of replicas')
- 88 else:
- 89 r_stop = [None] * replica
- 90
- 91 if 'r_step' in kwargs:
- 92 r_step = kwargs.get('r_step')
- 93 else:
- 94 r_step = 1
- 95
- 96 print('Read reweighting factors from', prefix[:-1], ',',
- 97 replica, 'replica', end='')
- 98
- 99 if names is None:
- 100 rep_names = []
- 101 for entry in ls:
- 102 truncated_entry = entry
- 103 suffixes = [".dat", ".rwms", ".ms1"]
- 104 for suffix in suffixes:
- 105 if truncated_entry.endswith(suffix):
- 106 truncated_entry = truncated_entry[0:-len(suffix)]
- 107 idx = truncated_entry.index('r')
- 108 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
- 109 else:
- 110 rep_names = names
- 111
- 112 print_err = 0
- 113 if 'print_err' in kwargs:
- 114 print_err = 1
- 115 print()
+ 15def _find_files(path, prefix, postfix, ext, known_files=[]):
+ 16 found = []
+ 17 files = []
+ 18
+ 19 if postfix != "":
+ 20 if postfix[-1] != ".":
+ 21 postfix = postfix + "."
+ 22 if postfix[0] != ".":
+ 23 postfix = "." + postfix
+ 24
+ 25 if ext[0] == ".":
+ 26 ext = ext[1:]
+ 27
+ 28 pattern = prefix + "*" + postfix + ext
+ 29
+ 30 for (dirpath, dirnames, filenames) in os.walk(path + "/"):
+ 31 found.extend(filenames)
+ 32 break
+ 33
+ 34 if known_files != []:
+ 35 for kf in known_files:
+ 36 if kf not in found:
+ 37 raise FileNotFoundError("Given file " + kf + " does not exist!")
+ 38
+ 39 return known_files
+ 40
+ 41 if not found:
+ 42 raise FileNotFoundError(f"Error, directory '{path}' not found")
+ 43
+ 44 for f in found:
+ 45 if fnmatch.fnmatch(f, pattern):
+ 46 files.append(f)
+ 47
+ 48 if files == []:
+ 49 raise Exception("No files found after pattern filter!")
+ 50
+ 51 files.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
+ 52 return files
+ 53
+ 54
+ 55def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
+ 56 """Read rwms format from given folder structure. Returns a list of length nrw
+ 57
+ 58 Parameters
+ 59 ----------
+ 60 path : str
+ 61 path that contains the data files
+ 62 prefix : str
+ 63 all files in path that start with prefix are considered as input files.
+ 64 May be used together postfix to consider only special file endings.
+ 65 Prefix is ignored, if the keyword 'files' is used.
+ 66 version : str
+ 67 version of openQCD, default 2.0
+ 68 names : list
+ 69 list of names that is assigned to the data according according
+ 70 to the order in the file list. Use careful, if you do not provide file names!
+ 71 r_start : list
+ 72 list which contains the first config to be read for each replicum
+ 73 r_stop : list
+ 74 list which contains the last config to be read for each replicum
+ 75 r_step : int
+ 76 integer that defines a fixed step size between two measurements (in units of configs)
+ 77 If not given, r_step=1 is assumed.
+ 78 postfix : str
+ 79 postfix of the file to read, e.g. '.ms1' for openQCD-files
+ 80 files : list
+ 81 list which contains the filenames to be read. No automatic detection of
+ 82 files performed if given.
+ 83 print_err : bool
+ 84 Print additional information that is useful for debugging.
+ 85
+ 86 Returns
+ 87 -------
+ 88 rwms : Obs
+ 89 Reweighting factors read
+ 90 """
+ 91 known_oqcd_versions = ['1.4', '1.6', '2.0']
+ 92 if not (version in known_oqcd_versions):
+ 93 raise Exception('Unknown openQCD version defined!')
+ 94 print("Working with openQCD version " + version)
+ 95 if 'postfix' in kwargs:
+ 96 postfix = kwargs.get('postfix')
+ 97 else:
+ 98 postfix = ''
+ 99
+ 100 if 'files' in kwargs:
+ 101 known_files = kwargs.get('files')
+ 102 else:
+ 103 known_files = []
+ 104
+ 105 ls = _find_files(path, prefix, postfix, 'dat', known_files=known_files)
+ 106
+ 107 replica = len(ls)
+ 108
+ 109 if 'r_start' in kwargs:
+ 110 r_start = kwargs.get('r_start')
+ 111 if len(r_start) != replica:
+ 112 raise Exception('r_start does not match number of replicas')
+ 113 r_start = [o if o else None for o in r_start]
+ 114 else:
+ 115 r_start = [None] * replica
116
- 117 deltas = []
- 118
- 119 configlist = []
- 120 r_start_index = []
- 121 r_stop_index = []
- 122
- 123 for rep in range(replica):
- 124 tmp_array = []
- 125 with open(path + '/' + ls[rep], 'rb') as fp:
- 126
- 127 t = fp.read(4) # number of reweighting factors
- 128 if rep == 0:
- 129 nrw = struct.unpack('i', t)[0]
- 130 if version == '2.0':
- 131 nrw = int(nrw / 2)
- 132 for k in range(nrw):
- 133 deltas.append([])
- 134 else:
- 135 if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
- 136 raise Exception('Error: different number of reweighting factors for replicum', rep)
- 137
- 138 for k in range(nrw):
- 139 tmp_array.append([])
- 140
- 141 # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
- 142 nfct = []
- 143 if version in ['1.6', '2.0']:
- 144 for i in range(nrw):
- 145 t = fp.read(4)
- 146 nfct.append(struct.unpack('i', t)[0])
- 147 else:
- 148 for i in range(nrw):
- 149 nfct.append(1)
- 150
- 151 nsrc = []
- 152 for i in range(nrw):
- 153 t = fp.read(4)
- 154 nsrc.append(struct.unpack('i', t)[0])
- 155 if version == '2.0':
- 156 if not struct.unpack('i', fp.read(4))[0] == 0:
- 157 print('something is wrong!')
- 158
- 159 configlist.append([])
- 160 while True:
- 161 t = fp.read(4)
- 162 if len(t) < 4:
- 163 break
- 164 config_no = struct.unpack('i', t)[0]
- 165 configlist[-1].append(config_no)
- 166 for i in range(nrw):
- 167 if (version == '2.0'):
- 168 tmpd = _read_array_openQCD2(fp)
- 169 tmpd = _read_array_openQCD2(fp)
- 170 tmp_rw = tmpd['arr']
- 171 tmp_nfct = 1.0
- 172 for j in range(tmpd['n'][0]):
- 173 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
- 174 if print_err:
- 175 print(config_no, i, j,
- 176 np.mean(np.exp(-np.asarray(tmp_rw[j]))),
- 177 np.std(np.exp(-np.asarray(tmp_rw[j]))))
- 178 print('Sources:',
- 179 np.exp(-np.asarray(tmp_rw[j])))
- 180 print('Partial factor:', tmp_nfct)
- 181 elif version == '1.6' or version == '1.4':
- 182 tmp_nfct = 1.0
- 183 for j in range(nfct[i]):
- 184 t = fp.read(8 * nsrc[i])
- 185 t = fp.read(8 * nsrc[i])
- 186 tmp_rw = struct.unpack('d' * nsrc[i], t)
- 187 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
- 188 if print_err:
- 189 print(config_no, i, j,
- 190 np.mean(np.exp(-np.asarray(tmp_rw))),
- 191 np.std(np.exp(-np.asarray(tmp_rw))))
- 192 print('Sources:', np.exp(-np.asarray(tmp_rw)))
- 193 print('Partial factor:', tmp_nfct)
- 194 tmp_array[i].append(tmp_nfct)
- 195
- 196 diffmeas = configlist[-1][-1] - configlist[-1][-2]
- 197 configlist[-1] = [item // diffmeas for item in configlist[-1]]
- 198 if configlist[-1][0] > 1 and diffmeas > 1:
- 199 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
- 200 offset = configlist[-1][0] - 1
- 201 configlist[-1] = [item - offset for item in configlist[-1]]
- 202
- 203 if r_start[rep] is None:
- 204 r_start_index.append(0)
- 205 else:
- 206 try:
- 207 r_start_index.append(configlist[-1].index(r_start[rep]))
- 208 except ValueError:
- 209 raise Exception('Config %d not in file with range [%d, %d]' % (
- 210 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
- 211
- 212 if r_stop[rep] is None:
- 213 r_stop_index.append(len(configlist[-1]) - 1)
- 214 else:
- 215 try:
- 216 r_stop_index.append(configlist[-1].index(r_stop[rep]))
- 217 except ValueError:
- 218 raise Exception('Config %d not in file with range [%d, %d]' % (
- 219 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
- 220
- 221 for k in range(nrw):
- 222 deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
- 223
- 224 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
- 225 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
- 226 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
- 227 if np.any([step != 1 for step in stepsizes]):
- 228 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
- 229
- 230 print(',', nrw, 'reweighting factors with', nsrc, 'sources')
- 231 result = []
- 232 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
- 233
- 234 for t in range(nrw):
- 235 result.append(Obs(deltas[t], rep_names, idl=idl))
- 236 return result
- 237
- 238
- 239def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
- 240 """Extract t0 from given .ms.dat files. Returns t0 as Obs.
- 241
- 242 It is assumed that all boundary effects have
- 243 sufficiently decayed at x0=xmin.
- 244 The data around the zero crossing of t^2<E> - 0.3
- 245 is fitted with a linear function
- 246 from which the exact root is extracted.
- 247
- 248 It is assumed that one measurement is performed for each config.
- 249 If this is not the case, the resulting idl, as well as the handling
- 250 of r_start, r_stop and r_step is wrong and the user has to correct
- 251 this in the resulting observable.
- 252
- 253 Parameters
- 254 ----------
- 255 path : str
- 256 Path to .ms.dat files
- 257 prefix : str
- 258 Ensemble prefix
- 259 dtr_read : int
- 260 Determines how many trajectories should be skipped
- 261 when reading the ms.dat files.
- 262 Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
- 263 xmin : int
- 264 First timeslice where the boundary
- 265 effects have sufficiently decayed.
- 266 spatial_extent : int
- 267 spatial extent of the lattice, required for normalization.
- 268 fit_range : int
- 269 Number of data points left and right of the zero
- 270 crossing to be included in the linear fit. (Default: 5)
- 271 r_start : list
- 272 list which contains the first config to be read for each replicum.
- 273 r_stop : list
- 274 list which contains the last config to be read for each replicum.
- 275 r_step : int
- 276 integer that defines a fixed step size between two measurements (in units of configs)
- 277 If not given, r_step=1 is assumed.
- 278 plaquette : bool
- 279 If true extract the plaquette estimate of t0 instead.
- 280 names : list
- 281 list of names that is assigned to the data according according
- 282 to the order in the file list. Use careful, if you do not provide file names!
- 283 files : list
- 284 list which contains the filenames to be read. No automatic detection of
- 285 files performed if given.
- 286 plot_fit : bool
- 287 If true, the fit for the extraction of t0 is shown together with the data.
- 288 assume_thermalization : bool
- 289 If True: If the first record divided by the distance between two measurements is larger than
- 290 1, it is assumed that this is due to thermalization and the first measurement belongs
- 291 to the first config (default).
- 292 If False: The config numbers are assumed to be traj_number // difference
- 293
- 294 Returns
- 295 -------
- 296 t0 : Obs
- 297 Extracted t0
- 298 """
- 299
- 300 ls = []
- 301 for (dirpath, dirnames, filenames) in os.walk(path):
- 302 ls.extend(filenames)
- 303 break
- 304
- 305 if not ls:
- 306 raise Exception('Error, directory not found')
- 307
- 308 if 'files' in kwargs:
- 309 ls = kwargs.get('files')
- 310 else:
- 311 for exc in ls:
- 312 if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'):
- 313 ls = list(set(ls) - set([exc]))
- 314 if len(ls) > 1:
- 315 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
- 316 replica = len(ls)
- 317
- 318 if 'r_start' in kwargs:
- 319 r_start = kwargs.get('r_start')
- 320 if len(r_start) != replica:
- 321 raise Exception('r_start does not match number of replicas')
- 322 r_start = [o if o else None for o in r_start]
- 323 else:
- 324 r_start = [None] * replica
- 325
- 326 if 'r_stop' in kwargs:
- 327 r_stop = kwargs.get('r_stop')
- 328 if len(r_stop) != replica:
- 329 raise Exception('r_stop does not match number of replicas')
- 330 else:
- 331 r_stop = [None] * replica
+ 117 if 'r_stop' in kwargs:
+ 118 r_stop = kwargs.get('r_stop')
+ 119 if len(r_stop) != replica:
+ 120 raise Exception('r_stop does not match number of replicas')
+ 121 else:
+ 122 r_stop = [None] * replica
+ 123
+ 124 if 'r_step' in kwargs:
+ 125 r_step = kwargs.get('r_step')
+ 126 else:
+ 127 r_step = 1
+ 128
+ 129 print('Read reweighting factors from', prefix[:-1], ',',
+ 130 replica, 'replica', end='')
+ 131
+ 132 if names is None:
+ 133 rep_names = []
+ 134 for entry in ls:
+ 135 truncated_entry = entry
+ 136 suffixes = [".dat", ".rwms", ".ms1"]
+ 137 for suffix in suffixes:
+ 138 if truncated_entry.endswith(suffix):
+ 139 truncated_entry = truncated_entry[0:-len(suffix)]
+ 140 idx = truncated_entry.index('r')
+ 141 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
+ 142 else:
+ 143 rep_names = names
+ 144
+ 145 print_err = 0
+ 146 if 'print_err' in kwargs:
+ 147 print_err = 1
+ 148 print()
+ 149
+ 150 deltas = []
+ 151
+ 152 configlist = []
+ 153 r_start_index = []
+ 154 r_stop_index = []
+ 155
+ 156 for rep in range(replica):
+ 157 tmp_array = []
+ 158 with open(path + '/' + ls[rep], 'rb') as fp:
+ 159
+ 160 t = fp.read(4) # number of reweighting factors
+ 161 if rep == 0:
+ 162 nrw = struct.unpack('i', t)[0]
+ 163 if version == '2.0':
+ 164 nrw = int(nrw / 2)
+ 165 for k in range(nrw):
+ 166 deltas.append([])
+ 167 else:
+ 168 if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
+ 169 raise Exception('Error: different number of reweighting factors for replicum', rep)
+ 170
+ 171 for k in range(nrw):
+ 172 tmp_array.append([])
+ 173
+ 174 # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
+ 175 nfct = []
+ 176 if version in ['1.6', '2.0']:
+ 177 for i in range(nrw):
+ 178 t = fp.read(4)
+ 179 nfct.append(struct.unpack('i', t)[0])
+ 180 else:
+ 181 for i in range(nrw):
+ 182 nfct.append(1)
+ 183
+ 184 nsrc = []
+ 185 for i in range(nrw):
+ 186 t = fp.read(4)
+ 187 nsrc.append(struct.unpack('i', t)[0])
+ 188 if version == '2.0':
+ 189 if not struct.unpack('i', fp.read(4))[0] == 0:
+ 190 raise Exception("You are using the input for openQCD version 2.0, this is not correct.")
+ 191
+ 192 configlist.append([])
+ 193 while True:
+ 194 t = fp.read(4)
+ 195 if len(t) < 4:
+ 196 break
+ 197 config_no = struct.unpack('i', t)[0]
+ 198 configlist[-1].append(config_no)
+ 199 for i in range(nrw):
+ 200 if (version == '2.0'):
+ 201 tmpd = _read_array_openQCD2(fp)
+ 202 tmpd = _read_array_openQCD2(fp)
+ 203 tmp_rw = tmpd['arr']
+ 204 tmp_nfct = 1.0
+ 205 for j in range(tmpd['n'][0]):
+ 206 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
+ 207 if print_err:
+ 208 print(config_no, i, j,
+ 209 np.mean(np.exp(-np.asarray(tmp_rw[j]))),
+ 210 np.std(np.exp(-np.asarray(tmp_rw[j]))))
+ 211 print('Sources:',
+ 212 np.exp(-np.asarray(tmp_rw[j])))
+ 213 print('Partial factor:', tmp_nfct)
+ 214 elif version == '1.6' or version == '1.4':
+ 215 tmp_nfct = 1.0
+ 216 for j in range(nfct[i]):
+ 217 t = fp.read(8 * nsrc[i])
+ 218 t = fp.read(8 * nsrc[i])
+ 219 tmp_rw = struct.unpack('d' * nsrc[i], t)
+ 220 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
+ 221 if print_err:
+ 222 print(config_no, i, j,
+ 223 np.mean(np.exp(-np.asarray(tmp_rw))),
+ 224 np.std(np.exp(-np.asarray(tmp_rw))))
+ 225 print('Sources:', np.exp(-np.asarray(tmp_rw)))
+ 226 print('Partial factor:', tmp_nfct)
+ 227 tmp_array[i].append(tmp_nfct)
+ 228
+ 229 diffmeas = configlist[-1][-1] - configlist[-1][-2]
+ 230 configlist[-1] = [item // diffmeas for item in configlist[-1]]
+ 231 if configlist[-1][0] > 1 and diffmeas > 1:
+ 232 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
+ 233 offset = configlist[-1][0] - 1
+ 234 configlist[-1] = [item - offset for item in configlist[-1]]
+ 235
+ 236 if r_start[rep] is None:
+ 237 r_start_index.append(0)
+ 238 else:
+ 239 try:
+ 240 r_start_index.append(configlist[-1].index(r_start[rep]))
+ 241 except ValueError:
+ 242 raise Exception('Config %d not in file with range [%d, %d]' % (
+ 243 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
+ 244
+ 245 if r_stop[rep] is None:
+ 246 r_stop_index.append(len(configlist[-1]) - 1)
+ 247 else:
+ 248 try:
+ 249 r_stop_index.append(configlist[-1].index(r_stop[rep]))
+ 250 except ValueError:
+ 251 raise Exception('Config %d not in file with range [%d, %d]' % (
+ 252 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
+ 253
+ 254 for k in range(nrw):
+ 255 deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
+ 256
+ 257 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
+ 258 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
+ 259 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
+ 260 if np.any([step != 1 for step in stepsizes]):
+ 261 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
+ 262
+ 263 print(',', nrw, 'reweighting factors with', nsrc, 'sources')
+ 264 result = []
+ 265 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
+ 266
+ 267 for t in range(nrw):
+ 268 result.append(Obs(deltas[t], rep_names, idl=idl))
+ 269 return result
+ 270
+ 271
+ 272def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
+ 273 """Extract t0 from given .ms.dat files. Returns t0 as Obs.
+ 274
+ 275 It is assumed that all boundary effects have
+ 276 sufficiently decayed at x0=xmin.
+ 277 The data around the zero crossing of t^2<E> - 0.3
+ 278 is fitted with a linear function
+ 279 from which the exact root is extracted.
+ 280
+ 281 It is assumed that one measurement is performed for each config.
+ 282 If this is not the case, the resulting idl, as well as the handling
+ 283 of r_start, r_stop and r_step is wrong and the user has to correct
+ 284 this in the resulting observable.
+ 285
+ 286 Parameters
+ 287 ----------
+ 288 path : str
+ 289 Path to .ms.dat files
+ 290 prefix : str
+ 291 Ensemble prefix
+ 292 dtr_read : int
+ 293 Determines how many trajectories should be skipped
+ 294 when reading the ms.dat files.
+ 295 Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
+ 296 xmin : int
+ 297 First timeslice where the boundary
+ 298 effects have sufficiently decayed.
+ 299 spatial_extent : int
+ 300 spatial extent of the lattice, required for normalization.
+ 301 fit_range : int
+ 302 Number of data points left and right of the zero
+ 303 crossing to be included in the linear fit. (Default: 5)
+ 304 r_start : list
+ 305 list which contains the first config to be read for each replicum.
+ 306 r_stop : list
+ 307 list which contains the last config to be read for each replicum.
+ 308 r_step : int
+ 309 integer that defines a fixed step size between two measurements (in units of configs)
+ 310 If not given, r_step=1 is assumed.
+ 311 plaquette : bool
+ 312 If true extract the plaquette estimate of t0 instead.
+ 313 names : list
+ 314 list of names that is assigned to the data according according
+ 315 to the order in the file list. Use careful, if you do not provide file names!
+ 316 files : list
+ 317 list which contains the filenames to be read. No automatic detection of
+ 318 files performed if given.
+ 319 plot_fit : bool
+ 320 If true, the fit for the extraction of t0 is shown together with the data.
+ 321 assume_thermalization : bool
+ 322 If True: If the first record divided by the distance between two measurements is larger than
+ 323 1, it is assumed that this is due to thermalization and the first measurement belongs
+ 324 to the first config (default).
+ 325 If False: The config numbers are assumed to be traj_number // difference
+ 326
+ 327 Returns
+ 328 -------
+ 329 t0 : Obs
+ 330 Extracted t0
+ 331 """
332
- 333 if 'r_step' in kwargs:
- 334 r_step = kwargs.get('r_step')
+ 333 if 'files' in kwargs:
+ 334 known_files = kwargs.get('files')
335 else:
- 336 r_step = 1
+ 336 known_files = []
337
- 338 print('Extract t0 from', prefix, ',', replica, 'replica')
+ 338 ls = _find_files(path, prefix, 'ms', 'dat', known_files=known_files)
339
- 340 if 'names' in kwargs:
- 341 rep_names = kwargs.get('names')
- 342 else:
- 343 rep_names = []
- 344 for entry in ls:
- 345 truncated_entry = entry.split('.')[0]
- 346 idx = truncated_entry.index('r')
- 347 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
- 348
- 349 Ysum = []
- 350
- 351 configlist = []
- 352 r_start_index = []
- 353 r_stop_index = []
- 354
- 355 for rep in range(replica):
+ 340 replica = len(ls)
+ 341
+ 342 if 'r_start' in kwargs:
+ 343 r_start = kwargs.get('r_start')
+ 344 if len(r_start) != replica:
+ 345 raise Exception('r_start does not match number of replicas')
+ 346 r_start = [o if o else None for o in r_start]
+ 347 else:
+ 348 r_start = [None] * replica
+ 349
+ 350 if 'r_stop' in kwargs:
+ 351 r_stop = kwargs.get('r_stop')
+ 352 if len(r_stop) != replica:
+ 353 raise Exception('r_stop does not match number of replicas')
+ 354 else:
+ 355 r_stop = [None] * replica
356
- 357 with open(path + '/' + ls[rep], 'rb') as fp:
- 358 t = fp.read(12)
- 359 header = struct.unpack('iii', t)
- 360 if rep == 0:
- 361 dn = header[0]
- 362 nn = header[1]
- 363 tmax = header[2]
- 364 elif dn != header[0] or nn != header[1] or tmax != header[2]:
- 365 raise Exception('Replica parameters do not match.')
- 366
- 367 t = fp.read(8)
- 368 if rep == 0:
- 369 eps = struct.unpack('d', t)[0]
- 370 print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
- 371 elif eps != struct.unpack('d', t)[0]:
- 372 raise Exception('Values for eps do not match among replica.')
- 373
- 374 Ysl = []
- 375
- 376 configlist.append([])
- 377 while True:
- 378 t = fp.read(4)
- 379 if (len(t) < 4):
- 380 break
- 381 nc = struct.unpack('i', t)[0]
- 382 configlist[-1].append(nc)
- 383
- 384 t = fp.read(8 * tmax * (nn + 1))
- 385 if kwargs.get('plaquette'):
- 386 if nc % dtr_read == 0:
- 387 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
- 388 t = fp.read(8 * tmax * (nn + 1))
- 389 if not kwargs.get('plaquette'):
- 390 if nc % dtr_read == 0:
- 391 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
- 392 t = fp.read(8 * tmax * (nn + 1))
- 393
- 394 Ysum.append([])
- 395 for i, item in enumerate(Ysl):
- 396 Ysum[-1].append([np.mean(item[current + xmin:
- 397 current + tmax - xmin])
- 398 for current in range(0, len(item), tmax)])
+ 357 if 'r_step' in kwargs:
+ 358 r_step = kwargs.get('r_step')
+ 359 else:
+ 360 r_step = 1
+ 361
+ 362 print('Extract t0 from', prefix, ',', replica, 'replica')
+ 363
+ 364 if 'names' in kwargs:
+ 365 rep_names = kwargs.get('names')
+ 366 else:
+ 367 rep_names = []
+ 368 for entry in ls:
+ 369 truncated_entry = entry.split('.')[0]
+ 370 idx = truncated_entry.index('r')
+ 371 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
+ 372
+ 373 Ysum = []
+ 374
+ 375 configlist = []
+ 376 r_start_index = []
+ 377 r_stop_index = []
+ 378
+ 379 for rep in range(replica):
+ 380
+ 381 with open(path + '/' + ls[rep], 'rb') as fp:
+ 382 t = fp.read(12)
+ 383 header = struct.unpack('iii', t)
+ 384 if rep == 0:
+ 385 dn = header[0]
+ 386 nn = header[1]
+ 387 tmax = header[2]
+ 388 elif dn != header[0] or nn != header[1] or tmax != header[2]:
+ 389 raise Exception('Replica parameters do not match.')
+ 390
+ 391 t = fp.read(8)
+ 392 if rep == 0:
+ 393 eps = struct.unpack('d', t)[0]
+ 394 print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
+ 395 elif eps != struct.unpack('d', t)[0]:
+ 396 raise Exception('Values for eps do not match among replica.')
+ 397
+ 398 Ysl = []
399
- 400 diffmeas = configlist[-1][-1] - configlist[-1][-2]
- 401 configlist[-1] = [item // diffmeas for item in configlist[-1]]
- 402 if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
- 403 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
- 404 offset = configlist[-1][0] - 1
- 405 configlist[-1] = [item - offset for item in configlist[-1]]
- 406
- 407 if r_start[rep] is None:
- 408 r_start_index.append(0)
- 409 else:
- 410 try:
- 411 r_start_index.append(configlist[-1].index(r_start[rep]))
- 412 except ValueError:
- 413 raise Exception('Config %d not in file with range [%d, %d]' % (
- 414 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
- 415
- 416 if r_stop[rep] is None:
- 417 r_stop_index.append(len(configlist[-1]) - 1)
- 418 else:
- 419 try:
- 420 r_stop_index.append(configlist[-1].index(r_stop[rep]))
- 421 except ValueError:
- 422 raise Exception('Config %d not in file with range [%d, %d]' % (
- 423 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
- 424
- 425 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
- 426 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
- 427 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
- 428 if np.any([step != 1 for step in stepsizes]):
- 429 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
+ 400 configlist.append([])
+ 401 while True:
+ 402 t = fp.read(4)
+ 403 if (len(t) < 4):
+ 404 break
+ 405 nc = struct.unpack('i', t)[0]
+ 406 configlist[-1].append(nc)
+ 407
+ 408 t = fp.read(8 * tmax * (nn + 1))
+ 409 if kwargs.get('plaquette'):
+ 410 if nc % dtr_read == 0:
+ 411 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
+ 412 t = fp.read(8 * tmax * (nn + 1))
+ 413 if not kwargs.get('plaquette'):
+ 414 if nc % dtr_read == 0:
+ 415 Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
+ 416 t = fp.read(8 * tmax * (nn + 1))
+ 417
+ 418 Ysum.append([])
+ 419 for i, item in enumerate(Ysl):
+ 420 Ysum[-1].append([np.mean(item[current + xmin:
+ 421 current + tmax - xmin])
+ 422 for current in range(0, len(item), tmax)])
+ 423
+ 424 diffmeas = configlist[-1][-1] - configlist[-1][-2]
+ 425 configlist[-1] = [item // diffmeas for item in configlist[-1]]
+ 426 if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
+ 427 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
+ 428 offset = configlist[-1][0] - 1
+ 429 configlist[-1] = [item - offset for item in configlist[-1]]
430
- 431 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
- 432 t2E_dict = {}
- 433 for n in range(nn + 1):
- 434 samples = []
- 435 for nrep, rep in enumerate(Ysum):
- 436 samples.append([])
- 437 for cnfg in rep:
- 438 samples[-1].append(cnfg[n])
- 439 samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
- 440 new_obs = Obs(samples, rep_names, idl=idl)
- 441 t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
- 442
- 443 zero_crossing = np.argmax(np.array(
- 444 [o.value for o in t2E_dict.values()]) > 0.0)
- 445
- 446 x = list(t2E_dict.keys())[zero_crossing - fit_range:
- 447 zero_crossing + fit_range]
- 448 y = list(t2E_dict.values())[zero_crossing - fit_range:
- 449 zero_crossing + fit_range]
- 450 [o.gamma_method() for o in y]
- 451
- 452 fit_result = fit_lin(x, y)
- 453
- 454 if kwargs.get('plot_fit'):
- 455 plt.figure()
- 456 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
- 457 ax0 = plt.subplot(gs[0])
- 458 xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
- 459 ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
- 460 [o.gamma_method() for o in ymore]
- 461 ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
- 462 xplot = np.linspace(np.min(x), np.max(x))
- 463 yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
- 464 [yi.gamma_method() for yi in yplot]
- 465 ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
- 466 retval = (-fit_result[0] / fit_result[1])
- 467 retval.gamma_method()
- 468 ylim = ax0.get_ylim()
- 469 ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
- 470 ax0.set_ylim(ylim)
- 471 ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
- 472 xlim = ax0.get_xlim()
- 473
- 474 fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
- 475 residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
- 476 ax1 = plt.subplot(gs[1])
- 477 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
- 478 ax1.tick_params(direction='out')
- 479 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
- 480 ax1.axhline(y=0.0, ls='--', color='k')
- 481 ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
- 482 ax1.set_xlim(xlim)
- 483 ax1.set_ylabel('Residuals')
- 484 ax1.set_xlabel(r'$t/a^2$')
- 485
- 486 plt.draw()
- 487 return -fit_result[0] / fit_result[1]
- 488
- 489
- 490def _parse_array_openQCD2(d, n, size, wa, quadrupel=False):
- 491 arr = []
- 492 if d == 2:
- 493 for i in range(n[0]):
- 494 tmp = wa[i * n[1]:(i + 1) * n[1]]
- 495 if quadrupel:
- 496 tmp2 = []
- 497 for j in range(0, len(tmp), 2):
- 498 tmp2.append(tmp[j])
- 499 arr.append(tmp2)
- 500 else:
- 501 arr.append(np.asarray(tmp))
- 502
- 503 else:
- 504 raise Exception('Only two-dimensional arrays supported!')
- 505
- 506 return arr
- 507
- 508
- 509def _read_array_openQCD2(fp):
- 510 t = fp.read(4)
- 511 d = struct.unpack('i', t)[0]
- 512 t = fp.read(4 * d)
- 513 n = struct.unpack('%di' % (d), t)
- 514 t = fp.read(4)
- 515 size = struct.unpack('i', t)[0]
- 516 if size == 4:
- 517 types = 'i'
- 518 elif size == 8:
- 519 types = 'd'
- 520 elif size == 16:
- 521 types = 'dd'
- 522 else:
- 523 raise Exception("Type for size '" + str(size) + "' not known.")
- 524 m = n[0]
- 525 for i in range(1, d):
- 526 m *= n[i]
- 527
- 528 t = fp.read(m * size)
- 529 tmp = struct.unpack('%d%s' % (m, types), t)
- 530
- 531 arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True)
- 532 return {'d': d, 'n': n, 'size': size, 'arr': arr}
- 533
- 534
- 535def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
- 536 """Read the topologial charge based on openQCD gradient flow measurements.
- 537
- 538 Parameters
- 539 ----------
- 540 path : str
- 541 path of the measurement files
- 542 prefix : str
- 543 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
- 544 Ignored if file names are passed explicitly via keyword files.
- 545 c : double
- 546 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
- 547 dtr_cnfg : int
- 548 (optional) parameter that specifies the number of measurements
- 549 between two configs.
- 550 If it is not set, the distance between two measurements
- 551 in the file is assumed to be the distance between two configurations.
- 552 steps : int
- 553 (optional) Distance between two configurations in units of trajectories /
- 554 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- 555 version : str
- 556 Either openQCD or sfqcd, depending on the data.
- 557 L : int
- 558 spatial length of the lattice in L/a.
- 559 HAS to be set if version != sfqcd, since openQCD does not provide
- 560 this in the header
- 561 r_start : list
- 562 list which contains the first config to be read for each replicum.
- 563 r_stop : list
- 564 list which contains the last config to be read for each replicum.
- 565 files : list
- 566 specify the exact files that need to be read
- 567 from path, practical if e.g. only one replicum is needed
- 568 postfix : str
- 569 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
- 570 names : list
- 571 Alternative labeling for replicas/ensembles.
- 572 Has to have the appropriate length.
- 573 Zeuthen_flow : bool
- 574 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
- 575 for version=='sfqcd' If False, the Wilson flow is used.
- 576 integer_charge : bool
- 577 If True, the charge is rounded towards the nearest integer on each config.
- 578
- 579 Returns
- 580 -------
- 581 result : Obs
- 582 Read topological charge
- 583 """
- 584
- 585 return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
- 586
- 587
- 588def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
- 589 """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
- 590
- 591 Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step.
- 592
- 593 Parameters
- 594 ----------
- 595 path : str
- 596 path of the measurement files
- 597 prefix : str
- 598 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
- 599 Ignored if file names are passed explicitly via keyword files.
- 600 c : double
- 601 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
- 602 dtr_cnfg : int
- 603 (optional) parameter that specifies the number of measurements
- 604 between two configs.
- 605 If it is not set, the distance between two measurements
- 606 in the file is assumed to be the distance between two configurations.
- 607 steps : int
- 608 (optional) Distance between two configurations in units of trajectories /
- 609 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- 610 r_start : list
- 611 list which contains the first config to be read for each replicum.
- 612 r_stop : list
- 613 list which contains the last config to be read for each replicum.
- 614 files : list
- 615 specify the exact files that need to be read
- 616 from path, practical if e.g. only one replicum is needed
- 617 names : list
- 618 Alternative labeling for replicas/ensembles.
- 619 Has to have the appropriate length.
- 620 postfix : str
- 621 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
- 622 Zeuthen_flow : bool
- 623 (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
- 624 """
- 625
- 626 if c != 0.3:
- 627 raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
- 628
- 629 plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
- 630 C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
- 631 L = plaq.tag["L"]
- 632 T = plaq.tag["T"]
- 633
- 634 if T != L:
- 635 raise Exception("The required lattice norm is only implemented for T=L at the moment.")
- 636
- 637 if Zeuthen_flow is not True:
- 638 raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
- 639
- 640 t = (c * L) ** 2 / 8
- 641
- 642 normdict = {4: 0.012341170468270,
- 643 6: 0.010162691462430,
- 644 8: 0.009031614807931,
- 645 10: 0.008744966371393,
- 646 12: 0.008650917856809,
- 647 14: 8.611154391267955E-03,
- 648 16: 0.008591758449508,
- 649 20: 0.008575359627103,
- 650 24: 0.008569387847540,
- 651 28: 8.566803713382559E-03,
- 652 32: 0.008565541650006,
- 653 40: 8.564480684962046E-03,
- 654 48: 8.564098025073460E-03,
- 655 64: 8.563853943383087E-03}
- 656
- 657 return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
- 658
- 659
- 660def _read_flow_obs(path, prefix, c, dtr_cnfg=1, version="openQCD", obspos=0, sum_t=True, **kwargs):
- 661 """Read a flow observable based on openQCD gradient flow measurements.
- 662
- 663 Parameters
- 664 ----------
- 665 path : str
- 666 path of the measurement files
- 667 prefix : str
- 668 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
- 669 Ignored if file names are passed explicitly via keyword files.
- 670 c : double
- 671 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
- 672 dtr_cnfg : int
- 673 (optional) parameter that specifies the number of measurements
- 674 between two configs.
- 675 If it is not set, the distance between two measurements
- 676 in the file is assumed to be the distance between two configurations.
- 677 steps : int
- 678 (optional) Distance between two configurations in units of trajectories /
- 679 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- 680 version : str
- 681 Either openQCD or sfqcd, depending on the data.
- 682 obspos : int
- 683 position of the obeservable in the measurement file. Only relevant for sfqcd files.
- 684 sum_t : bool
- 685 If true sum over all timeslices, if false only take the value at T/2.
- 686 L : int
- 687 spatial length of the lattice in L/a.
- 688 HAS to be set if version != sfqcd, since openQCD does not provide
- 689 this in the header
- 690 r_start : list
- 691 list which contains the first config to be read for each replicum.
- 692 r_stop : list
- 693 list which contains the last config to be read for each replicum.
- 694 files : list
- 695 specify the exact files that need to be read
- 696 from path, practical if e.g. only one replicum is needed
- 697 names : list
- 698 Alternative labeling for replicas/ensembles.
- 699 Has to have the appropriate length.
- 700 postfix : str
- 701 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
- 702 Zeuthen_flow : bool
- 703 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
- 704 for version=='sfqcd' If False, the Wilson flow is used.
- 705 integer_charge : bool
- 706 If True, the charge is rounded towards the nearest integer on each config.
- 707
- 708 Returns
- 709 -------
- 710 result : Obs
- 711 flow observable specified
- 712 """
- 713 known_versions = ["openQCD", "sfqcd"]
- 714
- 715 if version not in known_versions:
- 716 raise Exception("Unknown openQCD version.")
- 717 if "steps" in kwargs:
- 718 steps = kwargs.get("steps")
- 719 if version == "sfqcd":
- 720 if "L" in kwargs:
- 721 supposed_L = kwargs.get("L")
- 722 else:
- 723 supposed_L = None
- 724 postfix = ".gfms.dat"
- 725 else:
- 726 if "L" not in kwargs:
- 727 raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
- 728 else:
- 729 L = kwargs.get("L")
- 730 postfix = ".ms.dat"
+ 431 if r_start[rep] is None:
+ 432 r_start_index.append(0)
+ 433 else:
+ 434 try:
+ 435 r_start_index.append(configlist[-1].index(r_start[rep]))
+ 436 except ValueError:
+ 437 raise Exception('Config %d not in file with range [%d, %d]' % (
+ 438 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
+ 439
+ 440 if r_stop[rep] is None:
+ 441 r_stop_index.append(len(configlist[-1]) - 1)
+ 442 else:
+ 443 try:
+ 444 r_stop_index.append(configlist[-1].index(r_stop[rep]))
+ 445 except ValueError:
+ 446 raise Exception('Config %d not in file with range [%d, %d]' % (
+ 447 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
+ 448
+ 449 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
+ 450 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
+ 451 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
+ 452 if np.any([step != 1 for step in stepsizes]):
+ 453 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
+ 454
+ 455 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
+ 456 t2E_dict = {}
+ 457 for n in range(nn + 1):
+ 458 samples = []
+ 459 for nrep, rep in enumerate(Ysum):
+ 460 samples.append([])
+ 461 for cnfg in rep:
+ 462 samples[-1].append(cnfg[n])
+ 463 samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
+ 464 new_obs = Obs(samples, rep_names, idl=idl)
+ 465 t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
+ 466
+ 467 zero_crossing = np.argmax(np.array(
+ 468 [o.value for o in t2E_dict.values()]) > 0.0)
+ 469
+ 470 x = list(t2E_dict.keys())[zero_crossing - fit_range:
+ 471 zero_crossing + fit_range]
+ 472 y = list(t2E_dict.values())[zero_crossing - fit_range:
+ 473 zero_crossing + fit_range]
+ 474 [o.gamma_method() for o in y]
+ 475
+ 476 fit_result = fit_lin(x, y)
+ 477
+ 478 if kwargs.get('plot_fit'):
+ 479 plt.figure()
+ 480 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
+ 481 ax0 = plt.subplot(gs[0])
+ 482 xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
+ 483 ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
+ 484 [o.gamma_method() for o in ymore]
+ 485 ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
+ 486 xplot = np.linspace(np.min(x), np.max(x))
+ 487 yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
+ 488 [yi.gamma_method() for yi in yplot]
+ 489 ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
+ 490 retval = (-fit_result[0] / fit_result[1])
+ 491 retval.gamma_method()
+ 492 ylim = ax0.get_ylim()
+ 493 ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
+ 494 ax0.set_ylim(ylim)
+ 495 ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
+ 496 xlim = ax0.get_xlim()
+ 497
+ 498 fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
+ 499 residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
+ 500 ax1 = plt.subplot(gs[1])
+ 501 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
+ 502 ax1.tick_params(direction='out')
+ 503 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
+ 504 ax1.axhline(y=0.0, ls='--', color='k')
+ 505 ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
+ 506 ax1.set_xlim(xlim)
+ 507 ax1.set_ylabel('Residuals')
+ 508 ax1.set_xlabel(r'$t/a^2$')
+ 509
+ 510 plt.draw()
+ 511 return -fit_result[0] / fit_result[1]
+ 512
+ 513
+ 514def _parse_array_openQCD2(d, n, size, wa, quadrupel=False):
+ 515 arr = []
+ 516 if d == 2:
+ 517 for i in range(n[0]):
+ 518 tmp = wa[i * n[1]:(i + 1) * n[1]]
+ 519 if quadrupel:
+ 520 tmp2 = []
+ 521 for j in range(0, len(tmp), 2):
+ 522 tmp2.append(tmp[j])
+ 523 arr.append(tmp2)
+ 524 else:
+ 525 arr.append(np.asarray(tmp))
+ 526
+ 527 else:
+ 528 raise Exception('Only two-dimensional arrays supported!')
+ 529
+ 530 return arr
+ 531
+ 532
+ 533def _read_array_openQCD2(fp):
+ 534 t = fp.read(4)
+ 535 d = struct.unpack('i', t)[0]
+ 536 t = fp.read(4 * d)
+ 537 n = struct.unpack('%di' % (d), t)
+ 538 t = fp.read(4)
+ 539 size = struct.unpack('i', t)[0]
+ 540 if size == 4:
+ 541 types = 'i'
+ 542 elif size == 8:
+ 543 types = 'd'
+ 544 elif size == 16:
+ 545 types = 'dd'
+ 546 else:
+ 547 raise Exception("Type for size '" + str(size) + "' not known.")
+ 548 m = n[0]
+ 549 for i in range(1, d):
+ 550 m *= n[i]
+ 551
+ 552 t = fp.read(m * size)
+ 553 tmp = struct.unpack('%d%s' % (m, types), t)
+ 554
+ 555 arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True)
+ 556 return {'d': d, 'n': n, 'size': size, 'arr': arr}
+ 557
+ 558
+ 559def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
+ 560 """Read the topologial charge based on openQCD gradient flow measurements.
+ 561
+ 562 Parameters
+ 563 ----------
+ 564 path : str
+ 565 path of the measurement files
+ 566 prefix : str
+ 567 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
+ 568 Ignored if file names are passed explicitly via keyword files.
+ 569 c : double
+ 570 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
+ 571 dtr_cnfg : int
+ 572 (optional) parameter that specifies the number of measurements
+ 573 between two configs.
+ 574 If it is not set, the distance between two measurements
+ 575 in the file is assumed to be the distance between two configurations.
+ 576 steps : int
+ 577 (optional) Distance between two configurations in units of trajectories /
+ 578 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+ 579 version : str
+ 580 Either openQCD or sfqcd, depending on the data.
+ 581 L : int
+ 582 spatial length of the lattice in L/a.
+ 583 HAS to be set if version != sfqcd, since openQCD does not provide
+ 584 this in the header
+ 585 r_start : list
+ 586 list which contains the first config to be read for each replicum.
+ 587 r_stop : list
+ 588 list which contains the last config to be read for each replicum.
+ 589 files : list
+ 590 specify the exact files that need to be read
+ 591 from path, practical if e.g. only one replicum is needed
+ 592 postfix : str
+ 593 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+ 594 names : list
+ 595 Alternative labeling for replicas/ensembles.
+ 596 Has to have the appropriate length.
+ 597 Zeuthen_flow : bool
+ 598 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+ 599 for version=='sfqcd' If False, the Wilson flow is used.
+ 600 integer_charge : bool
+ 601 If True, the charge is rounded towards the nearest integer on each config.
+ 602
+ 603 Returns
+ 604 -------
+ 605 result : Obs
+ 606 Read topological charge
+ 607 """
+ 608
+ 609 return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
+ 610
+ 611
+ 612def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
+ 613 """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
+ 614
+ 615 Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step.
+ 616
+ 617 Parameters
+ 618 ----------
+ 619 path : str
+ 620 path of the measurement files
+ 621 prefix : str
+ 622 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
+ 623 Ignored if file names are passed explicitly via keyword files.
+ 624 c : double
+ 625 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
+ 626 dtr_cnfg : int
+ 627 (optional) parameter that specifies the number of measurements
+ 628 between two configs.
+ 629 If it is not set, the distance between two measurements
+ 630 in the file is assumed to be the distance between two configurations.
+ 631 steps : int
+ 632 (optional) Distance between two configurations in units of trajectories /
+ 633 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+ 634 r_start : list
+ 635 list which contains the first config to be read for each replicum.
+ 636 r_stop : list
+ 637 list which contains the last config to be read for each replicum.
+ 638 files : list
+ 639 specify the exact files that need to be read
+ 640 from path, practical if e.g. only one replicum is needed
+ 641 names : list
+ 642 Alternative labeling for replicas/ensembles.
+ 643 Has to have the appropriate length.
+ 644 postfix : str
+ 645 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+ 646 Zeuthen_flow : bool
+ 647 (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
+ 648 """
+ 649
+ 650 if c != 0.3:
+ 651 raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
+ 652
+ 653 plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
+ 654 C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
+ 655 L = plaq.tag["L"]
+ 656 T = plaq.tag["T"]
+ 657
+ 658 if T != L:
+ 659 raise Exception("The required lattice norm is only implemented for T=L at the moment.")
+ 660
+ 661 if Zeuthen_flow is not True:
+ 662 raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
+ 663
+ 664 t = (c * L) ** 2 / 8
+ 665
+ 666 normdict = {4: 0.012341170468270,
+ 667 6: 0.010162691462430,
+ 668 8: 0.009031614807931,
+ 669 10: 0.008744966371393,
+ 670 12: 0.008650917856809,
+ 671 14: 8.611154391267955E-03,
+ 672 16: 0.008591758449508,
+ 673 20: 0.008575359627103,
+ 674 24: 0.008569387847540,
+ 675 28: 8.566803713382559E-03,
+ 676 32: 0.008565541650006,
+ 677 40: 8.564480684962046E-03,
+ 678 48: 8.564098025073460E-03,
+ 679 64: 8.563853943383087E-03}
+ 680
+ 681 return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
+ 682
+ 683
+ 684def _read_flow_obs(path, prefix, c, dtr_cnfg=1, version="openQCD", obspos=0, sum_t=True, **kwargs):
+ 685 """Read a flow observable based on openQCD gradient flow measurements.
+ 686
+ 687 Parameters
+ 688 ----------
+ 689 path : str
+ 690 path of the measurement files
+ 691 prefix : str
+ 692 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
+ 693 Ignored if file names are passed explicitly via keyword files.
+ 694 c : double
+ 695 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
+ 696 dtr_cnfg : int
+ 697 (optional) parameter that specifies the number of measurements
+ 698 between two configs.
+ 699 If it is not set, the distance between two measurements
+ 700 in the file is assumed to be the distance between two configurations.
+ 701 steps : int
+ 702 (optional) Distance between two configurations in units of trajectories /
+ 703 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+ 704 version : str
+ 705 Either openQCD or sfqcd, depending on the data.
+ 706 obspos : int
+ 707 position of the obeservable in the measurement file. Only relevant for sfqcd files.
+ 708 sum_t : bool
+ 709 If true sum over all timeslices, if false only take the value at T/2.
+ 710 L : int
+ 711 spatial length of the lattice in L/a.
+ 712 HAS to be set if version != sfqcd, since openQCD does not provide
+ 713 this in the header
+ 714 r_start : list
+ 715 list which contains the first config to be read for each replicum.
+ 716 r_stop : list
+ 717 list which contains the last config to be read for each replicum.
+ 718 files : list
+ 719 specify the exact files that need to be read
+ 720 from path, practical if e.g. only one replicum is needed
+ 721 names : list
+ 722 Alternative labeling for replicas/ensembles.
+ 723 Has to have the appropriate length.
+ 724 postfix : str
+ 725 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+ 726 Zeuthen_flow : bool
+ 727 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+ 728 for version=='sfqcd' If False, the Wilson flow is used.
+ 729 integer_charge : bool
+ 730 If True, the charge is rounded towards the nearest integer on each config.
731
- 732 if "postfix" in kwargs:
- 733 postfix = kwargs.get("postfix")
- 734
- 735 if "files" in kwargs:
- 736 files = kwargs.get("files")
- 737 postfix = ''
- 738 else:
- 739 found = []
- 740 files = []
- 741 for (dirpath, dirnames, filenames) in os.walk(path + "/"):
- 742 found.extend(filenames)
- 743 break
- 744 for f in found:
- 745 if fnmatch.fnmatch(f, prefix + "*" + postfix):
- 746 files.append(f)
- 747
- 748 files = sorted(files)
- 749
- 750 if 'r_start' in kwargs:
- 751 r_start = kwargs.get('r_start')
- 752 if len(r_start) != len(files):
- 753 raise Exception('r_start does not match number of replicas')
- 754 r_start = [o if o else None for o in r_start]
- 755 else:
- 756 r_start = [None] * len(files)
- 757
- 758 if 'r_stop' in kwargs:
- 759 r_stop = kwargs.get('r_stop')
- 760 if len(r_stop) != len(files):
- 761 raise Exception('r_stop does not match number of replicas')
- 762 else:
- 763 r_stop = [None] * len(files)
- 764 rep_names = []
+ 732 Returns
+ 733 -------
+ 734 result : Obs
+ 735 flow observable specified
+ 736 """
+ 737 known_versions = ["openQCD", "sfqcd"]
+ 738
+ 739 if version not in known_versions:
+ 740 raise Exception("Unknown openQCD version.")
+ 741 if "steps" in kwargs:
+ 742 steps = kwargs.get("steps")
+ 743 if version == "sfqcd":
+ 744 if "L" in kwargs:
+ 745 supposed_L = kwargs.get("L")
+ 746 else:
+ 747 supposed_L = None
+ 748 postfix = "gfms"
+ 749 else:
+ 750 if "L" not in kwargs:
+ 751 raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
+ 752 else:
+ 753 L = kwargs.get("L")
+ 754 postfix = "ms"
+ 755
+ 756 if "postfix" in kwargs:
+ 757 postfix = kwargs.get("postfix")
+ 758
+ 759 if "files" in kwargs:
+ 760 known_files = kwargs.get("files")
+ 761 else:
+ 762 known_files = []
+ 763
+ 764 files = _find_files(path, prefix, postfix, "dat", known_files=known_files)
765
- 766 zeuthen = kwargs.get('Zeuthen_flow', False)
- 767 if zeuthen and version not in ['sfqcd']:
- 768 raise Exception('Zeuthen flow can only be used for version==sfqcd')
- 769
- 770 r_start_index = []
- 771 r_stop_index = []
- 772 deltas = []
- 773 configlist = []
- 774 if not zeuthen:
- 775 obspos += 8
- 776 for rep, file in enumerate(files):
- 777 with open(path + "/" + file, "rb") as fp:
- 778
- 779 Q = []
- 780 traj_list = []
- 781 if version in ['sfqcd']:
- 782 t = fp.read(12)
- 783 header = struct.unpack('<iii', t)
- 784 zthfl = header[0] # Zeuthen flow -> if it's equal to 2 it means that the Zeuthen flow is also 'measured' (apart from the Wilson flow)
- 785 ncs = header[1] # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's
- 786 tmax = header[2] # lattice T/a
- 787
- 788 t = fp.read(12)
- 789 Ls = struct.unpack('<iii', t)
- 790 if (Ls[0] == Ls[1] and Ls[1] == Ls[2]):
- 791 L = Ls[0]
- 792 if not (supposed_L == L) and supposed_L:
- 793 raise Exception("It seems the length given in the header and by you contradict each other")
- 794 else:
- 795 raise Exception("Found more than one spatial length in header!")
- 796
- 797 t = fp.read(16)
- 798 header2 = struct.unpack('<dd', t)
- 799 tol = header2[0]
- 800 cmax = header2[1] # highest value of c used
- 801
- 802 if c > cmax:
- 803 raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol))
- 804
- 805 if (zthfl == 2):
- 806 nfl = 2 # number of flows
- 807 else:
- 808 nfl = 1
- 809 iobs = 8 * nfl # number of flow observables calculated
- 810
- 811 while True:
- 812 t = fp.read(4)
- 813 if (len(t) < 4):
- 814 break
- 815 traj_list.append(struct.unpack('i', t)[0]) # trajectory number when measurement was done
- 816
- 817 for j in range(ncs + 1):
- 818 for i in range(iobs):
- 819 t = fp.read(8 * tmax)
- 820 if (i == obspos): # determines the flow observable -> i=0 <-> Zeuthen flow
- 821 Q.append(struct.unpack('d' * tmax, t))
- 822
- 823 else:
- 824 t = fp.read(12)
- 825 header = struct.unpack('<iii', t)
- 826 # step size in integration steps "dnms"
- 827 dn = header[0]
- 828 # number of measurements, so "ntot"/dn
- 829 nn = header[1]
- 830 # lattice T/a
- 831 tmax = header[2]
+ 766 if 'r_start' in kwargs:
+ 767 r_start = kwargs.get('r_start')
+ 768 if len(r_start) != len(files):
+ 769 raise Exception('r_start does not match number of replicas')
+ 770 r_start = [o if o else None for o in r_start]
+ 771 else:
+ 772 r_start = [None] * len(files)
+ 773
+ 774 if 'r_stop' in kwargs:
+ 775 r_stop = kwargs.get('r_stop')
+ 776 if len(r_stop) != len(files):
+ 777 raise Exception('r_stop does not match number of replicas')
+ 778 else:
+ 779 r_stop = [None] * len(files)
+ 780 rep_names = []
+ 781
+ 782 zeuthen = kwargs.get('Zeuthen_flow', False)
+ 783 if zeuthen and version not in ['sfqcd']:
+ 784 raise Exception('Zeuthen flow can only be used for version==sfqcd')
+ 785
+ 786 r_start_index = []
+ 787 r_stop_index = []
+ 788 deltas = []
+ 789 configlist = []
+ 790 if not zeuthen:
+ 791 obspos += 8
+ 792 for rep, file in enumerate(files):
+ 793 with open(path + "/" + file, "rb") as fp:
+ 794
+ 795 Q = []
+ 796 traj_list = []
+ 797 if version in ['sfqcd']:
+ 798 t = fp.read(12)
+ 799 header = struct.unpack('<iii', t)
+ 800 zthfl = header[0] # Zeuthen flow -> if it's equal to 2 it means that the Zeuthen flow is also 'measured' (apart from the Wilson flow)
+ 801 ncs = header[1] # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's
+ 802 tmax = header[2] # lattice T/a
+ 803
+ 804 t = fp.read(12)
+ 805 Ls = struct.unpack('<iii', t)
+ 806 if (Ls[0] == Ls[1] and Ls[1] == Ls[2]):
+ 807 L = Ls[0]
+ 808 if not (supposed_L == L) and supposed_L:
+ 809 raise Exception("It seems the length given in the header and by you contradict each other")
+ 810 else:
+ 811 raise Exception("Found more than one spatial length in header!")
+ 812
+ 813 t = fp.read(16)
+ 814 header2 = struct.unpack('<dd', t)
+ 815 tol = header2[0]
+ 816 cmax = header2[1] # highest value of c used
+ 817
+ 818 if c > cmax:
+ 819 raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol))
+ 820
+ 821 if (zthfl == 2):
+ 822 nfl = 2 # number of flows
+ 823 else:
+ 824 nfl = 1
+ 825 iobs = 8 * nfl # number of flow observables calculated
+ 826
+ 827 while True:
+ 828 t = fp.read(4)
+ 829 if (len(t) < 4):
+ 830 break
+ 831 traj_list.append(struct.unpack('i', t)[0]) # trajectory number when measurement was done
832
- 833 t = fp.read(8)
- 834 eps = struct.unpack('d', t)[0]
- 835
- 836 while True:
- 837 t = fp.read(4)
- 838 if (len(t) < 4):
- 839 break
- 840 traj_list.append(struct.unpack('i', t)[0])
- 841 # Wsl
- 842 t = fp.read(8 * tmax * (nn + 1))
- 843 # Ysl
- 844 t = fp.read(8 * tmax * (nn + 1))
- 845 # Qsl, which is asked for in this method
- 846 t = fp.read(8 * tmax * (nn + 1))
- 847 # unpack the array of Qtops,
- 848 # on each timeslice t=0,...,tmax-1 and the
- 849 # measurement number in = 0...nn (see README.qcd1)
- 850 tmpd = struct.unpack('d' * tmax * (nn + 1), t)
- 851 Q.append(tmpd)
- 852
- 853 if len(np.unique(np.diff(traj_list))) != 1:
- 854 raise Exception("Irregularities in stepsize found")
- 855 else:
- 856 if 'steps' in kwargs:
- 857 if steps != traj_list[1] - traj_list[0]:
- 858 raise Exception("steps and the found stepsize are not the same")
- 859 else:
- 860 steps = traj_list[1] - traj_list[0]
- 861
- 862 configlist.append([tr // steps // dtr_cnfg for tr in traj_list])
- 863 if configlist[-1][0] > 1:
- 864 offset = configlist[-1][0] - 1
- 865 warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % (
- 866 offset, offset * steps))
- 867 configlist[-1] = [item - offset for item in configlist[-1]]
+ 833 for j in range(ncs + 1):
+ 834 for i in range(iobs):
+ 835 t = fp.read(8 * tmax)
+ 836 if (i == obspos): # determines the flow observable -> i=0 <-> Zeuthen flow
+ 837 Q.append(struct.unpack('d' * tmax, t))
+ 838
+ 839 else:
+ 840 t = fp.read(12)
+ 841 header = struct.unpack('<iii', t)
+ 842 # step size in integration steps "dnms"
+ 843 dn = header[0]
+ 844 # number of measurements, so "ntot"/dn
+ 845 nn = header[1]
+ 846 # lattice T/a
+ 847 tmax = header[2]
+ 848
+ 849 t = fp.read(8)
+ 850 eps = struct.unpack('d', t)[0]
+ 851
+ 852 while True:
+ 853 t = fp.read(4)
+ 854 if (len(t) < 4):
+ 855 break
+ 856 traj_list.append(struct.unpack('i', t)[0])
+ 857 # Wsl
+ 858 t = fp.read(8 * tmax * (nn + 1))
+ 859 # Ysl
+ 860 t = fp.read(8 * tmax * (nn + 1))
+ 861 # Qsl, which is asked for in this method
+ 862 t = fp.read(8 * tmax * (nn + 1))
+ 863 # unpack the array of Qtops,
+ 864 # on each timeslice t=0,...,tmax-1 and the
+ 865 # measurement number in = 0...nn (see README.qcd1)
+ 866 tmpd = struct.unpack('d' * tmax * (nn + 1), t)
+ 867 Q.append(tmpd)
868
- 869 if r_start[rep] is None:
- 870 r_start_index.append(0)
+ 869 if len(np.unique(np.diff(traj_list))) != 1:
+ 870 raise Exception("Irregularities in stepsize found")
871 else:
- 872 try:
- 873 r_start_index.append(configlist[-1].index(r_start[rep]))
- 874 except ValueError:
- 875 raise Exception('Config %d not in file with range [%d, %d]' % (
- 876 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
+ 872 if 'steps' in kwargs:
+ 873 if steps != traj_list[1] - traj_list[0]:
+ 874 raise Exception("steps and the found stepsize are not the same")
+ 875 else:
+ 876 steps = traj_list[1] - traj_list[0]
877
- 878 if r_stop[rep] is None:
- 879 r_stop_index.append(len(configlist[-1]) - 1)
- 880 else:
- 881 try:
- 882 r_stop_index.append(configlist[-1].index(r_stop[rep]))
- 883 except ValueError:
- 884 raise Exception('Config %d not in file with range [%d, %d]' % (
- 885 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
- 886
- 887 if version in ['sfqcd']:
- 888 cstepsize = cmax / ncs
- 889 index_aim = round(c / cstepsize)
- 890 else:
- 891 t_aim = (c * L) ** 2 / 8
- 892 index_aim = round(t_aim / eps / dn)
+ 878 configlist.append([tr // steps // dtr_cnfg for tr in traj_list])
+ 879 if configlist[-1][0] > 1:
+ 880 offset = configlist[-1][0] - 1
+ 881 warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % (
+ 882 offset, offset * steps))
+ 883 configlist[-1] = [item - offset for item in configlist[-1]]
+ 884
+ 885 if r_start[rep] is None:
+ 886 r_start_index.append(0)
+ 887 else:
+ 888 try:
+ 889 r_start_index.append(configlist[-1].index(r_start[rep]))
+ 890 except ValueError:
+ 891 raise Exception('Config %d not in file with range [%d, %d]' % (
+ 892 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
893
- 894 Q_sum = []
- 895 for i, item in enumerate(Q):
- 896 if sum_t is True:
- 897 Q_sum.append([sum(item[current:current + tmax])
- 898 for current in range(0, len(item), tmax)])
- 899 else:
- 900 Q_sum.append([item[int(tmax / 2)]])
- 901 Q_top = []
- 902 if version in ['sfqcd']:
- 903 for i in range(len(Q_sum) // (ncs + 1)):
- 904 Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0])
- 905 else:
- 906 for i in range(len(Q) // dtr_cnfg):
- 907 Q_top.append(Q_sum[dtr_cnfg * i][index_aim])
- 908 if len(Q_top) != len(traj_list) // dtr_cnfg:
- 909 raise Exception("qtops and traj_list dont have the same length")
- 910
- 911 if kwargs.get('integer_charge', False):
- 912 Q_top = [round(q) for q in Q_top]
- 913
- 914 truncated_file = file[:-len(postfix)]
- 915
- 916 if "names" not in kwargs:
- 917 try:
- 918 idx = truncated_file.index('r')
- 919 except Exception:
- 920 if "names" not in kwargs:
- 921 raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
- 922 ens_name = truncated_file[:idx]
- 923 rep_names.append(ens_name + '|' + truncated_file[idx:])
- 924 else:
- 925 names = kwargs.get("names")
- 926 rep_names = names
- 927 deltas.append(Q_top)
- 928
- 929 idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))]
- 930 deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))]
- 931 result = Obs(deltas, rep_names, idl=idl)
- 932 result.tag = {"T": tmax - 1,
- 933 "L": L}
- 934 return result
- 935
- 936
- 937def qtop_projection(qtop, target=0):
- 938 """Returns the projection to the topological charge sector defined by target.
- 939
- 940 Parameters
- 941 ----------
- 942 path : Obs
- 943 Topological charge.
- 944 target : int
- 945 Specifies the topological sector to be reweighted to (default 0)
- 946
- 947 Returns
- 948 -------
- 949 reto : Obs
- 950 projection to the topological charge sector defined by target
- 951 """
- 952 if qtop.reweighted:
- 953 raise Exception('You can not use a reweighted observable for reweighting!')
- 954
- 955 proj_qtop = []
- 956 for n in qtop.deltas:
- 957 proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
- 958
- 959 reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
- 960 return reto
- 961
+ 894 if r_stop[rep] is None:
+ 895 r_stop_index.append(len(configlist[-1]) - 1)
+ 896 else:
+ 897 try:
+ 898 r_stop_index.append(configlist[-1].index(r_stop[rep]))
+ 899 except ValueError:
+ 900 raise Exception('Config %d not in file with range [%d, %d]' % (
+ 901 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
+ 902
+ 903 if version in ['sfqcd']:
+ 904 cstepsize = cmax / ncs
+ 905 index_aim = round(c / cstepsize)
+ 906 else:
+ 907 t_aim = (c * L) ** 2 / 8
+ 908 index_aim = round(t_aim / eps / dn)
+ 909
+ 910 Q_sum = []
+ 911 for i, item in enumerate(Q):
+ 912 if sum_t is True:
+ 913 Q_sum.append([sum(item[current:current + tmax])
+ 914 for current in range(0, len(item), tmax)])
+ 915 else:
+ 916 Q_sum.append([item[int(tmax / 2)]])
+ 917 Q_top = []
+ 918 if version in ['sfqcd']:
+ 919 for i in range(len(Q_sum) // (ncs + 1)):
+ 920 Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0])
+ 921 else:
+ 922 for i in range(len(Q) // dtr_cnfg):
+ 923 Q_top.append(Q_sum[dtr_cnfg * i][index_aim])
+ 924 if len(Q_top) != len(traj_list) // dtr_cnfg:
+ 925 raise Exception("qtops and traj_list dont have the same length")
+ 926
+ 927 if kwargs.get('integer_charge', False):
+ 928 Q_top = [round(q) for q in Q_top]
+ 929
+ 930 truncated_file = file[:-len(postfix)]
+ 931
+ 932 if "names" not in kwargs:
+ 933 try:
+ 934 idx = truncated_file.index('r')
+ 935 except Exception:
+ 936 if "names" not in kwargs:
+ 937 raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
+ 938 ens_name = truncated_file[:idx]
+ 939 rep_names.append(ens_name + '|' + truncated_file[idx:])
+ 940 else:
+ 941 names = kwargs.get("names")
+ 942 rep_names = names
+ 943 deltas.append(Q_top)
+ 944
+ 945 idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))]
+ 946 deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))]
+ 947 result = Obs(deltas, rep_names, idl=idl)
+ 948 result.tag = {"T": tmax - 1,
+ 949 "L": L}
+ 950 return result
+ 951
+ 952
+ 953def qtop_projection(qtop, target=0):
+ 954 """Returns the projection to the topological charge sector defined by target.
+ 955
+ 956 Parameters
+ 957 ----------
+ 958 path : Obs
+ 959 Topological charge.
+ 960 target : int
+ 961 Specifies the topological sector to be reweighted to (default 0)
962
- 963def read_qtop_sector(path, prefix, c, target=0, **kwargs):
- 964 """Constructs reweighting factors to a specified topological sector.
- 965
- 966 Parameters
- 967 ----------
- 968 path : str
- 969 path of the measurement files
- 970 prefix : str
- 971 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
- 972 c : double
- 973 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
- 974 target : int
- 975 Specifies the topological sector to be reweighted to (default 0)
- 976 dtr_cnfg : int
- 977 (optional) parameter that specifies the number of trajectories
- 978 between two configs.
- 979 if it is not set, the distance between two measurements
- 980 in the file is assumed to be the distance between two configurations.
- 981 steps : int
- 982 (optional) Distance between two configurations in units of trajectories /
- 983 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- 984 version : str
- 985 version string of the openQCD (sfqcd) version used to create
- 986 the ensemble. Default is 2.0. May also be set to sfqcd.
- 987 L : int
- 988 spatial length of the lattice in L/a.
- 989 HAS to be set if version != sfqcd, since openQCD does not provide
- 990 this in the header
- 991 r_start : list
- 992 offset of the first ensemble, making it easier to match
- 993 later on with other Obs
- 994 r_stop : list
- 995 last configurations that need to be read (per replicum)
- 996 files : list
- 997 specify the exact files that need to be read
- 998 from path, practical if e.g. only one replicum is needed
- 999 names : list
-1000 Alternative labeling for replicas/ensembles.
-1001 Has to have the appropriate length
-1002 Zeuthen_flow : bool
-1003 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
-1004 for version=='sfqcd' If False, the Wilson flow is used.
-1005
-1006 Returns
-1007 -------
-1008 reto : Obs
-1009 projection to the topological charge sector defined by target
-1010 """
-1011
-1012 if not isinstance(target, int):
-1013 raise Exception("'target' has to be an integer.")
-1014
-1015 kwargs['integer_charge'] = True
-1016 qtop = read_qtop(path, prefix, c, **kwargs)
-1017
-1018 return qtop_projection(qtop, target=target)
-1019
-1020
-1021def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
-1022 """
-1023 Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
-1024
-1025 Parameters
-1026 ----------
-1027 path : str
-1028 The directory to search for the files in.
-1029 prefix : str
-1030 The prefix to match the files against.
-1031 qc : str
-1032 The quark combination extension to match the files against.
-1033 corr : str
-1034 The correlator to extract data for.
-1035 sep : str, optional
-1036 The separator to use when parsing the replika names.
-1037 **kwargs
-1038 Additional keyword arguments. The following keyword arguments are recognized:
-1039
-1040 - names (List[str]): A list of names to use for the replicas.
-1041
-1042 Returns
-1043 -------
-1044 Corr
-1045 A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
-1046 or
-1047 CObs
-1048 A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
-1049
-1050
-1051 Raises
-1052 ------
-1053 FileNotFoundError
-1054 If no files matching the specified prefix and quark combination extension are found in the specified directory.
-1055 IOError
-1056 If there is an error reading a file.
-1057 struct.error
-1058 If there is an error unpacking binary data.
-1059 """
-1060
-1061 found = []
-1062 files = []
-1063 names = []
-1064
-1065 if "names" in kwargs:
-1066 names = kwargs.get("names")
-1067
-1068 for (dirpath, dirnames, filenames) in os.walk(path + "/"):
-1069 found.extend(filenames)
-1070 break
-1071
-1072 for f in found:
-1073 if fnmatch.fnmatch(f, prefix + "*.ms5_xsf_" + qc + ".dat"):
-1074 files.append(f)
-1075 if "names" not in kwargs:
-1076 if not sep == "":
-1077 se = f.split(".")[0]
-1078 for s in f.split(".")[1:-1]:
-1079 se += "." + s
-1080 names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
-1081 else:
-1082 names.append(prefix)
-1083
-1084 names = sorted(names)
-1085 files = sorted(files)
-1086
-1087 cnfgs = []
-1088 realsamples = []
-1089 imagsamples = []
-1090 repnum = 0
-1091 for file in files:
-1092 with open(path + "/" + file, "rb") as fp:
+ 963 Returns
+ 964 -------
+ 965 reto : Obs
+ 966 projection to the topological charge sector defined by target
+ 967 """
+ 968 if qtop.reweighted:
+ 969 raise Exception('You can not use a reweighted observable for reweighting!')
+ 970
+ 971 proj_qtop = []
+ 972 for n in qtop.deltas:
+ 973 proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
+ 974
+ 975 reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
+ 976 return reto
+ 977
+ 978
+ 979def read_qtop_sector(path, prefix, c, target=0, **kwargs):
+ 980 """Constructs reweighting factors to a specified topological sector.
+ 981
+ 982 Parameters
+ 983 ----------
+ 984 path : str
+ 985 path of the measurement files
+ 986 prefix : str
+ 987 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
+ 988 c : double
+ 989 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
+ 990 target : int
+ 991 Specifies the topological sector to be reweighted to (default 0)
+ 992 dtr_cnfg : int
+ 993 (optional) parameter that specifies the number of trajectories
+ 994 between two configs.
+ 995 if it is not set, the distance between two measurements
+ 996 in the file is assumed to be the distance between two configurations.
+ 997 steps : int
+ 998 (optional) Distance between two configurations in units of trajectories /
+ 999 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+1000 version : str
+1001 version string of the openQCD (sfqcd) version used to create
+1002 the ensemble. Default is 2.0. May also be set to sfqcd.
+1003 L : int
+1004 spatial length of the lattice in L/a.
+1005 HAS to be set if version != sfqcd, since openQCD does not provide
+1006 this in the header
+1007 r_start : list
+1008 offset of the first ensemble, making it easier to match
+1009 later on with other Obs
+1010 r_stop : list
+1011 last configurations that need to be read (per replicum)
+1012 files : list
+1013 specify the exact files that need to be read
+1014 from path, practical if e.g. only one replicum is needed
+1015 names : list
+1016 Alternative labeling for replicas/ensembles.
+1017 Has to have the appropriate length
+1018 Zeuthen_flow : bool
+1019 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+1020 for version=='sfqcd' If False, the Wilson flow is used.
+1021
+1022 Returns
+1023 -------
+1024 reto : Obs
+1025 projection to the topological charge sector defined by target
+1026 """
+1027
+1028 if not isinstance(target, int):
+1029 raise Exception("'target' has to be an integer.")
+1030
+1031 kwargs['integer_charge'] = True
+1032 qtop = read_qtop(path, prefix, c, **kwargs)
+1033
+1034 return qtop_projection(qtop, target=target)
+1035
+1036
+1037def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
+1038 """
+1039 Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
+1040
+1041 Parameters
+1042 ----------
+1043 path : str
+1044 The directory to search for the files in.
+1045 prefix : str
+1046 The prefix to match the files against.
+1047 qc : str
+1048 The quark combination extension to match the files against.
+1049 corr : str
+1050 The correlator to extract data for.
+1051 sep : str, optional
+1052 The separator to use when parsing the replika names.
+1053 **kwargs
+1054 Additional keyword arguments. The following keyword arguments are recognized:
+1055
+1056 - names (List[str]): A list of names to use for the replicas.
+1057
+1058 Returns
+1059 -------
+1060 Corr
+1061 A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
+1062 or
+1063 CObs
+1064 A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
+1065
+1066
+1067 Raises
+1068 ------
+1069 FileNotFoundError
+1070 If no files matching the specified prefix and quark combination extension are found in the specified directory.
+1071 IOError
+1072 If there is an error reading a file.
+1073 struct.error
+1074 If there is an error unpacking binary data.
+1075 """
+1076
+1077 # found = []
+1078 files = []
+1079 names = []
+1080
+1081 # test if the input is correct
+1082 if qc not in ['dd', 'ud', 'du', 'uu']:
+1083 raise Exception("Unknown quark conbination!")
+1084
+1085 if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
+1086 raise Exception("Unknown correlator!")
+1087
+1088 if "files" in kwargs:
+1089 known_files = kwargs.get("files")
+1090 else:
+1091 known_files = []
+1092 files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
1093
-1094 t = fp.read(8)
-1095 kappa = struct.unpack('d', t)[0]
-1096 t = fp.read(8)
-1097 csw = struct.unpack('d', t)[0]
-1098 t = fp.read(8)
-1099 dF = struct.unpack('d', t)[0]
-1100 t = fp.read(8)
-1101 zF = struct.unpack('d', t)[0]
-1102
-1103 t = fp.read(4)
-1104 tmax = struct.unpack('i', t)[0]
-1105 t = fp.read(4)
-1106 bnd = struct.unpack('i', t)[0]
-1107
-1108 placesBI = ["gS", "gP",
-1109 "gA", "gV",
-1110 "gVt", "lA",
-1111 "lV", "lVt",
-1112 "lT", "lTt"]
-1113 placesBB = ["g1", "l1"]
-1114
-1115 # the chunks have the following structure:
-1116 # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
-1117
-1118 chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
-1119 packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
-1120 cnfgs.append([])
-1121 realsamples.append([])
-1122 imagsamples.append([])
-1123 for t in range(tmax):
-1124 realsamples[repnum].append([])
-1125 imagsamples[repnum].append([])
-1126
-1127 while True:
-1128 cnfgt = fp.read(chunksize)
-1129 if not cnfgt:
-1130 break
-1131 asascii = struct.unpack(packstr, cnfgt)
-1132 cnfg = asascii[0]
-1133 cnfgs[repnum].append(cnfg)
-1134
-1135 if corr not in placesBB:
-1136 tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
-1137 else:
-1138 tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
+1094 if "names" in kwargs:
+1095 names = kwargs.get("names")
+1096 else:
+1097 for f in files:
+1098 if not sep == "":
+1099 se = f.split(".")[0]
+1100 for s in f.split(".")[1:-2]:
+1101 se += "." + s
+1102 names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
+1103 else:
+1104 names.append(prefix)
+1105
+1106 names = sorted(names)
+1107 files = sorted(files)
+1108
+1109 cnfgs = []
+1110 realsamples = []
+1111 imagsamples = []
+1112 repnum = 0
+1113 for file in files:
+1114 with open(path + "/" + file, "rb") as fp:
+1115
+1116 t = fp.read(8)
+1117 kappa = struct.unpack('d', t)[0]
+1118 t = fp.read(8)
+1119 csw = struct.unpack('d', t)[0]
+1120 t = fp.read(8)
+1121 dF = struct.unpack('d', t)[0]
+1122 t = fp.read(8)
+1123 zF = struct.unpack('d', t)[0]
+1124
+1125 t = fp.read(4)
+1126 tmax = struct.unpack('i', t)[0]
+1127 t = fp.read(4)
+1128 bnd = struct.unpack('i', t)[0]
+1129
+1130 placesBI = ["gS", "gP",
+1131 "gA", "gV",
+1132 "gVt", "lA",
+1133 "lV", "lVt",
+1134 "lT", "lTt"]
+1135 placesBB = ["g1", "l1"]
+1136
+1137 # the chunks have the following structure:
+1138 # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
1139
-1140 corrres = [[], []]
-1141 for i in range(len(tmpcorr)):
-1142 corrres[i % 2].append(tmpcorr[i])
-1143 for t in range(int(len(tmpcorr) / 2)):
-1144 realsamples[repnum][t].append(corrres[0][t])
-1145 for t in range(int(len(tmpcorr) / 2)):
-1146 imagsamples[repnum][t].append(corrres[1][t])
-1147 repnum += 1
+1140 chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
+1141 packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
+1142 cnfgs.append([])
+1143 realsamples.append([])
+1144 imagsamples.append([])
+1145 for t in range(tmax):
+1146 realsamples[repnum].append([])
+1147 imagsamples[repnum].append([])
1148
-1149 s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
-1150 for rep in range(1, repnum):
-1151 s += ", " + str(len(realsamples[rep][t]))
-1152 s += " samples"
-1153 print(s)
-1154 print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
-1155
-1156 # we have the data now... but we need to re format the whole thing and put it into Corr objects.
-1157
-1158 compObs = []
-1159
-1160 for t in range(int(len(tmpcorr) / 2)):
-1161 compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
-1162 Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
-1163
-1164 if len(compObs) == 1:
-1165 return compObs[0]
-1166 else:
-1167 return Corr(compObs)
+1149 while True:
+1150 cnfgt = fp.read(chunksize)
+1151 if not cnfgt:
+1152 break
+1153 asascii = struct.unpack(packstr, cnfgt)
+1154 cnfg = asascii[0]
+1155 cnfgs[repnum].append(cnfg)
+1156
+1157 if corr not in placesBB:
+1158 tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
+1159 else:
+1160 tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
+1161
+1162 corrres = [[], []]
+1163 for i in range(len(tmpcorr)):
+1164 corrres[i % 2].append(tmpcorr[i])
+1165 for t in range(int(len(tmpcorr) / 2)):
+1166 realsamples[repnum][t].append(corrres[0][t])
+1167 for t in range(int(len(tmpcorr) / 2)):
+1168 imagsamples[repnum][t].append(corrres[1][t])
+1169 repnum += 1
+1170
+1171 s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
+1172 for rep in range(1, repnum):
+1173 s += ", " + str(len(realsamples[rep][t]))
+1174 s += " samples"
+1175 print(s)
+1176 print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
+1177
+1178 # we have the data now... but we need to re format the whole thing and put it into Corr objects.
+1179
+1180 compObs = []
+1181
+1182 for t in range(int(len(tmpcorr) / 2)):
+1183 compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
+1184 Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
+1185
+1186 if len(compObs) == 1:
+1187 return compObs[0]
+1188 else:
+1189 return Corr(compObs)
@@ -1276,228 +1298,221 @@
- 16def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
- 17 """Read rwms format from given folder structure. Returns a list of length nrw
- 18
- 19 Parameters
- 20 ----------
- 21 path : str
- 22 path that contains the data files
- 23 prefix : str
- 24 all files in path that start with prefix are considered as input files.
- 25 May be used together postfix to consider only special file endings.
- 26 Prefix is ignored, if the keyword 'files' is used.
- 27 version : str
- 28 version of openQCD, default 2.0
- 29 names : list
- 30 list of names that is assigned to the data according according
- 31 to the order in the file list. Use careful, if you do not provide file names!
- 32 r_start : list
- 33 list which contains the first config to be read for each replicum
- 34 r_stop : list
- 35 list which contains the last config to be read for each replicum
- 36 r_step : int
- 37 integer that defines a fixed step size between two measurements (in units of configs)
- 38 If not given, r_step=1 is assumed.
- 39 postfix : str
- 40 postfix of the file to read, e.g. '.ms1' for openQCD-files
- 41 files : list
- 42 list which contains the filenames to be read. No automatic detection of
- 43 files performed if given.
- 44 print_err : bool
- 45 Print additional information that is useful for debugging.
- 46
- 47 Returns
- 48 -------
- 49 rwms : Obs
- 50 Reweighting factors read
- 51 """
- 52 known_oqcd_versions = ['1.4', '1.6', '2.0']
- 53 if not (version in known_oqcd_versions):
- 54 raise Exception('Unknown openQCD version defined!')
- 55 print("Working with openQCD version " + version)
- 56 if 'postfix' in kwargs:
- 57 postfix = kwargs.get('postfix')
- 58 else:
- 59 postfix = ''
- 60 ls = []
- 61 for (dirpath, dirnames, filenames) in os.walk(path):
- 62 ls.extend(filenames)
- 63 break
- 64
- 65 if not ls:
- 66 raise Exception(f"Error, directory '{path}' not found")
- 67 if 'files' in kwargs:
- 68 ls = kwargs.get('files')
- 69 else:
- 70 for exc in ls:
- 71 if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'):
- 72 ls = list(set(ls) - set([exc]))
- 73 if len(ls) > 1:
- 74 ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
- 75 replica = len(ls)
- 76
- 77 if 'r_start' in kwargs:
- 78 r_start = kwargs.get('r_start')
- 79 if len(r_start) != replica:
- 80 raise Exception('r_start does not match number of replicas')
- 81 r_start = [o if o else None for o in r_start]
- 82 else:
- 83 r_start = [None] * replica
- 84
- 85 if 'r_stop' in kwargs:
- 86 r_stop = kwargs.get('r_stop')
- 87 if len(r_stop) != replica:
- 88 raise Exception('r_stop does not match number of replicas')
- 89 else:
- 90 r_stop = [None] * replica
- 91
- 92 if 'r_step' in kwargs:
- 93 r_step = kwargs.get('r_step')
- 94 else:
- 95 r_step = 1
- 96
- 97 print('Read reweighting factors from', prefix[:-1], ',',
- 98 replica, 'replica', end='')
- 99
-100 if names is None:
-101 rep_names = []
-102 for entry in ls:
-103 truncated_entry = entry
-104 suffixes = [".dat", ".rwms", ".ms1"]
-105 for suffix in suffixes:
-106 if truncated_entry.endswith(suffix):
-107 truncated_entry = truncated_entry[0:-len(suffix)]
-108 idx = truncated_entry.index('r')
-109 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
-110 else:
-111 rep_names = names
-112
-113 print_err = 0
-114 if 'print_err' in kwargs:
-115 print_err = 1
-116 print()
+ 56def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
+ 57 """Read rwms format from given folder structure. Returns a list of length nrw
+ 58
+ 59 Parameters
+ 60 ----------
+ 61 path : str
+ 62 path that contains the data files
+ 63 prefix : str
+ 64 all files in path that start with prefix are considered as input files.
+ 65 May be used together postfix to consider only special file endings.
+ 66 Prefix is ignored, if the keyword 'files' is used.
+ 67 version : str
+ 68 version of openQCD, default 2.0
+ 69 names : list
+ 70 list of names that is assigned to the data according according
+ 71 to the order in the file list. Use careful, if you do not provide file names!
+ 72 r_start : list
+ 73 list which contains the first config to be read for each replicum
+ 74 r_stop : list
+ 75 list which contains the last config to be read for each replicum
+ 76 r_step : int
+ 77 integer that defines a fixed step size between two measurements (in units of configs)
+ 78 If not given, r_step=1 is assumed.
+ 79 postfix : str
+ 80 postfix of the file to read, e.g. '.ms1' for openQCD-files
+ 81 files : list
+ 82 list which contains the filenames to be read. No automatic detection of
+ 83 files performed if given.
+ 84 print_err : bool
+ 85 Print additional information that is useful for debugging.
+ 86
+ 87 Returns
+ 88 -------
+ 89 rwms : Obs
+ 90 Reweighting factors read
+ 91 """
+ 92 known_oqcd_versions = ['1.4', '1.6', '2.0']
+ 93 if not (version in known_oqcd_versions):
+ 94 raise Exception('Unknown openQCD version defined!')
+ 95 print("Working with openQCD version " + version)
+ 96 if 'postfix' in kwargs:
+ 97 postfix = kwargs.get('postfix')
+ 98 else:
+ 99 postfix = ''
+100
+101 if 'files' in kwargs:
+102 known_files = kwargs.get('files')
+103 else:
+104 known_files = []
+105
+106 ls = _find_files(path, prefix, postfix, 'dat', known_files=known_files)
+107
+108 replica = len(ls)
+109
+110 if 'r_start' in kwargs:
+111 r_start = kwargs.get('r_start')
+112 if len(r_start) != replica:
+113 raise Exception('r_start does not match number of replicas')
+114 r_start = [o if o else None for o in r_start]
+115 else:
+116 r_start = [None] * replica
117
-118 deltas = []
-119
-120 configlist = []
-121 r_start_index = []
-122 r_stop_index = []
-123
-124 for rep in range(replica):
-125 tmp_array = []
-126 with open(path + '/' + ls[rep], 'rb') as fp:
-127
-128 t = fp.read(4) # number of reweighting factors
-129 if rep == 0:
-130 nrw = struct.unpack('i', t)[0]
-131 if version == '2.0':
-132 nrw = int(nrw / 2)
-133 for k in range(nrw):
-134 deltas.append([])
-135 else:
-136 if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
-137 raise Exception('Error: different number of reweighting factors for replicum', rep)
-138
-139 for k in range(nrw):
-140 tmp_array.append([])
-141
-142 # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
-143 nfct = []
-144 if version in ['1.6', '2.0']:
-145 for i in range(nrw):
-146 t = fp.read(4)
-147 nfct.append(struct.unpack('i', t)[0])
-148 else:
-149 for i in range(nrw):
-150 nfct.append(1)
-151
-152 nsrc = []
-153 for i in range(nrw):
-154 t = fp.read(4)
-155 nsrc.append(struct.unpack('i', t)[0])
-156 if version == '2.0':
-157 if not struct.unpack('i', fp.read(4))[0] == 0:
-158 print('something is wrong!')
-159
-160 configlist.append([])
-161 while True:
-162 t = fp.read(4)
-163 if len(t) < 4:
-164 break
-165 config_no = struct.unpack('i', t)[0]
-166 configlist[-1].append(config_no)
-167 for i in range(nrw):
-168 if (version == '2.0'):
-169 tmpd = _read_array_openQCD2(fp)
-170 tmpd = _read_array_openQCD2(fp)
-171 tmp_rw = tmpd['arr']
-172 tmp_nfct = 1.0
-173 for j in range(tmpd['n'][0]):
-174 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
-175 if print_err:
-176 print(config_no, i, j,
-177 np.mean(np.exp(-np.asarray(tmp_rw[j]))),
-178 np.std(np.exp(-np.asarray(tmp_rw[j]))))
-179 print('Sources:',
-180 np.exp(-np.asarray(tmp_rw[j])))
-181 print('Partial factor:', tmp_nfct)
-182 elif version == '1.6' or version == '1.4':
-183 tmp_nfct = 1.0
-184 for j in range(nfct[i]):
-185 t = fp.read(8 * nsrc[i])
-186 t = fp.read(8 * nsrc[i])
-187 tmp_rw = struct.unpack('d' * nsrc[i], t)
-188 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
-189 if print_err:
-190 print(config_no, i, j,
-191 np.mean(np.exp(-np.asarray(tmp_rw))),
-192 np.std(np.exp(-np.asarray(tmp_rw))))
-193 print('Sources:', np.exp(-np.asarray(tmp_rw)))
-194 print('Partial factor:', tmp_nfct)
-195 tmp_array[i].append(tmp_nfct)
-196
-197 diffmeas = configlist[-1][-1] - configlist[-1][-2]
-198 configlist[-1] = [item // diffmeas for item in configlist[-1]]
-199 if configlist[-1][0] > 1 and diffmeas > 1:
-200 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
-201 offset = configlist[-1][0] - 1
-202 configlist[-1] = [item - offset for item in configlist[-1]]
-203
-204 if r_start[rep] is None:
-205 r_start_index.append(0)
-206 else:
-207 try:
-208 r_start_index.append(configlist[-1].index(r_start[rep]))
-209 except ValueError:
-210 raise Exception('Config %d not in file with range [%d, %d]' % (
-211 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
-212
-213 if r_stop[rep] is None:
-214 r_stop_index.append(len(configlist[-1]) - 1)
-215 else:
-216 try:
-217 r_stop_index.append(configlist[-1].index(r_stop[rep]))
-218 except ValueError:
-219 raise Exception('Config %d not in file with range [%d, %d]' % (
-220 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
-221
-222 for k in range(nrw):
-223 deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
-224
-225 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
-226 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
-227 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
-228 if np.any([step != 1 for step in stepsizes]):
-229 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
-230
-231 print(',', nrw, 'reweighting factors with', nsrc, 'sources')
-232 result = []
-233 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
-234
-235 for t in range(nrw):
-236 result.append(Obs(deltas[t], rep_names, idl=idl))
-237 return result
+118 if 'r_stop' in kwargs:
+119 r_stop = kwargs.get('r_stop')
+120 if len(r_stop) != replica:
+121 raise Exception('r_stop does not match number of replicas')
+122 else:
+123 r_stop = [None] * replica
+124
+125 if 'r_step' in kwargs:
+126 r_step = kwargs.get('r_step')
+127 else:
+128 r_step = 1
+129
+130 print('Read reweighting factors from', prefix[:-1], ',',
+131 replica, 'replica', end='')
+132
+133 if names is None:
+134 rep_names = []
+135 for entry in ls:
+136 truncated_entry = entry
+137 suffixes = [".dat", ".rwms", ".ms1"]
+138 for suffix in suffixes:
+139 if truncated_entry.endswith(suffix):
+140 truncated_entry = truncated_entry[0:-len(suffix)]
+141 idx = truncated_entry.index('r')
+142 rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
+143 else:
+144 rep_names = names
+145
+146 print_err = 0
+147 if 'print_err' in kwargs:
+148 print_err = 1
+149 print()
+150
+151 deltas = []
+152
+153 configlist = []
+154 r_start_index = []
+155 r_stop_index = []
+156
+157 for rep in range(replica):
+158 tmp_array = []
+159 with open(path + '/' + ls[rep], 'rb') as fp:
+160
+161 t = fp.read(4) # number of reweighting factors
+162 if rep == 0:
+163 nrw = struct.unpack('i', t)[0]
+164 if version == '2.0':
+165 nrw = int(nrw / 2)
+166 for k in range(nrw):
+167 deltas.append([])
+168 else:
+169 if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
+170 raise Exception('Error: different number of reweighting factors for replicum', rep)
+171
+172 for k in range(nrw):
+173 tmp_array.append([])
+174
+175 # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
+176 nfct = []
+177 if version in ['1.6', '2.0']:
+178 for i in range(nrw):
+179 t = fp.read(4)
+180 nfct.append(struct.unpack('i', t)[0])
+181 else:
+182 for i in range(nrw):
+183 nfct.append(1)
+184
+185 nsrc = []
+186 for i in range(nrw):
+187 t = fp.read(4)
+188 nsrc.append(struct.unpack('i', t)[0])
+189 if version == '2.0':
+190 if not struct.unpack('i', fp.read(4))[0] == 0:
+191 raise Exception("You are using the input for openQCD version 2.0, this is not correct.")
+192
+193 configlist.append([])
+194 while True:
+195 t = fp.read(4)
+196 if len(t) < 4:
+197 break
+198 config_no = struct.unpack('i', t)[0]
+199 configlist[-1].append(config_no)
+200 for i in range(nrw):
+201 if (version == '2.0'):
+202 tmpd = _read_array_openQCD2(fp)
+203 tmpd = _read_array_openQCD2(fp)
+204 tmp_rw = tmpd['arr']
+205 tmp_nfct = 1.0
+206 for j in range(tmpd['n'][0]):
+207 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
+208 if print_err:
+209 print(config_no, i, j,
+210 np.mean(np.exp(-np.asarray(tmp_rw[j]))),
+211 np.std(np.exp(-np.asarray(tmp_rw[j]))))
+212 print('Sources:',
+213 np.exp(-np.asarray(tmp_rw[j])))
+214 print('Partial factor:', tmp_nfct)
+215 elif version == '1.6' or version == '1.4':
+216 tmp_nfct = 1.0
+217 for j in range(nfct[i]):
+218 t = fp.read(8 * nsrc[i])
+219 t = fp.read(8 * nsrc[i])
+220 tmp_rw = struct.unpack('d' * nsrc[i], t)
+221 tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
+222 if print_err:
+223 print(config_no, i, j,
+224 np.mean(np.exp(-np.asarray(tmp_rw))),
+225 np.std(np.exp(-np.asarray(tmp_rw))))
+226 print('Sources:', np.exp(-np.asarray(tmp_rw)))
+227 print('Partial factor:', tmp_nfct)
+228 tmp_array[i].append(tmp_nfct)
+229
+230 diffmeas = configlist[-1][-1] - configlist[-1][-2]
+231 configlist[-1] = [item // diffmeas for item in configlist[-1]]
+232 if configlist[-1][0] > 1 and diffmeas > 1:
+233 warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
+234 offset = configlist[-1][0] - 1
+235 configlist[-1] = [item - offset for item in configlist[-1]]
+236
+237 if r_start[rep] is None:
+238 r_start_index.append(0)
+239 else:
+240 try:
+241 r_start_index.append(configlist[-1].index(r_start[rep]))
+242 except ValueError:
+243 raise Exception('Config %d not in file with range [%d, %d]' % (
+244 r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
+245
+246 if r_stop[rep] is None:
+247 r_stop_index.append(len(configlist[-1]) - 1)
+248 else:
+249 try:
+250 r_stop_index.append(configlist[-1].index(r_stop[rep]))
+251 except ValueError:
+252 raise Exception('Config %d not in file with range [%d, %d]' % (
+253 r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
+254
+255 for k in range(nrw):
+256 deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
+257
+258 if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
+259 raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
+260 stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
+261 if np.any([step != 1 for step in stepsizes]):
+262 warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
+263
+264 print(',', nrw, 'reweighting factors with', nsrc, 'sources')
+265 result = []
+266 idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
+267
+268 for t in range(nrw):
+269 result.append(Obs(deltas[t], rep_names, idl=idl))
+270 return result
536def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
-537 """Read the topologial charge based on openQCD gradient flow measurements.
-538
-539 Parameters
-540 ----------
-541 path : str
-542 path of the measurement files
-543 prefix : str
-544 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
-545 Ignored if file names are passed explicitly via keyword files.
-546 c : double
-547 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
-548 dtr_cnfg : int
-549 (optional) parameter that specifies the number of measurements
-550 between two configs.
-551 If it is not set, the distance between two measurements
-552 in the file is assumed to be the distance between two configurations.
-553 steps : int
-554 (optional) Distance between two configurations in units of trajectories /
-555 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
-556 version : str
-557 Either openQCD or sfqcd, depending on the data.
-558 L : int
-559 spatial length of the lattice in L/a.
-560 HAS to be set if version != sfqcd, since openQCD does not provide
-561 this in the header
-562 r_start : list
-563 list which contains the first config to be read for each replicum.
-564 r_stop : list
-565 list which contains the last config to be read for each replicum.
-566 files : list
-567 specify the exact files that need to be read
-568 from path, practical if e.g. only one replicum is needed
-569 postfix : str
-570 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
-571 names : list
-572 Alternative labeling for replicas/ensembles.
-573 Has to have the appropriate length.
-574 Zeuthen_flow : bool
-575 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
-576 for version=='sfqcd' If False, the Wilson flow is used.
-577 integer_charge : bool
-578 If True, the charge is rounded towards the nearest integer on each config.
-579
-580 Returns
-581 -------
-582 result : Obs
-583 Read topological charge
-584 """
-585
-586 return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
+ 560def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
+561 """Read the topologial charge based on openQCD gradient flow measurements.
+562
+563 Parameters
+564 ----------
+565 path : str
+566 path of the measurement files
+567 prefix : str
+568 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
+569 Ignored if file names are passed explicitly via keyword files.
+570 c : double
+571 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
+572 dtr_cnfg : int
+573 (optional) parameter that specifies the number of measurements
+574 between two configs.
+575 If it is not set, the distance between two measurements
+576 in the file is assumed to be the distance between two configurations.
+577 steps : int
+578 (optional) Distance between two configurations in units of trajectories /
+579 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+580 version : str
+581 Either openQCD or sfqcd, depending on the data.
+582 L : int
+583 spatial length of the lattice in L/a.
+584 HAS to be set if version != sfqcd, since openQCD does not provide
+585 this in the header
+586 r_start : list
+587 list which contains the first config to be read for each replicum.
+588 r_stop : list
+589 list which contains the last config to be read for each replicum.
+590 files : list
+591 specify the exact files that need to be read
+592 from path, practical if e.g. only one replicum is needed
+593 postfix : str
+594 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+595 names : list
+596 Alternative labeling for replicas/ensembles.
+597 Has to have the appropriate length.
+598 Zeuthen_flow : bool
+599 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+600 for version=='sfqcd' If False, the Wilson flow is used.
+601 integer_charge : bool
+602 If True, the charge is rounded towards the nearest integer on each config.
+603
+604 Returns
+605 -------
+606 result : Obs
+607 Read topological charge
+608 """
+609
+610 return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
589def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
-590 """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
-591
-592 Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step.
-593
-594 Parameters
-595 ----------
-596 path : str
-597 path of the measurement files
-598 prefix : str
-599 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
-600 Ignored if file names are passed explicitly via keyword files.
-601 c : double
-602 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
-603 dtr_cnfg : int
-604 (optional) parameter that specifies the number of measurements
-605 between two configs.
-606 If it is not set, the distance between two measurements
-607 in the file is assumed to be the distance between two configurations.
-608 steps : int
-609 (optional) Distance between two configurations in units of trajectories /
-610 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
-611 r_start : list
-612 list which contains the first config to be read for each replicum.
-613 r_stop : list
-614 list which contains the last config to be read for each replicum.
-615 files : list
-616 specify the exact files that need to be read
-617 from path, practical if e.g. only one replicum is needed
-618 names : list
-619 Alternative labeling for replicas/ensembles.
-620 Has to have the appropriate length.
-621 postfix : str
-622 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
-623 Zeuthen_flow : bool
-624 (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
-625 """
-626
-627 if c != 0.3:
-628 raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
-629
-630 plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
-631 C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
-632 L = plaq.tag["L"]
-633 T = plaq.tag["T"]
-634
-635 if T != L:
-636 raise Exception("The required lattice norm is only implemented for T=L at the moment.")
-637
-638 if Zeuthen_flow is not True:
-639 raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
-640
-641 t = (c * L) ** 2 / 8
-642
-643 normdict = {4: 0.012341170468270,
-644 6: 0.010162691462430,
-645 8: 0.009031614807931,
-646 10: 0.008744966371393,
-647 12: 0.008650917856809,
-648 14: 8.611154391267955E-03,
-649 16: 0.008591758449508,
-650 20: 0.008575359627103,
-651 24: 0.008569387847540,
-652 28: 8.566803713382559E-03,
-653 32: 0.008565541650006,
-654 40: 8.564480684962046E-03,
-655 48: 8.564098025073460E-03,
-656 64: 8.563853943383087E-03}
-657
-658 return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
+ 613def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
+614 """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
+615
+616 Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step.
+617
+618 Parameters
+619 ----------
+620 path : str
+621 path of the measurement files
+622 prefix : str
+623 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
+624 Ignored if file names are passed explicitly via keyword files.
+625 c : double
+626 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
+627 dtr_cnfg : int
+628 (optional) parameter that specifies the number of measurements
+629 between two configs.
+630 If it is not set, the distance between two measurements
+631 in the file is assumed to be the distance between two configurations.
+632 steps : int
+633 (optional) Distance between two configurations in units of trajectories /
+634 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+635 r_start : list
+636 list which contains the first config to be read for each replicum.
+637 r_stop : list
+638 list which contains the last config to be read for each replicum.
+639 files : list
+640 specify the exact files that need to be read
+641 from path, practical if e.g. only one replicum is needed
+642 names : list
+643 Alternative labeling for replicas/ensembles.
+644 Has to have the appropriate length.
+645 postfix : str
+646 postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
+647 Zeuthen_flow : bool
+648 (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
+649 """
+650
+651 if c != 0.3:
+652 raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
+653
+654 plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
+655 C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
+656 L = plaq.tag["L"]
+657 T = plaq.tag["T"]
+658
+659 if T != L:
+660 raise Exception("The required lattice norm is only implemented for T=L at the moment.")
+661
+662 if Zeuthen_flow is not True:
+663 raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
+664
+665 t = (c * L) ** 2 / 8
+666
+667 normdict = {4: 0.012341170468270,
+668 6: 0.010162691462430,
+669 8: 0.009031614807931,
+670 10: 0.008744966371393,
+671 12: 0.008650917856809,
+672 14: 8.611154391267955E-03,
+673 16: 0.008591758449508,
+674 20: 0.008575359627103,
+675 24: 0.008569387847540,
+676 28: 8.566803713382559E-03,
+677 32: 0.008565541650006,
+678 40: 8.564480684962046E-03,
+679 48: 8.564098025073460E-03,
+680 64: 8.563853943383087E-03}
+681
+682 return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
938def qtop_projection(qtop, target=0):
-939 """Returns the projection to the topological charge sector defined by target.
-940
-941 Parameters
-942 ----------
-943 path : Obs
-944 Topological charge.
-945 target : int
-946 Specifies the topological sector to be reweighted to (default 0)
-947
-948 Returns
-949 -------
-950 reto : Obs
-951 projection to the topological charge sector defined by target
-952 """
-953 if qtop.reweighted:
-954 raise Exception('You can not use a reweighted observable for reweighting!')
-955
-956 proj_qtop = []
-957 for n in qtop.deltas:
-958 proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
-959
-960 reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
-961 return reto
+ 954def qtop_projection(qtop, target=0):
+955 """Returns the projection to the topological charge sector defined by target.
+956
+957 Parameters
+958 ----------
+959 path : Obs
+960 Topological charge.
+961 target : int
+962 Specifies the topological sector to be reweighted to (default 0)
+963
+964 Returns
+965 -------
+966 reto : Obs
+967 projection to the topological charge sector defined by target
+968 """
+969 if qtop.reweighted:
+970 raise Exception('You can not use a reweighted observable for reweighting!')
+971
+972 proj_qtop = []
+973 for n in qtop.deltas:
+974 proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
+975
+976 reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
+977 return reto
964def read_qtop_sector(path, prefix, c, target=0, **kwargs):
- 965 """Constructs reweighting factors to a specified topological sector.
- 966
- 967 Parameters
- 968 ----------
- 969 path : str
- 970 path of the measurement files
- 971 prefix : str
- 972 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
- 973 c : double
- 974 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
- 975 target : int
- 976 Specifies the topological sector to be reweighted to (default 0)
- 977 dtr_cnfg : int
- 978 (optional) parameter that specifies the number of trajectories
- 979 between two configs.
- 980 if it is not set, the distance between two measurements
- 981 in the file is assumed to be the distance between two configurations.
- 982 steps : int
- 983 (optional) Distance between two configurations in units of trajectories /
- 984 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- 985 version : str
- 986 version string of the openQCD (sfqcd) version used to create
- 987 the ensemble. Default is 2.0. May also be set to sfqcd.
- 988 L : int
- 989 spatial length of the lattice in L/a.
- 990 HAS to be set if version != sfqcd, since openQCD does not provide
- 991 this in the header
- 992 r_start : list
- 993 offset of the first ensemble, making it easier to match
- 994 later on with other Obs
- 995 r_stop : list
- 996 last configurations that need to be read (per replicum)
- 997 files : list
- 998 specify the exact files that need to be read
- 999 from path, practical if e.g. only one replicum is needed
-1000 names : list
-1001 Alternative labeling for replicas/ensembles.
-1002 Has to have the appropriate length
-1003 Zeuthen_flow : bool
-1004 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
-1005 for version=='sfqcd' If False, the Wilson flow is used.
-1006
-1007 Returns
-1008 -------
-1009 reto : Obs
-1010 projection to the topological charge sector defined by target
-1011 """
-1012
-1013 if not isinstance(target, int):
-1014 raise Exception("'target' has to be an integer.")
-1015
-1016 kwargs['integer_charge'] = True
-1017 qtop = read_qtop(path, prefix, c, **kwargs)
-1018
-1019 return qtop_projection(qtop, target=target)
+ 980def read_qtop_sector(path, prefix, c, target=0, **kwargs):
+ 981 """Constructs reweighting factors to a specified topological sector.
+ 982
+ 983 Parameters
+ 984 ----------
+ 985 path : str
+ 986 path of the measurement files
+ 987 prefix : str
+ 988 prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
+ 989 c : double
+ 990 Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
+ 991 target : int
+ 992 Specifies the topological sector to be reweighted to (default 0)
+ 993 dtr_cnfg : int
+ 994 (optional) parameter that specifies the number of trajectories
+ 995 between two configs.
+ 996 if it is not set, the distance between two measurements
+ 997 in the file is assumed to be the distance between two configurations.
+ 998 steps : int
+ 999 (optional) Distance between two configurations in units of trajectories /
+1000 cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
+1001 version : str
+1002 version string of the openQCD (sfqcd) version used to create
+1003 the ensemble. Default is 2.0. May also be set to sfqcd.
+1004 L : int
+1005 spatial length of the lattice in L/a.
+1006 HAS to be set if version != sfqcd, since openQCD does not provide
+1007 this in the header
+1008 r_start : list
+1009 offset of the first ensemble, making it easier to match
+1010 later on with other Obs
+1011 r_stop : list
+1012 last configurations that need to be read (per replicum)
+1013 files : list
+1014 specify the exact files that need to be read
+1015 from path, practical if e.g. only one replicum is needed
+1016 names : list
+1017 Alternative labeling for replicas/ensembles.
+1018 Has to have the appropriate length
+1019 Zeuthen_flow : bool
+1020 (optional) If True, the Zeuthen flow is used for Qtop. Only possible
+1021 for version=='sfqcd' If False, the Wilson flow is used.
+1022
+1023 Returns
+1024 -------
+1025 reto : Obs
+1026 projection to the topological charge sector defined by target
+1027 """
+1028
+1029 if not isinstance(target, int):
+1030 raise Exception("'target' has to be an integer.")
+1031
+1032 kwargs['integer_charge'] = True
+1033 qtop = read_qtop(path, prefix, c, **kwargs)
+1034
+1035 return qtop_projection(qtop, target=target)
1022def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
-1023 """
-1024 Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
-1025
-1026 Parameters
-1027 ----------
-1028 path : str
-1029 The directory to search for the files in.
-1030 prefix : str
-1031 The prefix to match the files against.
-1032 qc : str
-1033 The quark combination extension to match the files against.
-1034 corr : str
-1035 The correlator to extract data for.
-1036 sep : str, optional
-1037 The separator to use when parsing the replika names.
-1038 **kwargs
-1039 Additional keyword arguments. The following keyword arguments are recognized:
-1040
-1041 - names (List[str]): A list of names to use for the replicas.
-1042
-1043 Returns
-1044 -------
-1045 Corr
-1046 A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
-1047 or
-1048 CObs
-1049 A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
-1050
-1051
-1052 Raises
-1053 ------
-1054 FileNotFoundError
-1055 If no files matching the specified prefix and quark combination extension are found in the specified directory.
-1056 IOError
-1057 If there is an error reading a file.
-1058 struct.error
-1059 If there is an error unpacking binary data.
-1060 """
-1061
-1062 found = []
-1063 files = []
-1064 names = []
-1065
-1066 if "names" in kwargs:
-1067 names = kwargs.get("names")
-1068
-1069 for (dirpath, dirnames, filenames) in os.walk(path + "/"):
-1070 found.extend(filenames)
-1071 break
-1072
-1073 for f in found:
-1074 if fnmatch.fnmatch(f, prefix + "*.ms5_xsf_" + qc + ".dat"):
-1075 files.append(f)
-1076 if "names" not in kwargs:
-1077 if not sep == "":
-1078 se = f.split(".")[0]
-1079 for s in f.split(".")[1:-1]:
-1080 se += "." + s
-1081 names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
-1082 else:
-1083 names.append(prefix)
-1084
-1085 names = sorted(names)
-1086 files = sorted(files)
-1087
-1088 cnfgs = []
-1089 realsamples = []
-1090 imagsamples = []
-1091 repnum = 0
-1092 for file in files:
-1093 with open(path + "/" + file, "rb") as fp:
+ 1038def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
+1039 """
+1040 Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
+1041
+1042 Parameters
+1043 ----------
+1044 path : str
+1045 The directory to search for the files in.
+1046 prefix : str
+1047 The prefix to match the files against.
+1048 qc : str
+1049 The quark combination extension to match the files against.
+1050 corr : str
+1051 The correlator to extract data for.
+1052 sep : str, optional
+1053 The separator to use when parsing the replika names.
+1054 **kwargs
+1055 Additional keyword arguments. The following keyword arguments are recognized:
+1056
+1057 - names (List[str]): A list of names to use for the replicas.
+1058
+1059 Returns
+1060 -------
+1061 Corr
+1062 A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
+1063 or
+1064 CObs
+1065 A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
+1066
+1067
+1068 Raises
+1069 ------
+1070 FileNotFoundError
+1071 If no files matching the specified prefix and quark combination extension are found in the specified directory.
+1072 IOError
+1073 If there is an error reading a file.
+1074 struct.error
+1075 If there is an error unpacking binary data.
+1076 """
+1077
+1078 # found = []
+1079 files = []
+1080 names = []
+1081
+1082 # test if the input is correct
+1083 if qc not in ['dd', 'ud', 'du', 'uu']:
+1084 raise Exception("Unknown quark conbination!")
+1085
+1086 if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
+1087 raise Exception("Unknown correlator!")
+1088
+1089 if "files" in kwargs:
+1090 known_files = kwargs.get("files")
+1091 else:
+1092 known_files = []
+1093 files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
1094
-1095 t = fp.read(8)
-1096 kappa = struct.unpack('d', t)[0]
-1097 t = fp.read(8)
-1098 csw = struct.unpack('d', t)[0]
-1099 t = fp.read(8)
-1100 dF = struct.unpack('d', t)[0]
-1101 t = fp.read(8)
-1102 zF = struct.unpack('d', t)[0]
-1103
-1104 t = fp.read(4)
-1105 tmax = struct.unpack('i', t)[0]
-1106 t = fp.read(4)
-1107 bnd = struct.unpack('i', t)[0]
-1108
-1109 placesBI = ["gS", "gP",
-1110 "gA", "gV",
-1111 "gVt", "lA",
-1112 "lV", "lVt",
-1113 "lT", "lTt"]
-1114 placesBB = ["g1", "l1"]
-1115
-1116 # the chunks have the following structure:
-1117 # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
-1118
-1119 chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
-1120 packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
-1121 cnfgs.append([])
-1122 realsamples.append([])
-1123 imagsamples.append([])
-1124 for t in range(tmax):
-1125 realsamples[repnum].append([])
-1126 imagsamples[repnum].append([])
-1127
-1128 while True:
-1129 cnfgt = fp.read(chunksize)
-1130 if not cnfgt:
-1131 break
-1132 asascii = struct.unpack(packstr, cnfgt)
-1133 cnfg = asascii[0]
-1134 cnfgs[repnum].append(cnfg)
-1135
-1136 if corr not in placesBB:
-1137 tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
-1138 else:
-1139 tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
+1095 if "names" in kwargs:
+1096 names = kwargs.get("names")
+1097 else:
+1098 for f in files:
+1099 if not sep == "":
+1100 se = f.split(".")[0]
+1101 for s in f.split(".")[1:-2]:
+1102 se += "." + s
+1103 names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
+1104 else:
+1105 names.append(prefix)
+1106
+1107 names = sorted(names)
+1108 files = sorted(files)
+1109
+1110 cnfgs = []
+1111 realsamples = []
+1112 imagsamples = []
+1113 repnum = 0
+1114 for file in files:
+1115 with open(path + "/" + file, "rb") as fp:
+1116
+1117 t = fp.read(8)
+1118 kappa = struct.unpack('d', t)[0]
+1119 t = fp.read(8)
+1120 csw = struct.unpack('d', t)[0]
+1121 t = fp.read(8)
+1122 dF = struct.unpack('d', t)[0]
+1123 t = fp.read(8)
+1124 zF = struct.unpack('d', t)[0]
+1125
+1126 t = fp.read(4)
+1127 tmax = struct.unpack('i', t)[0]
+1128 t = fp.read(4)
+1129 bnd = struct.unpack('i', t)[0]
+1130
+1131 placesBI = ["gS", "gP",
+1132 "gA", "gV",
+1133 "gVt", "lA",
+1134 "lV", "lVt",
+1135 "lT", "lTt"]
+1136 placesBB = ["g1", "l1"]
+1137
+1138 # the chunks have the following structure:
+1139 # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
1140
-1141 corrres = [[], []]
-1142 for i in range(len(tmpcorr)):
-1143 corrres[i % 2].append(tmpcorr[i])
-1144 for t in range(int(len(tmpcorr) / 2)):
-1145 realsamples[repnum][t].append(corrres[0][t])
-1146 for t in range(int(len(tmpcorr) / 2)):
-1147 imagsamples[repnum][t].append(corrres[1][t])
-1148 repnum += 1
+1141 chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
+1142 packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
+1143 cnfgs.append([])
+1144 realsamples.append([])
+1145 imagsamples.append([])
+1146 for t in range(tmax):
+1147 realsamples[repnum].append([])
+1148 imagsamples[repnum].append([])
1149
-1150 s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
-1151 for rep in range(1, repnum):
-1152 s += ", " + str(len(realsamples[rep][t]))
-1153 s += " samples"
-1154 print(s)
-1155 print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
-1156
-1157 # we have the data now... but we need to re format the whole thing and put it into Corr objects.
-1158
-1159 compObs = []
-1160
-1161 for t in range(int(len(tmpcorr) / 2)):
-1162 compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
-1163 Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
-1164
-1165 if len(compObs) == 1:
-1166 return compObs[0]
-1167 else:
-1168 return Corr(compObs)
+1150 while True:
+1151 cnfgt = fp.read(chunksize)
+1152 if not cnfgt:
+1153 break
+1154 asascii = struct.unpack(packstr, cnfgt)
+1155 cnfg = asascii[0]
+1156 cnfgs[repnum].append(cnfg)
+1157
+1158 if corr not in placesBB:
+1159 tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
+1160 else:
+1161 tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
+1162
+1163 corrres = [[], []]
+1164 for i in range(len(tmpcorr)):
+1165 corrres[i % 2].append(tmpcorr[i])
+1166 for t in range(int(len(tmpcorr) / 2)):
+1167 realsamples[repnum][t].append(corrres[0][t])
+1168 for t in range(int(len(tmpcorr) / 2)):
+1169 imagsamples[repnum][t].append(corrres[1][t])
+1170 repnum += 1
+1171
+1172 s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
+1173 for rep in range(1, repnum):
+1174 s += ", " + str(len(realsamples[rep][t]))
+1175 s += " samples"
+1176 print(s)
+1177 print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
+1178
+1179 # we have the data now... but we need to re format the whole thing and put it into Corr objects.
+1180
+1181 compObs = []
+1182
+1183 for t in range(int(len(tmpcorr) / 2)):
+1184 compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
+1185 Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
+1186
+1187 if len(compObs) == 1:
+1188 return compObs[0]
+1189 else:
+1190 return Corr(compObs)