pyerrors.input.openQCD

View Source
  0import os
  1import fnmatch
  2import re
  3import struct
  4import warnings
  5import numpy as np  # Thinly-wrapped numpy
  6import matplotlib.pyplot as plt
  7from matplotlib import gridspec
  8from ..obs import Obs
  9from ..fits import fit_lin
 10
 11
 12def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
 13    """Read rwms format from given folder structure. Returns a list of length nrw
 14
 15    Parameters
 16    ----------
 17    path : str
 18        path that contains the data files
 19    prefix : str
 20        all files in path that start with prefix are considered as input files.
 21        May be used together postfix to consider only special file endings.
 22        Prefix is ignored, if the keyword 'files' is used.
 23    version : str
 24        version of openQCD, default 2.0
 25    names : list
 26        list of names that is assigned to the data according according
 27        to the order in the file list. Use careful, if you do not provide file names!
 28    r_start : list
 29        list which contains the first config to be read for each replicum
 30    r_stop : list
 31        list which contains the last config to be read for each replicum
 32    r_step : int
 33        integer that defines a fixed step size between two measurements (in units of configs)
 34        If not given, r_step=1 is assumed.
 35    postfix : str
 36        postfix of the file to read, e.g. '.ms1' for openQCD-files
 37    files : list
 38        list which contains the filenames to be read. No automatic detection of
 39        files performed if given.
 40    print_err : bool
 41        Print additional information that is useful for debugging.
 42    """
 43    known_oqcd_versions = ['1.4', '1.6', '2.0']
 44    if not (version in known_oqcd_versions):
 45        raise Exception('Unknown openQCD version defined!')
 46    print("Working with openQCD version " + version)
 47    if 'postfix' in kwargs:
 48        postfix = kwargs.get('postfix')
 49    else:
 50        postfix = ''
 51    ls = []
 52    for (dirpath, dirnames, filenames) in os.walk(path):
 53        ls.extend(filenames)
 54        break
 55
 56    if not ls:
 57        raise Exception('Error, directory not found')
 58    if 'files' in kwargs:
 59        ls = kwargs.get('files')
 60    else:
 61        for exc in ls:
 62            if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'):
 63                ls = list(set(ls) - set([exc]))
 64        if len(ls) > 1:
 65            ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
 66    replica = len(ls)
 67
 68    if 'r_start' in kwargs:
 69        r_start = kwargs.get('r_start')
 70        if len(r_start) != replica:
 71            raise Exception('r_start does not match number of replicas')
 72        r_start = [o if o else None for o in r_start]
 73    else:
 74        r_start = [None] * replica
 75
 76    if 'r_stop' in kwargs:
 77        r_stop = kwargs.get('r_stop')
 78        if len(r_stop) != replica:
 79            raise Exception('r_stop does not match number of replicas')
 80    else:
 81        r_stop = [None] * replica
 82
 83    if 'r_step' in kwargs:
 84        r_step = kwargs.get('r_step')
 85    else:
 86        r_step = 1
 87
 88    print('Read reweighting factors from', prefix[:-1], ',',
 89          replica, 'replica', end='')
 90
 91    if names is None:
 92        rep_names = []
 93        for entry in ls:
 94            truncated_entry = entry
 95            suffixes = [".dat", ".rwms", ".ms1"]
 96            for suffix in suffixes:
 97                if truncated_entry.endswith(suffix):
 98                    truncated_entry = truncated_entry[0:-len(suffix)]
 99            idx = truncated_entry.index('r')
