pyerrors.input.openQCD

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

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 read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
570def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
571    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
572
573    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.
574
575    Parameters
576    ----------
577    path : str
578        path of the measurement files
579    prefix : str
580        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
581        Ignored if file names are passed explicitly via keyword files.
582    c : double
583        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
584    dtr_cnfg : int
585        (optional) parameter that specifies the number of measurements
586        between two configs.
587        If it is not set, the distance between two measurements
588        in the file is assumed to be the distance between two configurations.
589    steps : int
590        (optional) Distance between two configurations in units of trajectories /
591         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
592    r_start : list
593        list which contains the first config to be read for each replicum.
594    r_stop : list
595        list which contains the last config to be read for each replicum.
596    files : list
597        specify the exact files that need to be read
598        from path, practical if e.g. only one replicum is needed
599    names : list
600        Alternative labeling for replicas/ensembles.
601        Has to have the appropriate length.
602    Zeuthen_flow : bool
603        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
604    """
605
606    if c != 0.3:
607        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
608
609    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)
610    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)
611    L = plaq.tag["L"]
612    T = plaq.tag["T"]
613
614    if T != L:
615        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
616
617    if Zeuthen_flow is not True:
618        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
619
620    t = (c * L) ** 2 / 8
621
622    normdict = {4: 0.012341170468270,
623                6: 0.010162691462430,
624                8: 0.009031614807931,
625                10: 0.008744966371393,
626                12: 0.008650917856809,
627                14: 8.611154391267955E-03,
628                16: 0.008591758449508,
629                20: 0.008575359627103,
630                24: 0.008569387847540,
631                28: 8.566803713382559E-03,
632                32: 0.008565541650006,
633                40: 8.564480684962046E-03,
634                48: 8.564098025073460E-03,
635                64: 8.563853943383087E-03}
636
637    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]

Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.

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.

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
  • 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 the coupling. If False, the Wilson flow is used.
def qtop_projection(qtop, target=0):
905def qtop_projection(qtop, target=0):
906    """Returns the projection to the topological charge sector defined by target.
907
908    Parameters
909    ----------
910    path : Obs
911        Topological charge.
912    target : int
913        Specifies the topological sector to be reweighted to (default 0)
914    """
915    if qtop.reweighted:
916        raise Exception('You can not use a reweighted observable for reweighting!')
917
918    proj_qtop = []
919    for n in qtop.deltas:
920        proj_qtop.append(np.array([1 if round(qtop.value + q) == target else 0 for q in qtop.deltas[n]]))
921
922    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
923    reto.is_merged = qtop.is_merged
924    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):
927def read_qtop_sector(path, prefix, c, target=0, **kwargs):
928    """Constructs reweighting factors to a specified topological sector.
929
930    Parameters
931    ----------
932    path : str
933        path of the measurement files
934    prefix : str
935        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
936    c : double
937        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
938    target : int
939        Specifies the topological sector to be reweighted to (default 0)
940    dtr_cnfg : int
941        (optional) parameter that specifies the number of trajectories
942        between two configs.
943        if it is not set, the distance between two measurements
944        in the file is assumed to be the distance between two configurations.
945    steps : int
946        (optional) Distance between two configurations in units of trajectories /
947         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
948    version : str
949        version string of the openQCD (sfqcd) version used to create
950        the ensemble. Default is 2.0. May also be set to sfqcd.
951    L : int
952        spatial length of the lattice in L/a.
953        HAS to be set if version != sfqcd, since openQCD does not provide
954        this in the header
955    r_start : list
956        offset of the first ensemble, making it easier to match
957        later on with other Obs
958    r_stop : list
959        last configurations that need to be read (per replicum)
960    files : list
961        specify the exact files that need to be read
962        from path, practical if e.g. only one replicum is needed
963    names : list
964        Alternative labeling for replicas/ensembles.
965        Has to have the appropriate length
966    Zeuthen_flow : bool
967        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
968        for version=='sfqcd' If False, the Wilson flow is used.
969    """
970
971    if not isinstance(target, int):
972        raise Exception("'target' has to be an integer.")
973
974    kwargs['integer_charge'] = True
975    qtop = read_qtop(path, prefix, c, **kwargs)
976
977    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.