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(f"Error, directory '{path}' 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    postfix : str
557        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
558    names : list
559        Alternative labeling for replicas/ensembles.
560        Has to have the appropriate length.
561    Zeuthen_flow : bool
562        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
563        for version=='sfqcd' If False, the Wilson flow is used.
564    integer_charge : bool
565        If True, the charge is rounded towards the nearest integer on each config.
566    """
567
568    return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
569
570
571def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
572    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
573
574    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.
575
576    Parameters
577    ----------
578    path : str
579        path of the measurement files
580    prefix : str
581        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
582        Ignored if file names are passed explicitly via keyword files.
583    c : double
584        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
585    dtr_cnfg : int
586        (optional) parameter that specifies the number of measurements
587        between two configs.
588        If it is not set, the distance between two measurements
589        in the file is assumed to be the distance between two configurations.
590    steps : int
591        (optional) Distance between two configurations in units of trajectories /
592         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
593    r_start : list
594        list which contains the first config to be read for each replicum.
595    r_stop : list
596        list which contains the last config to be read for each replicum.
597    files : list
598        specify the exact files that need to be read
599        from path, practical if e.g. only one replicum is needed
600    names : list
601        Alternative labeling for replicas/ensembles.
602        Has to have the appropriate length.
603    postfix : str
604        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
605    Zeuthen_flow : bool
606        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
607    """
608
609    if c != 0.3:
610        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
611
612    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)
613    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)
614    L = plaq.tag["L"]
615    T = plaq.tag["T"]
616
617    if T != L:
618        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
619
620    if Zeuthen_flow is not True:
621        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
622
623    t = (c * L) ** 2 / 8
624
625    normdict = {4: 0.012341170468270,
626                6: 0.010162691462430,
627                8: 0.009031614807931,
628                10: 0.008744966371393,
629                12: 0.008650917856809,
630                14: 8.611154391267955E-03,
631                16: 0.008591758449508,
632                20: 0.008575359627103,
633                24: 0.008569387847540,
634                28: 8.566803713382559E-03,
635                32: 0.008565541650006,
636                40: 8.564480684962046E-03,
637                48: 8.564098025073460E-03,
638                64: 8.563853943383087E-03}
639
640    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
641
642
643def _read_flow_obs(path, prefix, c, dtr_cnfg=1, version="openQCD", obspos=0, sum_t=True, **kwargs):
644    """Read a flow observable based on openQCD gradient flow measurements.
645
646    Parameters
647    ----------
648    path : str
649        path of the measurement files
650    prefix : str
651        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
652        Ignored if file names are passed explicitly via keyword files.
653    c : double
654        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
655    dtr_cnfg : int
656        (optional) parameter that specifies the number of measurements
657        between two configs.
658        If it is not set, the distance between two measurements
659        in the file is assumed to be the distance between two configurations.
660    steps : int
661        (optional) Distance between two configurations in units of trajectories /
662         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
663    version : str
664        Either openQCD or sfqcd, depending on the data.
665    obspos : int
666        position of the obeservable in the measurement file. Only relevant for sfqcd files.
667    sum_t : bool
668        If true sum over all timeslices, if false only take the value at T/2.
669    L : int
670        spatial length of the lattice in L/a.
671        HAS to be set if version != sfqcd, since openQCD does not provide
672        this in the header
673    r_start : list
674        list which contains the first config to be read for each replicum.
675    r_stop : list
676        list which contains the last config to be read for each replicum.
677    files : list
678        specify the exact files that need to be read
679        from path, practical if e.g. only one replicum is needed
680    names : list
681        Alternative labeling for replicas/ensembles.
682        Has to have the appropriate length.
683    postfix : str
684        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
685    Zeuthen_flow : bool
686        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
687        for version=='sfqcd' If False, the Wilson flow is used.
688    integer_charge : bool
689        If True, the charge is rounded towards the nearest integer on each config.
690    """
691    known_versions = ["openQCD", "sfqcd"]
692
693    if version not in known_versions:
694        raise Exception("Unknown openQCD version.")
695    if "steps" in kwargs:
696        steps = kwargs.get("steps")
697    if version == "sfqcd":
698        if "L" in kwargs:
699            supposed_L = kwargs.get("L")
700        else:
701            supposed_L = None
702        postfix = ".gfms.dat"
703    else:
704        if "L" not in kwargs:
705            raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
706        else:
707            L = kwargs.get("L")
708        postfix = ".ms.dat"
709
710    if "postfix" in kwargs:
711        postfix = kwargs.get("postfix")
712
713    if "files" in kwargs:
714        files = kwargs.get("files")
715        postfix = ''
716    else:
717        found = []
718        files = []
719        for (dirpath, dirnames, filenames) in os.walk(path + "/"):
720            found.extend(filenames)
721            break
722        for f in found:
723            if fnmatch.fnmatch(f, prefix + "*" + postfix):
724                files.append(f)
725
726        files = sorted(files)
727
728    if 'r_start' in kwargs:
729        r_start = kwargs.get('r_start')
730        if len(r_start) != len(files):
731            raise Exception('r_start does not match number of replicas')
732        r_start = [o if o else None for o in r_start]
733    else:
734        r_start = [None] * len(files)
735
736    if 'r_stop' in kwargs:
737        r_stop = kwargs.get('r_stop')
738        if len(r_stop) != len(files):
739            raise Exception('r_stop does not match number of replicas')
740    else:
741        r_stop = [None] * len(files)
742    rep_names = []
743
744    zeuthen = kwargs.get('Zeuthen_flow', False)
745    if zeuthen and version not in ['sfqcd']:
746        raise Exception('Zeuthen flow can only be used for version==sfqcd')
747
748    r_start_index = []
749    r_stop_index = []
750    deltas = []
751    configlist = []
752    if not zeuthen:
753        obspos += 8
754    for rep, file in enumerate(files):
755        with open(path + "/" + file, "rb") as fp:
756
757            Q = []
758            traj_list = []
759            if version in ['sfqcd']:
760                t = fp.read(12)
761                header = struct.unpack('<iii', t)
762                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)
763                ncs = header[1]  # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's
764                tmax = header[2]  # lattice T/a
765
766                t = fp.read(12)
767                Ls = struct.unpack('<iii', t)
768                if (Ls[0] == Ls[1] and Ls[1] == Ls[2]):
769                    L = Ls[0]
770                    if not (supposed_L == L) and supposed_L:
771                        raise Exception("It seems the length given in the header and by you contradict each other")
772                else:
773                    raise Exception("Found more than one spatial length in header!")
774
775                t = fp.read(16)
776                header2 = struct.unpack('<dd', t)
777                tol = header2[0]
778                cmax = header2[1]  # highest value of c used
779
780                if c > cmax:
781                    raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol))
782
783                if (zthfl == 2):
784                    nfl = 2  # number of flows
785                else:
786                    nfl = 1
787                iobs = 8 * nfl  # number of flow observables calculated
788
789                while True:
790                    t = fp.read(4)
791                    if (len(t) < 4):
792                        break
793                    traj_list.append(struct.unpack('i', t)[0])   # trajectory number when measurement was done
794
795                    for j in range(ncs + 1):
796                        for i in range(iobs):
797                            t = fp.read(8 * tmax)
798                            if (i == obspos):  # determines the flow observable -> i=0 <-> Zeuthen flow
799                                Q.append(struct.unpack('d' * tmax, t))
800
801            else:
802                t = fp.read(12)
803                header = struct.unpack('<iii', t)
804                # step size in integration steps "dnms"
805                dn = header[0]
806                # number of measurements, so "ntot"/dn
807                nn = header[1]
808                # lattice T/a
809                tmax = header[2]
810
811                t = fp.read(8)
812                eps = struct.unpack('d', t)[0]
813
814                while True:
815                    t = fp.read(4)
816                    if (len(t) < 4):
817                        break
818                    traj_list.append(struct.unpack('i', t)[0])
819                    # Wsl
820                    t = fp.read(8 * tmax * (nn + 1))
821                    # Ysl
822                    t = fp.read(8 * tmax * (nn + 1))
823                    # Qsl, which is asked for in this method
824                    t = fp.read(8 * tmax * (nn + 1))
825                    # unpack the array of Qtops,
826                    # on each timeslice t=0,...,tmax-1 and the
827                    # measurement number in = 0...nn (see README.qcd1)
828                    tmpd = struct.unpack('d' * tmax * (nn + 1), t)
829                    Q.append(tmpd)
830
831        if len(np.unique(np.diff(traj_list))) != 1:
832            raise Exception("Irregularities in stepsize found")
833        else:
834            if 'steps' in kwargs:
835                if steps != traj_list[1] - traj_list[0]:
836                    raise Exception("steps and the found stepsize are not the same")
837            else:
838                steps = traj_list[1] - traj_list[0]
839
840        configlist.append([tr // steps // dtr_cnfg for tr in traj_list])
841        if configlist[-1][0] > 1:
842            offset = configlist[-1][0] - 1
843            warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % (
844                offset, offset * steps))
845            configlist[-1] = [item - offset for item in configlist[-1]]
846
847        if r_start[rep] is None:
848            r_start_index.append(0)
849        else:
850            try:
851                r_start_index.append(configlist[-1].index(r_start[rep]))
852            except ValueError:
853                raise Exception('Config %d not in file with range [%d, %d]' % (
854                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
855
856        if r_stop[rep] is None:
857            r_stop_index.append(len(configlist[-1]) - 1)
858        else:
859            try:
860                r_stop_index.append(configlist[-1].index(r_stop[rep]))
861            except ValueError:
862                raise Exception('Config %d not in file with range [%d, %d]' % (
863                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
864
865        if version in ['sfqcd']:
866            cstepsize = cmax / ncs
867            index_aim = round(c / cstepsize)
868        else:
869            t_aim = (c * L) ** 2 / 8
870            index_aim = round(t_aim / eps / dn)
871
872        Q_sum = []
873        for i, item in enumerate(Q):
874            if sum_t is True:
875                Q_sum.append([sum(item[current:current + tmax])
876                             for current in range(0, len(item), tmax)])
877            else:
878                Q_sum.append([item[int(tmax / 2)]])
879        Q_top = []
880        if version in ['sfqcd']:
881            for i in range(len(Q_sum) // (ncs + 1)):
882                Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0])
883        else:
884            for i in range(len(Q) // dtr_cnfg):
885                Q_top.append(Q_sum[dtr_cnfg * i][index_aim])
886        if len(Q_top) != len(traj_list) // dtr_cnfg:
887            raise Exception("qtops and traj_list dont have the same length")
888
889        if kwargs.get('integer_charge', False):
890            Q_top = [round(q) for q in Q_top]
891
892        truncated_file = file[:-len(postfix)]
893
894        if "names" not in kwargs:
895            try:
896                idx = truncated_file.index('r')
897            except Exception:
898                if "names" not in kwargs:
899                    raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
900            ens_name = truncated_file[:idx]
901            rep_names.append(ens_name + '|' + truncated_file[idx:])
902        else:
903            names = kwargs.get("names")
904            rep_names = names
905        deltas.append(Q_top)
906
907    idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))]
908    deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))]
909    result = Obs(deltas, rep_names, idl=idl)
910    result.tag = {"T": tmax - 1,
911                  "L": L}
912    return result
913
914
915def qtop_projection(qtop, target=0):
916    """Returns the projection to the topological charge sector defined by target.
917
918    Parameters
919    ----------
920    path : Obs
921        Topological charge.
922    target : int
923        Specifies the topological sector to be reweighted to (default 0)
924    """
925    if qtop.reweighted:
926        raise Exception('You can not use a reweighted observable for reweighting!')
927
928    proj_qtop = []
929    for n in qtop.deltas:
930        proj_qtop.append(np.array([1 if round(qtop.value + q) == target else 0 for q in qtop.deltas[n]]))
931
932    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
933    reto.is_merged = qtop.is_merged
934    return reto
935
936
937def read_qtop_sector(path, prefix, c, target=0, **kwargs):
938    """Constructs reweighting factors to a specified topological sector.
939
940    Parameters
941    ----------
942    path : str
943        path of the measurement files
944    prefix : str
945        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
946    c : double
947        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
948    target : int
949        Specifies the topological sector to be reweighted to (default 0)
950    dtr_cnfg : int
951        (optional) parameter that specifies the number of trajectories
952        between two configs.
953        if it is not set, the distance between two measurements
954        in the file is assumed to be the distance between two configurations.
955    steps : int
956        (optional) Distance between two configurations in units of trajectories /
957         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
958    version : str
959        version string of the openQCD (sfqcd) version used to create
960        the ensemble. Default is 2.0. May also be set to sfqcd.
961    L : int
962        spatial length of the lattice in L/a.
963        HAS to be set if version != sfqcd, since openQCD does not provide
964        this in the header
965    r_start : list
966        offset of the first ensemble, making it easier to match
967        later on with other Obs
968    r_stop : list
969        last configurations that need to be read (per replicum)
970    files : list
971        specify the exact files that need to be read
972        from path, practical if e.g. only one replicum is needed
973    names : list
974        Alternative labeling for replicas/ensembles.
975        Has to have the appropriate length
976    Zeuthen_flow : bool
977        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
978        for version=='sfqcd' If False, the Wilson flow is used.
979    """
980
981    if not isinstance(target, int):
982        raise Exception("'target' has to be an integer.")
983
984    kwargs['integer_charge'] = True
985    qtop = read_qtop(path, prefix, c, **kwargs)
986
987    return qtop_projection(qtop, target=target)