100            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
101    else:
102        rep_names = names
103
104    print_err = 0
105    if 'print_err' in kwargs:
106        print_err = 1
107        print()
108
109    deltas = []
110
111    configlist = []
112    r_start_index = []
113    r_stop_index = []
114
115    for rep in range(replica):
116        tmp_array = []
117        with open(path + '/' + ls[rep], 'rb') as fp:
118
119            t = fp.read(4)  # number of reweighting factors
120            if rep == 0:
121                nrw = struct.unpack('i', t)[0]
122                if version == '2.0':
123                    nrw = int(nrw / 2)
124                for k in range(nrw):
125                    deltas.append([])
126            else:
127                if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
128                    raise Exception('Error: different number of reweighting factors for replicum', rep)
129
130            for k in range(nrw):
131                tmp_array.append([])
132
133            # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
134            nfct = []
135            if version in ['1.6', '2.0']:
136                for i in range(nrw):
137                    t = fp.read(4)
138                    nfct.append(struct.unpack('i', t)[0])
139            else:
140                for i in range(nrw):
141                    nfct.append(1)
142
143            nsrc = []
144            for i in range(nrw):
145                t = fp.read(4)
146                nsrc.append(struct.unpack('i', t)[0])
147            if version == '2.0':
148                if not struct.unpack('i', fp.read(4))[0] == 0:
149                    print('something is wrong!')
150
151            configlist.append([])
152            while 0 < 1:
153                t = fp.read(4)
154                if len(t) < 4:
155                    break
156                config_no = struct.unpack('i', t)[0]
157                configlist[-1].append(config_no)
158                for i in range(nrw):
159                    if(version == '2.0'):
160                        tmpd = _read_array_openQCD2(fp)
161                        tmpd = _read_array_openQCD2(fp)
162                        tmp_rw = tmpd['arr']
163                        tmp_nfct = 1.0
164                        for j in range(tmpd['n'][0]):
165                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
166                            if print_err:
167                                print(config_no, i, j,
168                                      np.mean(np.exp(-np.asarray(tmp_rw[j]))),
169                                      np.std(np.exp(-np.asarray(tmp_rw[j]))))
170                                print('Sources:',
171                                      np.exp(-np.asarray(tmp_rw[j])))
172                                print('Partial factor:', tmp_nfct)
173                    elif version == '1.6' or version == '1.4':
174                        tmp_nfct = 1.0
175                        for j in range(nfct[i]):
176                            t = fp.read(8 * nsrc[i])
177                            t = fp.read(8 * nsrc[i])
178                            tmp_rw = struct.unpack('d' * nsrc[i], t)
179                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
180                            if print_err:
181                                print(config_no, i, j,
182                                      np.mean(np.exp(-np.asarray(tmp_rw))),
183                                      np.std(np.exp(-np.asarray(tmp_rw))))
184                                print('Sources:', np.exp(-np.asarray(tmp_rw)))
185                                print('Partial factor:', tmp_nfct)
186                    tmp_array[i].append(tmp_nfct)
187
188            diffmeas = configlist[-1][-1] - configlist[-1][-2]
189            configlist[-1] = [item // diffmeas for item in configlist[-1]]
190            if configlist[-1][0] > 1 and diffmeas > 1:
191                warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
192                offset = configlist[-1][0] - 1
193                configlist[-1] = [item - offset for item in configlist[-1]]
194
195            if r_start[rep] is None:
196                r_start_index.append(0)
197            else:
198                try:
199                    r_start_index.append(configlist[-1].index(r_start[rep]))
200                except ValueError:
201                    raise Exception('Config %d not in file with range [%d, %d]' % (
202                        r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
203
204            if r_stop[rep] is None:
205                r_stop_index.append(len(configlist[-1]) - 1)
206            else:
207                try:
208                    r_stop_index.append(configlist[-1].index(r_stop[rep]))
209                except ValueError:
210                    raise Exception('Config %d not in file with range [%d, %d]' % (
211                        r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
212
213            for k in range(nrw):
214                deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
215
216    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
217        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
218    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
219    if np.any([step != 1 for step in stepsizes]):
220        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
221
222    print(',', nrw, 'reweighting factors with', nsrc, 'sources')
223    result = []
224    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
225
226    for t in range(nrw):
227        result.append(Obs(deltas[t], rep_names, idl=idl))
228    return result
229
230
231def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
232    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
233
234    It is assumed that all boundary effects have
235    sufficiently decayed at x0=xmin.
236    The data around the zero crossing of t^2<E> - 0.3
237    is fitted with a linear function
238    from which the exact root is extracted.
239
240    It is assumed that one measurement is performed for each config.
241    If this is not the case, the resulting idl, as well as the handling
242    of r_start, r_stop and r_step is wrong and the user has to correct
243    this in the resulting observable.
244
245    Parameters
246    ----------
247    path : str
248        Path to .ms.dat files
249    prefix : str
250        Ensemble prefix
251    dtr_read : int
252        Determines how many trajectories should be skipped
253        when reading the ms.dat files.
254        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
255    xmin : int
256        First timeslice where the boundary
257        effects have sufficiently decayed.
258    spatial_extent : int
259        spatial extent of the lattice, required for normalization.
260    fit_range : int
261        Number of data points left and right of the zero
262        crossing to be included in the linear fit. (Default: 5)
263    r_start : list
264        list which contains the first config to be read for each replicum.
265    r_stop : list
266        list which contains the last config to be read for each replicum.
267    r_step : int
268        integer that defines a fixed step size between two measurements (in units of configs)
269        If not given, r_step=1 is assumed.
270    plaquette : bool
271        If true extract the plaquette estimate of t0 instead.
272    names : list
273        list of names that is assigned to the data according according
274        to the order in the file list. Use careful, if you do not provide file names!
275    files : list
276        list which contains the filenames to be read. No automatic detection of
277        files performed if given.
278    plot_fit : bool
279        If true, the fit for the extraction of t0 is shown together with the data.
280    assume_thermalization : bool
281        If True: If the first record divided by the distance between two measurements is larger than
282        1, it is assumed that this is due to thermalization and the first measurement belongs
283        to the first config (default).
284        If False: The config numbers are assumed to be traj_number // difference
285    """
286
287    ls = []
288    for (dirpath, dirnames, filenames) in os.walk(path):
289        ls.extend(filenames)
290        break
291
292    if not ls:
293        raise Exception('Error, directory not found')
294
295    if 'files' in kwargs:
296        ls = kwargs.get('files')
297    else:
298        for exc in ls:
299            if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'):
300                ls = list(set(ls) - set([exc]))
301        if len(ls) > 1:
302            ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
303    replica = len(ls)
304
305    if 'r_start' in kwargs:
306        r_start = kwargs.get('r_start')
307        if len(r_start) != replica:
308            raise Exception('r_start does not match number of replicas')
309        r_start = [o if o else None for o in r_start]
310    else:
311        r_start = [None] * replica
312
313    if 'r_stop' in kwargs:
314        r_stop = kwargs.get('r_stop')
315        if len(r_stop) != replica:
316            raise Exception('r_stop does not match number of replicas')
317    else:
318        r_stop = [None] * replica
319
320    if 'r_step' in kwargs:
321        r_step = kwargs.get('r_step')
322    else:
323        r_step = 1
324
325    print('Extract t0 from', prefix, ',', replica, 'replica')
326
327    if 'names' in kwargs:
328        rep_names = kwargs.get('names')
329    else:
330        rep_names = []
331        for entry in ls:
332            truncated_entry = entry.split('.')[0]
333            idx = truncated_entry.index('r')
334            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
335
336    Ysum = []
337
338    configlist = []
339    r_start_index = []
340    r_stop_index = []
341
342    for rep in range(replica):
343
344        with open(path + '/' + ls[rep], 'rb') as fp:
345            t = fp.read(12)
346            header = struct.unpack('iii', t)
347            if rep == 0:
348                dn = header[0]
349                nn = header[1]
350                tmax = header[2]
351            elif dn != header[0] or nn != header[1] or tmax != header[2]:
352                raise Exception('Replica parameters do not match.')
353
354            t = fp.read(8)
355            if rep == 0:
356                eps = struct.unpack('d', t)[0]
357                print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
358            elif eps != struct.unpack('d', t)[0]:
359                raise Exception('Values for eps do not match among replica.')
360
361            Ysl = []
362
363            configlist.append([])
364            while 0 < 1:
365                t = fp.read(4)
366                if(len(t) < 4):
367                    break
368                nc = struct.unpack('i', t)[0]
369                configlist[-1].append(nc)
370
371                t = fp.read(8 * tmax * (nn + 1))
372                if kwargs.get('plaquette'):
373                    if nc % dtr_read == 0:
374                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
375                t = fp.read(8 * tmax * (nn + 1))
376                if not kwargs.get('plaquette'):
377                    if nc % dtr_read == 0:
378                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
379                t = fp.read(8 * tmax * (nn + 1))
380
381        Ysum.append([])
382        for i, item in enumerate(Ysl):
383            Ysum[-1].append([np.mean(item[current + xmin:
384                             current + tmax - xmin])
385                            for current in range(0, len(item), tmax)])
386
387        diffmeas = configlist[-1][-1] - configlist[-1][-2]
388        configlist[-1] = [item // diffmeas for item in configlist[-1]]
389        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
390            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
391            offset = configlist[-1][0] - 1
392            configlist[-1] = [item - offset for item in configlist[-1]]
393
394        if r_start[rep] is None:
395            r_start_index.append(0)
396        else:
397            try:
398                r_start_index.append(configlist[-1].index(r_start[rep]))
399            except ValueError:
400                raise Exception('Config %d not in file with range [%d, %d]' % (
401                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
402
403        if r_stop[rep] is None:
404            r_stop_index.append(len(configlist[-1]) - 1)
405        else:
406            try:
407                r_stop_index.append(configlist[-1].index(r_stop[rep]))
408            except ValueError:
409                raise Exception('Config %d not in file with range [%d, %d]' % (
410                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
411
412    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
413        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
414    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
415    if np.any([step != 1 for step in stepsizes]):
416        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
417
418    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
419    t2E_dict = {}
420    for n in range(nn + 1):
421        samples = []
422        for nrep, rep in enumerate(Ysum):
423            samples.append([])
424            for cnfg in rep:
425                samples[-1].append(cnfg[n])
426            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
427        new_obs = Obs(samples, rep_names, idl=idl)
428        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
429
430    zero_crossing = np.argmax(np.array(
431        [o.value for o in t2E_dict.values()]) > 0.0)
432
433    x = list(t2E_dict.keys())[zero_crossing - fit_range:
434                              zero_crossing + fit_range]
435    y = list(t2E_dict.values())[zero_crossing - fit_range:
436                                zero_crossing + fit_range]
437    [o.gamma_method() for o in y]
438
439    fit_result = fit_lin(x, y)
440
441    if kwargs.get('plot_fit'):
442        plt.figure()
443        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
444        ax0 = plt.subplot(gs[0])
445        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
446        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
447        [o.gamma_method() for o in ymore]
448        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
449        xplot = np.linspace(np.min(x), np.max(x))
450        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
451        [yi.gamma_method() for yi in yplot]
452        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
453        retval = (-fit_result[0] / fit_result[1])
454        retval.gamma_method()
455        ylim = ax0.get_ylim()
456        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
457        ax0.set_ylim(ylim)
458        ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
459        xlim = ax0.get_xlim()
460
461        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
462        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
463        ax1 = plt.subplot(gs[1])
464        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
465        ax1.tick_params(direction='out')
466        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
467        ax1.axhline(y=0.0, ls='--', color='k')
468        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
469        ax1.set_xlim(xlim)
470        ax1.set_ylabel('Residuals')
471        ax1.set_xlabel(r'$t/a^2$')
472
473        plt.draw()
474    return -fit_result[0] / fit_result[1]
475
476
477def _parse_array_openQCD2(d, n, size, wa, quadrupel=False):
478    arr = []
479    if d == 2:
480        for i in range(n[0]):
481            tmp = wa[i * n[1]:(i + 1) * n[1]]
482            if quadrupel:
483                tmp2 = []
484                for j in range(0, len(tmp), 2):
485                    tmp2.append(tmp[j])
486                arr.append(tmp2)
487            else:
488                arr.append(np.asarray(tmp))
489
490    else:
491        raise Exception('Only two-dimensional arrays supported!')
492
493    return arr
494
495
496def _read_array_openQCD2(fp):
497    t = fp.read(4)
498    d = struct.unpack('i', t)[0]
499    t = fp.read(4 * d)
500    n = struct.unpack('%di' % (d), t)
501    t = fp.read(4)
502    size = struct.unpack('i', t)[0]
503    if size == 4:
504        types = 'i'
505    elif size == 8:
506        types = 'd'
507    elif size == 16:
508        types = 'dd'
509    else:
510        raise Exception("Type for size '" + str(size) + "' not known.")
511    m = n[0]
512    for i in range(1, d):
513        m *= n[i]
514
515    t = fp.read(m * size)
516    tmp = struct.unpack('%d%s' % (m, types), t)
517
518    arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True)
519    return {'d': d, 'n': n, 'size': size, 'arr': arr}
520
521
522def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
523    """Read the topologial charge based on openQCD gradient flow measurements.
524
525    Parameters
526    ----------
527    path : str
528        path of the measurement files
529    prefix : str
530        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
531        Ignored if file names are passed explicitly via keyword files.
532    c : double
533        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
534    dtr_cnfg : int
535        (optional) parameter that specifies the number of measurements
536        between two configs.
537        If it is not set, the distance between two measurements
538        in the file is assumed to be the distance between two configurations.
539    steps : int
540        (optional) Distance between two configurations in units of trajectories /
541         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
542    version : str
543        Either openQCD or sfqcd, depending on the data.
544    L : int
545        spatial length of the lattice in L/a.
546        HAS to be set if version != sfqcd, since openQCD does not provide
547        this in the header
548    r_start : list
549        list which contains the first config to be read for each replicum.
550    r_stop : list
551        list which contains the last config to be read for each replicum.
552    files : list
553        specify the exact files that need to be read
554        from path, practical if e.g. only one replicum is needed
555    names : list
556        Alternative labeling for replicas/ensembles.
557        Has to have the appropriate length.
558    Zeuthen_flow : bool
559        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
560        for version=='sfqcd' If False, the Wilson flow is used.
561    integer_charge : bool
562        If True, the charge is rounded towards the nearest integer on each config.
563    """
564    known_versions = ["openQCD", "sfqcd"]
565
566    if version not in known_versions:
567        raise Exception("Unknown openQCD version.")
568    if "steps" in kwargs:
569        steps = kwargs.get("steps")
570    if version == "sfqcd":
571        if "L" in kwargs:
572            supposed_L = kwargs.get("L")
573        else:
574            supposed_L = None
575        postfix = ".gfms.dat"
576    else:
577        if "L" not in kwargs:
578            raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
579        else:
580            L = kwargs.get("L")
581        postfix = ".ms.dat"
582
583    if "files" in kwargs:
584        files = kwargs.get("files")
585        postfix = ''
586    else:
587        found = []
588        files = []
589        for (dirpath, dirnames, filenames) in os.walk(path + "/"):
590            found.extend(filenames)
591            break
592        for f in found:
593            if fnmatch.fnmatch(f, prefix + "*" + postfix):
594                files.append(f)
595
596    if 'r_start' in kwargs:
597        r_start = kwargs.get('r_start')
598        if len(r_start) != len(files):
599            raise Exception('r_start does not match number of replicas')
600        r_start = [o if o else None for o in r_start]
601    else:
602        r_start = [None] * len(files)
603
604    if 'r_stop' in kwargs:
605        r_stop = kwargs.get('r_stop')
606        if len(r_stop) != len(files):
607            raise Exception('r_stop does not match number of replicas')
608    else:
609        r_stop = [None] * len(files)
610    rep_names = []
611
612    zeuthen = kwargs.get('Zeuthen_flow', False)
613    if zeuthen and version not in ['sfqcd']:
614        raise Exception('Zeuthen flow can only be used for version==sfqcd')
615
616    r_start_index = []
617    r_stop_index = []
618    deltas = []
619    configlist = []
620    for rep, file in enumerate(files):
621        with open(path + "/" + file, "rb") as fp:
622
623            Q = []
624            traj_list = []
625            if version in ['sfqcd']:
626                if zeuthen:
627                    obspos = 0
628                else:
629                    obspos = 8
630                t = fp.read(12)
631                header = struct.unpack('<iii', t)
632                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)
633                ncs = header[1]  # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's
634                tmax = header[2]  # lattice T/a
635
636                t = fp.read(12)
637                Ls = struct.unpack('<iii', t)
638                if(Ls[0] == Ls[1] and Ls[1] == Ls[2]):
639                    L = Ls[0]
640                    if not (supposed_L == L) and supposed_L:
641                        raise Exception("It seems the length given in the header and by you contradict each other")
642                else:
643                    raise Exception("Found more than one spatial length in header!")
644
645                t = fp.read(16)
646                header2 = struct.unpack('<dd', t)
647                tol = header2[0]
648                cmax = header2[1]  # highest value of c used
649
650                if c > cmax:
651                    raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol))
652
653                if(zthfl == 2):
654                    nfl = 2  # number of flows
655                else:
656                    nfl = 1
657                iobs = 8 * nfl  # number of flow observables calculated
658
659                while 0 < 1:
660                    t = fp.read(4)
661                    if(len(t) < 4):
662                        break
663                    traj_list.append(struct.unpack('i', t)[0])   # trajectory number when measurement was done
664
665                    for j in range(ncs + 1):
666                        for i in range(iobs):
667                            t = fp.read(8 * tmax)
668                            if (i == obspos):  # determines the flow observable -> i=0 <-> Zeuthen flow
669                                Q.append(struct.unpack('d' * tmax, t))
670
671            else:
672                t = fp.read(12)
673                header = struct.unpack('<iii', t)
674                # step size in integration steps "dnms"
675                dn = header[0]
676                # number of measurements, so "ntot"/dn
677                nn = header[1]
678                # lattice T/a
679                tmax = header[2]
680
681                t = fp.read(8)
682                eps = struct.unpack('d', t)[0]
683
684                while 0 < 1:
685                    t = fp.read(4)
686                    if(len(t) < 4):
687                        break
688                    traj_list.append(struct.unpack('i', t)[0])
689                    # Wsl
690                    t = fp.read(8 * tmax * (nn + 1))
691                    # Ysl
692                    t = fp.read(8 * tmax * (nn + 1))
693                    # Qsl, which is asked for in this method
694                    t = fp.read(8 * tmax * (nn + 1))
695                    # unpack the array of Qtops,
696                    # on each timeslice t=0,...,tmax-1 and the
697                    # measurement number in = 0...nn (see README.qcd1)
698                    tmpd = struct.unpack('d' * tmax * (nn + 1), t)
699                    Q.append(tmpd)
700
701        if len(np.unique(np.diff(traj_list))) != 1:
702            raise Exception("Irregularities in stepsize found")
703        else:
704            if 'steps' in kwargs:
705                if steps != traj_list[1] - traj_list[0]:
706                    raise Exception("steps and the found stepsize are not the same")
707            else:
708                steps = traj_list[1] - traj_list[0]
709
710        configlist.append([tr // steps // dtr_cnfg for tr in traj_list])
711        if configlist[-1][0] > 1:
712            offset = configlist[-1][0] - 1
713            warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % (
714                offset, offset * steps))
715            configlist[-1] = [item - offset for item in configlist[-1]]
716
717        if r_start[rep] is None:
718            r_start_index.append(0)
719        else:
720            try:
721                r_start_index.append(configlist[-1].index(r_start[rep]))
722            except ValueError:
723                raise Exception('Config %d not in file with range [%d, %d]' % (
724                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
725
726        if r_stop[rep] is None:
727            r_stop_index.append(len(configlist[-1]) - 1)
728        else:
729            try:
730                r_stop_index.append(configlist[-1].index(r_stop[rep]))
731            except ValueError:
732                raise Exception('Config %d not in file with range [%d, %d]' % (
733                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
734
735        if version in ['sfqcd']:
736            cstepsize = cmax / ncs
737            index_aim = round(c / cstepsize)
738        else:
739            t_aim = (c * L) ** 2 / 8
740            index_aim = round(t_aim / eps / dn)
741
742        Q_sum = []
743        for i, item in enumerate(Q):
744            Q_sum.append([sum(item[current:current + tmax])
745                         for current in range(0, len(item), tmax)])
746        Q_top = []
747        if version in ['sfqcd']:
748            for i in range(len(Q_sum) // (ncs + 1)):
749                Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0])
750        else:
751            for i in range(len(Q) // dtr_cnfg):
752                Q_top.append(Q_sum[dtr_cnfg * i][index_aim])
753        if len(Q_top) != len(traj_list) // dtr_cnfg:
754            raise Exception("qtops and traj_list dont have the same length")
755
756        if kwargs.get('integer_charge', False):
757            Q_top = [round(q) for q in Q_top]
758
759        truncated_file = file[:-len(postfix)]
760
761        if "names" not in kwargs:
762            try:
763                idx = truncated_file.index('r')
764            except Exception:
765                if "names" not in kwargs:
766                    raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
767            ens_name = truncated_file[:idx]
768            rep_names.append(ens_name + '|' + truncated_file[idx:])
769        else:
770            names = kwargs.get("names")
771            rep_names = names
772        deltas.append(Q_top)
773
774    idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))]
775    deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))]
776    result = Obs(deltas, rep_names, idl=idl)
777    return result
778
779
780def qtop_projection(qtop, target=0):
781    """Returns the projection to the topological charge sector defined by target.
782
783    Parameters
784    ----------
785    path : Obs
786        Topological charge.
787    target : int
788        Specifies the topological sector to be reweighted to (default 0)
789    """
790    if qtop.reweighted:
791        raise Exception('You can not use a reweighted observable for reweighting!')
792
793    proj_qtop = []
794    for n in qtop.deltas:
795        proj_qtop.append(np.array([1 if round(qtop.value + q) == target else 0 for q in qtop.deltas[n]]))
796
797    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
798    reto.is_merged = qtop.is_merged
799    return reto
800
801
802def read_qtop_sector(path, prefix, c, target=0, **kwargs):
803    """Constructs reweighting factors to a specified topological sector.
804
805    Parameters
806    ----------
807    path : str
808        path of the measurement files
809    prefix : str
810        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
811    c : double
812        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
813    target : int
814        Specifies the topological sector to be reweighted to (default 0)
815    dtr_cnfg : int
816        (optional) parameter that specifies the number of trajectories
817        between two configs.
818        if it is not set, the distance between two measurements
819        in the file is assumed to be the distance between two configurations.
820    steps : int
821        (optional) Distance between two configurations in units of trajectories /
822         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
823    version : str
824        version string of the openQCD (sfqcd) version used to create
825        the ensemble. Default is 2.0. May also be set to sfqcd.
826    L : int
827        spatial length of the lattice in L/a.
828        HAS to be set if version != sfqcd, since openQCD does not provide
829        this in the header
830    r_start : list
831        offset of the first ensemble, making it easier to match
832        later on with other Obs
833    r_stop : list
834        last configurations that need to be read (per replicum)
835    files : list
836        specify the exact files that need to be read
837        from path, practical if e.g. only one replicum is needed
838    names : list
839        Alternative labeling for replicas/ensembles.
840        Has to have the appropriate length
841    Zeuthen_flow : bool
842        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
843        for version=='sfqcd' If False, the Wilson flow is used.
844    """
845
846    if not isinstance(target, int):
847        raise Exception("'target' has to be an integer.")
848
849    kwargs['integer_charge'] = True
850    qtop = read_qtop(path, prefix, c, **kwargs)
851
852    return qtop_projection(qtop, target=target)
#   def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
View Source
 13def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
 14    """Read rwms format from given folder structure. Returns a list of length nrw
 15
 16    Parameters
 17    ----------
 18    path : str
 19        path that contains the data files
 20    prefix : str
 21        all files in path that start with prefix are considered as input files.
 22        May be used together postfix to consider only special file endings.
 23        Prefix is ignored, if the keyword 'files' is used.
 24    version : str
 25        version of openQCD, default 2.0
 26    names : list
 27        list of names that is assigned to the data according according
 28        to the order in the file list. Use careful, if you do not provide file names!
 29    r_start : list
 30        list which contains the first config to be read for each replicum
 31    r_stop : list
 32        list which contains the last config to be read for each replicum
 33    r_step : int
 34        integer that defines a fixed step size between two measurements (in units of configs)
 35        If not given, r_step=1 is assumed.
 36    postfix : str
 37        postfix of the file to read, e.g. '.ms1' for openQCD-files
 38    files : list
 39        list which contains the filenames to be read. No automatic detection of
 40        files performed if given.
 41    print_err : bool
 42        Print additional information that is useful for debugging.
 43    """
 44    known_oqcd_versions = ['1.4', '1.6', '2.0']
 45    if not (version in known_oqcd_versions):
 46        raise Exception('Unknown openQCD version defined!')
 47    print("Working with openQCD version " + version)
 48    if 'postfix' in kwargs:
 49        postfix = kwargs.get('postfix')
 50    else:
 51        postfix = ''
 52    ls = []
 53    for (dirpath, dirnames, filenames) in os.walk(path):
 54        ls.extend(filenames)
 55        break
 56
 57    if not ls:
 58        raise Exception('Error, directory not found')
 59    if 'files' in kwargs:
 60        ls = kwargs.get('files')
 61    else:
 62        for exc in ls:
 63            if not fnmatch.fnmatch(exc, prefix + '*' + postfix + '.dat'):
 64                ls = list(set(ls) - set([exc]))
 65        if len(ls) > 1:
 66            ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
 67    replica = len(ls)
 68
 69    if 'r_start' in kwargs:
 70        r_start = kwargs.get('r_start')
 71        if len(r_start) != replica:
 72            raise Exception('r_start does not match number of replicas')
 73        r_start = [o if o else None for o in r_start]
 74    else:
 75        r_start = [None] * replica
 76
 77    if 'r_stop' in kwargs:
 78        r_stop = kwargs.get('r_stop')
 79        if len(r_stop) != replica:
 80            raise Exception('r_stop does not match number of replicas')
 81    else:
 82        r_stop = [None] * replica
 83
 84    if 'r_step' in kwargs:
 85        r_step = kwargs.get('r_step')
 86    else:
 87        r_step = 1
 88
 89    print('Read reweighting factors from', prefix[:-1], ',',
 90          replica, 'replica', end='')
 91
 92    if names is None:
 93        rep_names = []
 94        for entry in ls:
 95            truncated_entry = entry
 96            suffixes = [".dat", ".rwms", ".ms1"]
 97            for suffix in suffixes:
 98                if truncated_entry.endswith(suffix):
 99                    truncated_entry = truncated_entry[0:-len(suffix)]
100            idx = truncated_entry.index('r')
101            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
102    else:
103        rep_names = names
104
105    print_err = 0
106    if 'print_err' in kwargs:
107        print_err = 1
108        print()
109
110    deltas = []
111
112    configlist = []
113    r_start_index = []
114    r_stop_index = []
115
116    for rep in range(replica):
117        tmp_array = []
118        with open(path + '/' + ls[rep], 'rb') as fp:
119
120            t = fp.read(4)  # number of reweighting factors
121            if rep == 0:
122                nrw = struct.unpack('i', t)[0]
123                if version == '2.0':
124                    nrw = int(nrw / 2)
125                for k in range(nrw):
126                    deltas.append([])
127            else:
128                if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
129                    raise Exception('Error: different number of reweighting factors for replicum', rep)
130
131            for k in range(nrw):
132                tmp_array.append([])
133
134            # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
135            nfct = []
136            if version in ['1.6', '2.0']:
137                for i in range(nrw):
138                    t = fp.read(4)
139                    nfct.append(struct.unpack('i', t)[0])
140            else:
141                for i in range(nrw):
142                    nfct.append(1)
143
144            nsrc = []
145            for i in range(nrw):
146                t = fp.read(4)
147                nsrc.append(struct.unpack('i', t)[0])
148            if version == '2.0':
149                if not struct.unpack('i', fp.read(4))[0] == 0:
150                    print('something is wrong!')
151
152            configlist.append([])
153            while 0 < 1:
154                t = fp.read(4)
155                if len(t) < 4:
156                    break
157                config_no = struct.unpack('i', t)[0]
158                configlist[-1].append(config_no)
159                for i in range(nrw):
160                    if(version == '2.0'):
161                        tmpd = _read_array_openQCD2(fp)
162                        tmpd = _read_array_openQCD2(fp)
163                        tmp_rw = tmpd['arr']
164                        tmp_nfct = 1.0
165                        for j in range(tmpd['n'][0]):
166                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
167                            if print_err:
168                                print(config_no, i, j,
169                                      np.mean(np.exp(-np.asarray(tmp_rw[j]))),
170                                      np.std(np.exp(-np.asarray(tmp_rw[j]))))
171                                print('Sources:',
172                                      np.exp(-np.asarray(tmp_rw[j])))
173                                print('Partial factor:', tmp_nfct)
174                    elif version == '1.6' or version == '1.4':
175                        tmp_nfct = 1.0
176                        for j in range(nfct[i]):
177                            t = fp.read(8 * nsrc[i])
178                            t = fp.read(8 * nsrc[i])
179                            tmp_rw = struct.unpack('d' * nsrc[i], t)
180                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
181                            if print_err:
182                                print(config_no, i, j,
183                                      np.mean(np.exp(-np.asarray(tmp_rw))),
184                                      np.std(np.exp(-np.asarray(tmp_rw))))
185                                print('Sources:', np.exp(-np.asarray(tmp_rw)))
186                                print('Partial factor:', tmp_nfct)
187                    tmp_array[i].append(tmp_nfct)
188
189            diffmeas = configlist[-1][-1] - configlist[-1][-2]
190            configlist[-1] = [item // diffmeas for item in configlist[-1]]
191            if configlist[-1][0] > 1 and diffmeas > 1:
192                warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
193                offset = configlist[-1][0] - 1
194                configlist[-1] = [item - offset for item in configlist[-1]]
195
196            if r_start[rep] is None:
197                r_start_index.append(0)
198            else:
199                try:
200                    r_start_index.append(configlist[-1].index(r_start[rep]))
201                except ValueError:
202                    raise Exception('Config %d not in file with range [%d, %d]' % (
203                        r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
204
205            if r_stop[rep] is None:
206                r_stop_index.append(len(configlist[-1]) - 1)
207            else:
208                try:
209                    r_stop_index.append(configlist[-1].index(r_stop[rep]))
210                except ValueError:
211                    raise Exception('Config %d not in file with range [%d, %d]' % (
212                        r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
213
214            for k in range(nrw):
215                deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
216
217    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
218        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
219    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
220    if np.any([step != 1 for step in stepsizes]):
221        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
222
223    print(',', nrw, 'reweighting factors with', nsrc, 'sources')
224    result = []
225    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
226
227    for t in range(nrw):
228        result.append(Obs(deltas[t], rep_names, idl=idl))
229    return result

Read rwms format from given folder structure. Returns a list of length nrw

Parameters
  • path (str): path that contains the data files
  • prefix (str): all files in path that start with prefix are considered as input files. May be used together postfix to consider only special file endings. Prefix is ignored, if the keyword 'files' is used.
  • version (str): version of openQCD, default 2.0
  • names (list): list of names that is assigned to the data according according to the order in the file list. Use careful, if you do not provide file names!
  • r_start (list): list which contains the first config to be read for each replicum
  • r_stop (list): list which contains the last config to be read for each replicum
  • r_step (int): integer that defines a fixed step size between two measurements (in units of configs) If not given, r_step=1 is assumed.
  • postfix (str): postfix of the file to read, e.g. '.ms1' for openQCD-files
  • files (list): list which contains the filenames to be read. No automatic detection of files performed if given.
  • print_err (bool): Print additional information that is useful for debugging.
#   def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
View Source
232def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
233    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
234
235    It is assumed that all boundary effects have
236    sufficiently decayed at x0=xmin.
237    The data around the zero crossing of t^2<E> - 0.3
238    is fitted with a linear function
239    from which the exact root is extracted.
240
241    It is assumed that one measurement is performed for each config.
242    If this is not the case, the resulting idl, as well as the handling
243    of r_start, r_stop and r_step is wrong and the user has to correct
244    this in the resulting observable.
245
246    Parameters
247    ----------
248    path : str
249        Path to .ms.dat files
250    prefix : str
251        Ensemble prefix
252    dtr_read : int
253        Determines how many trajectories should be skipped
254        when reading the ms.dat files.
255        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
256    xmin : int
257        First timeslice where the boundary
258        effects have sufficiently decayed.
259    spatial_extent : int
260        spatial extent of the lattice, required for normalization.
261    fit_range : int
262        Number of data points left and right of the zero
263        crossing to be included in the linear fit. (Default: 5)
264    r_start : list
265        list which contains the first config to be read for each replicum.
266    r_stop : list
267        list which contains the last config to be read for each replicum.
268    r_step : int
269        integer that defines a fixed step size between two measurements (in units of configs)
270        If not given, r_step=1 is assumed.
271    plaquette : bool
272        If true extract the plaquette estimate of t0 instead.
273    names : list
274        list of names that is assigned to the data according according
275        to the order in the file list. Use careful, if you do not provide file names!
276    files : list
277        list which contains the filenames to be read. No automatic detection of
278        files performed if given.
279    plot_fit : bool
280        If true, the fit for the extraction of t0 is shown together with the data.
281    assume_thermalization : bool
282        If True: If the first record divided by the distance between two measurements is larger than
283        1, it is assumed that this is due to thermalization and the first measurement belongs
284        to the first config (default).
285        If False: The config numbers are assumed to be traj_number // difference
286    """
287
288    ls = []
289    for (dirpath, dirnames, filenames) in os.walk(path):
290        ls.extend(filenames)
291        break
292
293    if not ls:
294        raise Exception('Error, directory not found')
295
296    if 'files' in kwargs:
297        ls = kwargs.get('files')
298    else:
299        for exc in ls:
300            if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'):
301                ls = list(set(ls) - set([exc]))
302        if len(ls) > 1:
303            ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
304    replica = len(ls)
305
306    if 'r_start' in kwargs:
307        r_start = kwargs.get('r_start')
308        if len(r_start) != replica:
309            raise Exception('r_start does not match number of replicas')
310        r_start = [o if o else None for o in r_start]
311    else:
312        r_start = [None] * replica
313
314    if 'r_stop' in kwargs:
315        r_stop = kwargs.get('r_stop')
316        if len(r_stop) != replica:
317            raise Exception('r_stop does not match number of replicas')
318    else:
319        r_stop = [None] * replica
320
321    if 'r_step' in kwargs:
322        r_step = kwargs.get('r_step')
323    else:
324        r_step = 1
325
326    print('Extract t0 from', prefix, ',', replica, 'replica')
327
328    if 'names' in kwargs:
329        rep_names = kwargs.get('names')
330    else:
331        rep_names = []
332        for entry in ls:
333            truncated_entry = entry.split('.')[0]
334            idx = truncated_entry.index('r')
335            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
336
337    Ysum = []
338
339    configlist = []
340    r_start_index = []
341    r_stop_index = []
342
343    for rep in range(replica):
344
345        with open(path + '/' + ls[rep], 'rb') as fp:
346            t = fp.read(12)
347            header = struct.unpack('iii', t)
348            if rep == 0:
349                dn = header[0]
350                nn = header[1]
351                tmax = header[2]
352            elif dn != header[0] or nn != header[1] or tmax != header[2]:
353                raise Exception('Replica parameters do not match.')
354
355            t = fp.read(8)
356            if rep == 0:
357                eps = struct.unpack('d', t)[0]
358                print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
359            elif eps != struct.unpack('d', t)[0]:
360                raise Exception('Values for eps do not match among replica.')
361
362            Ysl = []
363
364            configlist.append([])
365            while 0 < 1:
366                t = fp.read(4)
367                if(len(t) < 4):
368                    break
369                nc = struct.unpack('i', t)[0]
370                configlist[-1].append(nc)
371
372                t = fp.read(8 * tmax * (nn + 1))
373                if kwargs.get('plaquette'):
374                    if nc % dtr_read == 0:
375                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
376                t = fp.read(8 * tmax * (nn + 1))
377                if not kwargs.get('plaquette'):
378                    if nc % dtr_read == 0:
379                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
380                t = fp.read(8 * tmax * (nn + 1))
381
382        Ysum.append([])
383        for i, item in enumerate(Ysl):
384            Ysum[-1].append([np.mean(item[current + xmin:
385                             current + tmax - xmin])
386                            for current in range(0, len(item), tmax)])
387
388        diffmeas = configlist[-1][-1] - configlist[-1][-2]
389        configlist[-1] = [item // diffmeas for item in configlist[-1]]
390        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
391            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
392            offset = configlist[-1][0] - 1
393            configlist[-1] = [item - offset for item in configlist[-1]]
394
395        if r_start[rep] is None:
396            r_start_index.append(0)
397        else:
398            try:
399                r_start_index.append(configlist[-1].index(r_start[rep]))
400            except ValueError:
401                raise Exception('Config %d not in file with range [%d, %d]' % (
402                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
403
404        if r_stop[rep] is None:
405            r_stop_index.append(len(configlist[-1]) - 1)
406        else:
407            try:
408                r_stop_index.append(configlist[-1].index(r_stop[rep]))
409            except ValueError:
410                raise Exception('Config %d not in file with range [%d, %d]' % (
411                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
412
413    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
414        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
415    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
416    if np.any([step != 1 for step in stepsizes]):
417        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
418
419    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
420    t2E_dict = {}
421    for n in range(nn + 1):
422        samples = []
423        for nrep, rep in enumerate(Ysum):
424            samples.append([])
425            for cnfg in rep:
426                samples[-1].append(cnfg[n])
427            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
428        new_obs = Obs(samples, rep_names, idl=idl)
429        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
430
431    zero_crossing = np.argmax(np.array(
432        [o.value for o in t2E_dict.values()]) > 0.0)
433
434    x = list(t2E_dict.keys())[zero_crossing - fit_range:
435                              zero_crossing + fit_range]
436    y = list(t2E_dict.values())[zero_crossing - fit_range:
437                                zero_crossing + fit_range]
438    [o.gamma_method() for o in y]
439
440    fit_result = fit_lin(x, y)
441
442    if kwargs.get('plot_fit'):
443        plt.figure()
444        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
445        ax0 = plt.subplot(gs[0])
446        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
447        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
448        [o.gamma_method() for o in ymore]
449        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
450        xplot = np.linspace(np.min(x), np.max(x))
451        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
452        [yi.gamma_method() for yi in yplot]
453        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
454        retval = (-fit_result[0] / fit_result[1])
455        retval.gamma_method()
456        ylim = ax0.get_ylim()
457        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
458        ax0.set_ylim(ylim)
459        ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
460        xlim = ax0.get_xlim()
461
462        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
463        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
464        ax1 = plt.subplot(gs[1])
465        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
466        ax1.tick_params(direction='out')
467        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
468        ax1.axhline(y=0.0, ls='--', color='k')
469        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
470        ax1.set_xlim(xlim)
471        ax1.set_ylabel('Residuals')
472        ax1.set_xlabel(r'$t/a^2$')
473
474        plt.draw()
475    return -fit_result[0] / fit_result[1]

Extract t0 from given .ms.dat files. Returns t0 as Obs.

It is assumed that all boundary effects have sufficiently decayed at x0=xmin. The data around the zero crossing of t^2 - 0.3 is fitted with a linear function from which the exact root is extracted.

It is assumed that one measurement is performed for each config. If this is not the case, the resulting idl, as well as the handling of r_start, r_stop and r_step is wrong and the user has to correct this in the resulting observable.

Parameters
  • path (str): Path to .ms.dat files
  • prefix (str): Ensemble prefix
  • dtr_read (int): Determines how many trajectories should be skipped when reading the ms.dat files. Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
  • xmin (int): First timeslice where the boundary effects have sufficiently decayed.
  • spatial_extent (int): spatial extent of the lattice, required for normalization.
  • fit_range (int): Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5)
  • r_start (list): list which contains the first config to be read for each replicum.
  • r_stop (list): list which contains the last config to be read for each replicum.
  • r_step (int): integer that defines a fixed step size between two measurements (in units of configs) If not given, r_step=1 is assumed.
  • plaquette (bool): If true extract the plaquette estimate of t0 instead.
  • names (list): list of names that is assigned to the data according according to the order in the file list. Use careful, if you do not provide file names!
  • files (list): list which contains the filenames to be read. No automatic detection of files performed if given.
  • plot_fit (bool): If true, the fit for the extraction of t0 is shown together with the data.
  • assume_thermalization (bool): If True: If the first record divided by the distance between two measurements is larger than 1, it is assumed that this is due to thermalization and the first measurement belongs to the first config (default). If False: The config numbers are assumed to be traj_number // difference
#   def read_qtop(path, prefix, c, dtr_cnfg=1, version='openQCD', **kwargs):
View Source
523def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
524    """Read the topologial charge based on openQCD gradient flow measurements.
525
526    Parameters
527    ----------
528    path : str
529        path of the measurement files
530    prefix : str
531        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
532        Ignored if file names are passed explicitly via keyword files.
533    c : double
534        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
535    dtr_cnfg : int
536        (optional) parameter that specifies the number of measurements
537        between two configs.
538        If it is not set, the distance between two measurements
539        in the file is assumed to be the distance between two configurations.
540    steps : int
541        (optional) Distance between two configurations in units of trajectories /
542         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
543    version : str
544        Either openQCD or sfqcd, depending on the data.
545    L : int
546        spatial length of the lattice in L/a.
547        HAS to be set if version != sfqcd, since openQCD does not provide
548        this in the header
549    r_start : list
550        list which contains the first config to be read for each replicum.
551    r_stop : list
552        list which contains the last config to be read for each replicum.
553    files : list
554        specify the exact files that need to be read
555        from path, practical if e.g. only one replicum is needed
556    names : list
557        Alternative labeling for replicas/ensembles.
558        Has to have the appropriate length.
559    Zeuthen_flow : bool
560        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
561        for version=='sfqcd' If False, the Wilson flow is used.
562    integer_charge : bool
563        If True, the charge is rounded towards the nearest integer on each config.
564    """
565    known_versions = ["openQCD", "sfqcd"]
566
567    if version not in known_versions:
568        raise Exception("Unknown openQCD version.")
569    if "steps" in kwargs:
570        steps = kwargs.get("steps")
571    if version == "sfqcd":
572        if "L" in kwargs:
573            supposed_L = kwargs.get("L")
574        else:
575            supposed_L = None
576        postfix = ".gfms.dat"
577    else:
578        if "L" not in kwargs:
579            raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
580        else:
581            L = kwargs.get("L")
582        postfix = ".ms.dat"
583
584    if "files" in kwargs:
585        files = kwargs.get("files")
586        postfix = ''
587    else:
588        found = []
589        files = []
590        for (dirpath, dirnames, filenames) in os.walk(path + "/"):
591            found.extend(filenames)
592            break
593        for f in found:
594            if fnmatch.fnmatch(f, prefix + "*" + postfix):
595                files.append(f)
596
597    if 'r_start' in kwargs:
598        r_start = kwargs.get('r_start')
599        if len(r_start) != len(files):
600            raise Exception('r_start does not match number of replicas')
601        r_start = [o if o else None for o in r_start]
602    else:
603        r_start = [None] * len(files)
604
605    if 'r_stop' in kwargs:
606        r_stop = kwargs.get('r_stop')
607        if len(r_stop) != len(files):
608            raise Exception('r_stop does not match number of replicas')
609    else:
610        r_stop = [None] * len(files)
611    rep_names = []
612
613    zeuthen = kwargs.get('Zeuthen_flow', False)
614    if zeuthen and version not in ['sfqcd']:
615        raise Exception('Zeuthen flow can only be used for version==sfqcd')
616
617    r_start_index = []
618    r_stop_index = []
619    deltas = []
620    configlist = []
621    for rep, file in enumerate(files):
622        with open(path + "/" + file, "rb") as fp:
623
624            Q = []
625            traj_list = []
626            if version in ['sfqcd']:
627                if zeuthen:
628                    obspos = 0
629                else:
630                    obspos = 8
631                t = fp.read(12)
632                header = struct.unpack('<iii', t)
633                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)
634                ncs = header[1]  # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's
635                tmax = header[2]  # lattice T/a
636
637                t = fp.read(12)
638                Ls = struct.unpack('<iii', t)
639                if(Ls[0] == Ls[1] and Ls[1] == Ls[2]):
640                    L = Ls[0]
641                    if not (supposed_L == L) and supposed_L:
642                        raise Exception("It seems the length given in the header and by you contradict each other")
643                else:
644                    raise Exception("Found more than one spatial length in header!")
645
646                t = fp.read(16)
647                header2 = struct.unpack('<dd', t)
648                tol = header2[0]
649                cmax = header2[1]  # highest value of c used
650
651                if c > cmax:
652                    raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol))
653
654                if(zthfl == 2):
655                    nfl = 2  # number of flows
656                else:
657                    nfl = 1
658                iobs = 8 * nfl  # number of flow observables calculated
659
660                while 0 < 1:
661                    t = fp.read(4)
662                    if(len(t) < 4):
663                        break
664                    traj_list.append(struct.unpack('i', t)[0])   # trajectory number when measurement was done
665
666                    for j in range(ncs + 1):
667                        for i in range(iobs):
668                            t = fp.read(8 * tmax)
669                            if (i == obspos):  # determines the flow observable -> i=0 <-> Zeuthen flow
670                                Q.append(struct.unpack('d' * tmax, t))
671
672            else:
673                t = fp.read(12)
674                header = struct.unpack('<iii', t)
675                # step size in integration steps "dnms"
676                dn = header[0]
677                # number of measurements, so "ntot"/dn
678                nn = header[1]
679                # lattice T/a
680                tmax = header[2]
681
682                t = fp.read(8)
683                eps = struct.unpack('d', t)[0]
684
685                while 0 < 1:
686                    t = fp.read(4)
687                    if(len(t) < 4):
688                        break
689                    traj_list.append(struct.unpack('i', t)[0])
690                    # Wsl
691                    t = fp.read(8 * tmax * (nn + 1))
692                    # Ysl
693                    t = fp.read(8 * tmax * (nn + 1))
694                    # Qsl, which is asked for in this method
695                    t = fp.read(8 * tmax * (nn + 1))
696                    # unpack the array of Qtops,
697                    # on each timeslice t=0,...,tmax-1 and the
698                    # measurement number in = 0...nn (see README.qcd1)
699                    tmpd = struct.unpack('d' * tmax * (nn + 1), t)
700                    Q.append(tmpd)
701
702        if len(np.unique(np.diff(traj_list))) != 1:
703            raise Exception("Irregularities in stepsize found")
704        else:
705            if 'steps' in kwargs:
706                if steps != traj_list[1] - traj_list[0]:
707                    raise Exception("steps and the found stepsize are not the same")
708            else:
709                steps = traj_list[1] - traj_list[0]
710
711        configlist.append([tr // steps // dtr_cnfg for tr in traj_list])
712        if configlist[-1][0] > 1:
713            offset = configlist[-1][0] - 1
714            warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % (
715                offset, offset * steps))
716            configlist[-1] = [item - offset for item in configlist[-1]]
717
718        if r_start[rep] is None:
719            r_start_index.append(0)
720        else:
721            try:
722                r_start_index.append(configlist[-1].index(r_start[rep]))
723            except ValueError:
724                raise Exception('Config %d not in file with range [%d, %d]' % (
725                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
726
727        if r_stop[rep] is None:
728            r_stop_index.append(len(configlist[-1]) - 1)
729        else:
730            try:
731                r_stop_index.append(configlist[-1].index(r_stop[rep]))
732            except ValueError:
733                raise Exception('Config %d not in file with range [%d, %d]' % (
734                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
735
736        if version in ['sfqcd']:
737            cstepsize = cmax / ncs
738            index_aim = round(c / cstepsize)
739        else:
740            t_aim = (c * L) ** 2 / 8
741            index_aim = round(t_aim / eps / dn)
742
743        Q_sum = []
744        for i, item in enumerate(Q):
745            Q_sum.append([sum(item[current:current + tmax])
746                         for current in range(0, len(item), tmax)])
747        Q_top = []
748        if version in ['sfqcd']:
749            for i in range(len(Q_sum) // (ncs + 1)):
750                Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0])
751        else:
752            for i in range(len(Q) // dtr_cnfg):
753                Q_top.append(Q_sum[dtr_cnfg * i][index_aim])
754        if len(Q_top) != len(traj_list) // dtr_cnfg:
755            raise Exception("qtops and traj_list dont have the same length")
756
757        if kwargs.get('integer_charge', False):
758            Q_top = [round(q) for q in Q_top]
759
760        truncated_file = file[:-len(postfix)]
761
762        if "names" not in kwargs:
763            try:
764                idx = truncated_file.index('r')
765            except Exception:
766                if "names" not in kwargs:
767                    raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
768            ens_name = truncated_file[:idx]
769            rep_names.append(ens_name + '|' + truncated_file[idx:])
770        else:
771            names = kwargs.get("names")
772            rep_names = names
773        deltas.append(Q_top)
774
775    idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))]
776    deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))]
777    result = Obs(deltas, rep_names, idl=idl)
778    return result

Read the topologial charge based on openQCD gradient flow measurements.

Parameters
  • path (str): path of the measurement files
  • prefix (str): prefix of the measurement files, e.g. _id0_r0.ms.dat. Ignored if file names are passed explicitly via keyword files.
  • c (double): Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
  • dtr_cnfg (int): (optional) parameter that specifies the number of measurements between two configs. If it is not set, the distance between two measurements in the file is assumed to be the distance between two configurations.
  • steps (int): (optional) Distance between two configurations in units of trajectories / cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
  • version (str): Either openQCD or sfqcd, depending on the data.
  • L (int): spatial length of the lattice in L/a. HAS to be set if version != sfqcd, since openQCD does not provide this in the header
  • r_start (list): list which contains the first config to be read for each replicum.
  • r_stop (list): list which contains the last config to be read for each replicum.
  • files (list): specify the exact files that need to be read from path, practical if e.g. only one replicum is needed
  • names (list): Alternative labeling for replicas/ensembles. Has to have the appropriate length.
  • Zeuthen_flow (bool): (optional) If True, the Zeuthen flow is used for Qtop. Only possible for version=='sfqcd' If False, the Wilson flow is used.
  • integer_charge (bool): If True, the charge is rounded towards the nearest integer on each config.
#   def qtop_projection(qtop, target=0):
View Source
781def qtop_projection(qtop, target=0):
782    """Returns the projection to the topological charge sector defined by target.
783
784    Parameters
785    ----------
786    path : Obs
787        Topological charge.
788    target : int
789        Specifies the topological sector to be reweighted to (default 0)
790    """
791    if qtop.reweighted:
792        raise Exception('You can not use a reweighted observable for reweighting!')
793
794    proj_qtop = []
795    for n in qtop.deltas:
796        proj_qtop.append(np.array([1 if round(qtop.value + q) == target else 0 for q in qtop.deltas[n]]))
797
798    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
799    reto.is_merged = qtop.is_merged
800    return reto

Returns the projection to the topological charge sector defined by target.

Parameters
  • path (Obs): Topological charge.
  • target (int): Specifies the topological sector to be reweighted to (default 0)
#   def read_qtop_sector(path, prefix, c, target=0, **kwargs):
View Source
803def read_qtop_sector(path, prefix, c, target=0, **kwargs):
804    """Constructs reweighting factors to a specified topological sector.
805
806    Parameters
807    ----------
808    path : str
809        path of the measurement files
810    prefix : str
811        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
812    c : double
813        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
814    target : int
815        Specifies the topological sector to be reweighted to (default 0)
816    dtr_cnfg : int
817        (optional) parameter that specifies the number of trajectories
818        between two configs.
819        if it is not set, the distance between two measurements
820        in the file is assumed to be the distance between two configurations.
821    steps : int
822        (optional) Distance between two configurations in units of trajectories /
823         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
824    version : str
825        version string of the openQCD (sfqcd) version used to create
826        the ensemble. Default is 2.0. May also be set to sfqcd.
827    L : int
828        spatial length of the lattice in L/a.
829        HAS to be set if version != sfqcd, since openQCD does not provide
830        this in the header
831    r_start : list
832        offset of the first ensemble, making it easier to match
833        later on with other Obs
834    r_stop : list
835        last configurations that need to be read (per replicum)
836    files : list
837        specify the exact files that need to be read
838        from path, practical if e.g. only one replicum is needed
839    names : list
840        Alternative labeling for replicas/ensembles.
841        Has to have the appropriate length
842    Zeuthen_flow : bool
843        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
844        for version=='sfqcd' If False, the Wilson flow is used.
845    """
846
847    if not isinstance(target, int):
848        raise Exception("'target' has to be an integer.")
849
850    kwargs['integer_charge'] = True
851    qtop = read_qtop(path, prefix, c, **kwargs)
852
853    return qtop_projection(qtop, target=target)

Constructs reweighting factors to a specified topological sector.

Parameters
  • path (str): path of the measurement files
  • prefix (str): prefix of the measurement files, e.g. _id0_r0.ms.dat
  • c (double): Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
  • target (int): Specifies the topological sector to be reweighted to (default 0)
  • dtr_cnfg (int): (optional) parameter that specifies the number of trajectories between two configs. if it is not set, the distance between two measurements in the file is assumed to be the distance between two configurations.
  • steps (int): (optional) Distance between two configurations in units of trajectories / cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
  • version (str): version string of the openQCD (sfqcd) version used to create the ensemble. Default is 2.0. May also be set to sfqcd.
  • L (int): spatial length of the lattice in L/a. HAS to be set if version != sfqcd, since openQCD does not provide this in the header
  • r_start (list): offset of the first ensemble, making it easier to match later on with other Obs
  • r_stop (list): last configurations that need to be read (per replicum)
  • files (list): specify the exact files that need to be read from path, practical if e.g. only one replicum is needed
  • names (list): Alternative labeling for replicas/ensembles. Has to have the appropriate length
  • Zeuthen_flow (bool): (optional) If True, the Zeuthen flow is used for Qtop. Only possible for version=='sfqcd' If False, the Wilson flow is used.