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 0 < 1:
154                t = fp.read(4)
155                if len(t) < 4:
156                    break
157                config_no = struct.unpack('i', t)[0]
158                configlist[-1].append(config_no)
159                for i in range(nrw):
160                    if(version == '2.0'):
161                        tmpd = _read_array_openQCD2(fp)
162                        tmpd = _read_array_openQCD2(fp)
163                        tmp_rw = tmpd['arr']
164                        tmp_nfct = 1.0
165                        for j in range(tmpd['n'][0]):
166                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
167                            if print_err:
168                                print(config_no, i, j,
169                                      np.mean(np.exp(-np.asarray(tmp_rw[j]))),
170                                      np.std(np.exp(-np.asarray(tmp_rw[j]))))
171                                print('Sources:',
172                                      np.exp(-np.asarray(tmp_rw[j])))
173                                print('Partial factor:', tmp_nfct)
174                    elif version == '1.6' or version == '1.4':
175                        tmp_nfct = 1.0
176                        for j in range(nfct[i]):
177                            t = fp.read(8 * nsrc[i])
178                            t = fp.read(8 * nsrc[i])
179                            tmp_rw = struct.unpack('d' * nsrc[i], t)
180                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
181                            if print_err:
182                                print(config_no, i, j,
183                                      np.mean(np.exp(-np.asarray(tmp_rw))),
184                                      np.std(np.exp(-np.asarray(tmp_rw))))
185                                print('Sources:', np.exp(-np.asarray(tmp_rw)))
186                                print('Partial factor:', tmp_nfct)
187                    tmp_array[i].append(tmp_nfct)
188
189            diffmeas = configlist[-1][-1] - configlist[-1][-2]
190            configlist[-1] = [item // diffmeas for item in configlist[-1]]
191            if configlist[-1][0] > 1 and diffmeas > 1:
192                warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
193                offset = configlist[-1][0] - 1
194                configlist[-1] = [item - offset for item in configlist[-1]]
195
196            if r_start[rep] is None:
197                r_start_index.append(0)
198            else:
199                try:
200                    r_start_index.append(configlist[-1].index(r_start[rep]))
201                except ValueError:
202                    raise Exception('Config %d not in file with range [%d, %d]' % (
203                        r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
204
205            if r_stop[rep] is None:
206                r_stop_index.append(len(configlist[-1]) - 1)
207            else:
208                try:
209                    r_stop_index.append(configlist[-1].index(r_stop[rep]))
210                except ValueError:
211                    raise Exception('Config %d not in file with range [%d, %d]' % (
212                        r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
213
214            for k in range(nrw):
215                deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
216
217    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
218        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
219    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
220    if np.any([step != 1 for step in stepsizes]):
221        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
222
223    print(',', nrw, 'reweighting factors with', nsrc, 'sources')
224    result = []
225    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
226
227    for t in range(nrw):
228        result.append(Obs(deltas[t], rep_names, idl=idl))
229    return result
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 0 < 1:
366                t = fp.read(4)
367                if(len(t) < 4):
368                    break
369                nc = struct.unpack('i', t)[0]
370                configlist[-1].append(nc)
371
372                t = fp.read(8 * tmax * (nn + 1))
373                if kwargs.get('plaquette'):
374                    if nc % dtr_read == 0:
375                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
376                t = fp.read(8 * tmax * (nn + 1))
377                if not kwargs.get('plaquette'):
378                    if nc % dtr_read == 0:
379                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
380                t = fp.read(8 * tmax * (nn + 1))
381
382        Ysum.append([])
383        for i, item in enumerate(Ysl):
384            Ysum[-1].append([np.mean(item[current + xmin:
385                             current + tmax - xmin])
386                            for current in range(0, len(item), tmax)])
387
388        diffmeas = configlist[-1][-1] - configlist[-1][-2]
389        configlist[-1] = [item // diffmeas for item in configlist[-1]]
390        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
391            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
392            offset = configlist[-1][0] - 1
393            configlist[-1] = [item - offset for item in configlist[-1]]
394
395        if r_start[rep] is None:
396            r_start_index.append(0)
397        else:
398            try:
399                r_start_index.append(configlist[-1].index(r_start[rep]))
400            except ValueError:
401                raise Exception('Config %d not in file with range [%d, %d]' % (
402                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
403
404        if r_stop[rep] is None:
405            r_stop_index.append(len(configlist[-1]) - 1)
406        else:
407            try:
408                r_stop_index.append(configlist[-1].index(r_stop[rep]))
409            except ValueError:
410                raise Exception('Config %d not in file with range [%d, %d]' % (
411                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
412
413    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
414        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
415    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
416    if np.any([step != 1 for step in stepsizes]):
417        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
418
419    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
420    t2E_dict = {}
421    for n in range(nn + 1):
422        samples = []
423        for nrep, rep in enumerate(Ysum):
424            samples.append([])
425            for cnfg in rep:
426                samples[-1].append(cnfg[n])
427            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
428        new_obs = Obs(samples, rep_names, idl=idl)
429        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
430
431    zero_crossing = np.argmax(np.array(
432        [o.value for o in t2E_dict.values()]) > 0.0)
433
434    x = list(t2E_dict.keys())[zero_crossing - fit_range:
435                              zero_crossing + fit_range]
436    y = list(t2E_dict.values())[zero_crossing - fit_range:
437                                zero_crossing + fit_range]
438    [o.gamma_method() for o in y]
439
440    fit_result = fit_lin(x, y)
441
442    if kwargs.get('plot_fit'):
443        plt.figure()
444        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
445        ax0 = plt.subplot(gs[0])
446        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
447        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
448        [o.gamma_method() for o in ymore]
449        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
450        xplot = np.linspace(np.min(x), np.max(x))
451        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
452        [yi.gamma_method() for yi in yplot]
453        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
454        retval = (-fit_result[0] / fit_result[1])
455        retval.gamma_method()
456        ylim = ax0.get_ylim()
457        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
458        ax0.set_ylim(ylim)
459        ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
460        xlim = ax0.get_xlim()
461
462        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
463        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
464        ax1 = plt.subplot(gs[1])
465        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
466        ax1.tick_params(direction='out')
467        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
468        ax1.axhline(y=0.0, ls='--', color='k')
469        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
470        ax1.set_xlim(xlim)
471        ax1.set_ylabel('Residuals')
472        ax1.set_xlabel(r'$t/a^2$')
473
474        plt.draw()
475    return -fit_result[0] / fit_result[1]
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    known_versions = ["openQCD", "sfqcd"]
566
567    if version not in known_versions:
568        raise Exception("Unknown openQCD version.")
569    if "steps" in kwargs:
570        steps = kwargs.get("steps")
571    if version == "sfqcd":
572        if "L" in kwargs:
573            supposed_L = kwargs.get("L")
574        else:
575            supposed_L = None
576        postfix = ".gfms.dat"
577    else:
578        if "L" not in kwargs:
579            raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
580        else:
581            L = kwargs.get("L")
582        postfix = ".ms.dat"
583
584    if "files" in kwargs:
585        files = kwargs.get("files")
586        postfix = ''
587    else:
588        found = []
589        files = []
590        for (dirpath, dirnames, filenames) in os.walk(path + "/"):
591            found.extend(filenames)
592            break
593        for f in found:
594            if fnmatch.fnmatch(f, prefix + "*" + postfix):
595                files.append(f)
596
597    if 'r_start' in kwargs:
598        r_start = kwargs.get('r_start')
599        if len(r_start) != len(files):
600            raise Exception('r_start does not match number of replicas')
601        r_start = [o if o else None for o in r_start]
602    else:
603        r_start = [None] * len(files)
604
605    if 'r_stop' in kwargs:
606        r_stop = kwargs.get('r_stop')
607        if len(r_stop) != len(files):
608            raise Exception('r_stop does not match number of replicas')
609    else:
610        r_stop = [None] * len(files)
611    rep_names = []
612
613    zeuthen = kwargs.get('Zeuthen_flow', False)
614    if zeuthen and version not in ['sfqcd']:
615        raise Exception('Zeuthen flow can only be used for version==sfqcd')
616
617    r_start_index = []
618    r_stop_index = []
619    deltas = []
620    configlist = []
621    for rep, file in enumerate(files):
622        with open(path + "/" + file, "rb") as fp:
623
624            Q = []
625            traj_list = []
626            if version in ['sfqcd']:
627                if zeuthen:
628                    obspos = 0
629                else:
630                    obspos = 8
631                t = fp.read(12)
632                header = struct.unpack('<iii', t)
633                zthfl = header[0]  # Zeuthen flow -> if it's equal to 2 it means that the Zeuthen flow is also 'measured' (apart from the Wilson flow)
634                ncs = header[1]  # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's
635                tmax = header[2]  # lattice T/a
636
637                t = fp.read(12)
638                Ls = struct.unpack('<iii', t)
639                if(Ls[0] == Ls[1] and Ls[1] == Ls[2]):
640                    L = Ls[0]
641                    if not (supposed_L == L) and supposed_L:
642                        raise Exception("It seems the length given in the header and by you contradict each other")
643                else:
644                    raise Exception("Found more than one spatial length in header!")
645
646                t = fp.read(16)
647                header2 = struct.unpack('<dd', t)
648                tol = header2[0]
649                cmax = header2[1]  # highest value of c used
650
651                if c > cmax:
652                    raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol))
653
654                if(zthfl == 2):
655                    nfl = 2  # number of flows
656                else:
657                    nfl = 1
658                iobs = 8 * nfl  # number of flow observables calculated
659
660                while 0 < 1:
661                    t = fp.read(4)
662                    if(len(t) < 4):
663                        break
664                    traj_list.append(struct.unpack('i', t)[0])   # trajectory number when measurement was done
665
666                    for j in range(ncs + 1):
667                        for i in range(iobs):
668                            t = fp.read(8 * tmax)
669                            if (i == obspos):  # determines the flow observable -> i=0 <-> Zeuthen flow
670                                Q.append(struct.unpack('d' * tmax, t))
671
672            else:
673                t = fp.read(12)
674                header = struct.unpack('<iii', t)
675                # step size in integration steps "dnms"
676                dn = header[0]
677                # number of measurements, so "ntot"/dn
678                nn = header[1]
679                # lattice T/a
680                tmax = header[2]
681
682                t = fp.read(8)
683                eps = struct.unpack('d', t)[0]
684
685                while 0 < 1:
686                    t = fp.read(4)
687                    if(len(t) < 4):
688                        break
689                    traj_list.append(struct.unpack('i', t)[0])
690                    # Wsl
691                    t = fp.read(8 * tmax * (nn + 1))
692                    # Ysl
693                    t = fp.read(8 * tmax * (nn + 1))
694                    # Qsl, which is asked for in this method
695                    t = fp.read(8 * tmax * (nn + 1))
696                    # unpack the array of Qtops,
697                    # on each timeslice t=0,...,tmax-1 and the
698                    # measurement number in = 0...nn (see README.qcd1)
699                    tmpd = struct.unpack('d' * tmax * (nn + 1), t)
700                    Q.append(tmpd)
701
702        if len(np.unique(np.diff(traj_list))) != 1:
703            raise Exception("Irregularities in stepsize found")
704        else:
705            if 'steps' in kwargs:
706                if steps != traj_list[1] - traj_list[0]:
707                    raise Exception("steps and the found stepsize are not the same")
708            else:
709                steps = traj_list[1] - traj_list[0]
710
711        configlist.append([tr // steps // dtr_cnfg for tr in traj_list])
712        if configlist[-1][0] > 1:
713            offset = configlist[-1][0] - 1
714            warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % (
715                offset, offset * steps))
716            configlist[-1] = [item - offset for item in configlist[-1]]
717
718        if r_start[rep] is None:
719            r_start_index.append(0)
720        else:
721            try:
722                r_start_index.append(configlist[-1].index(r_start[rep]))
723            except ValueError:
724                raise Exception('Config %d not in file with range [%d, %d]' % (
725                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
726
727        if r_stop[rep] is None:
728            r_stop_index.append(len(configlist[-1]) - 1)
729        else:
730            try:
731                r_stop_index.append(configlist[-1].index(r_stop[rep]))
732            except ValueError:
733                raise Exception('Config %d not in file with range [%d, %d]' % (
734                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
735
736        if version in ['sfqcd']:
737            cstepsize = cmax / ncs
738            index_aim = round(c / cstepsize)
739        else:
740            t_aim = (c * L) ** 2 / 8
741            index_aim = round(t_aim / eps / dn)
742
743        Q_sum = []
744        for i, item in enumerate(Q):
745            Q_sum.append([sum(item[current:current + tmax])
746                         for current in range(0, len(item), tmax)])
747        Q_top = []
748        if version in ['sfqcd']:
749            for i in range(len(Q_sum) // (ncs + 1)):
750                Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0])
751        else:
752            for i in range(len(Q) // dtr_cnfg):
753                Q_top.append(Q_sum[dtr_cnfg * i][index_aim])
754        if len(Q_top) != len(traj_list) // dtr_cnfg:
755            raise Exception("qtops and traj_list dont have the same length")
756
757        if kwargs.get('integer_charge', False):
758            Q_top = [round(q) for q in Q_top]
759
760        truncated_file = file[:-len(postfix)]
761
762        if "names" not in kwargs:
763            try:
764                idx = truncated_file.index('r')
765            except Exception:
766                if "names" not in kwargs:
767                    raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
768            ens_name = truncated_file[:idx]
769            rep_names.append(ens_name + '|' + truncated_file[idx:])
770        else:
771            names = kwargs.get("names")
772            rep_names = names
773        deltas.append(Q_top)
774
775    idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))]
776    deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))]
777    result = Obs(deltas, rep_names, idl=idl)
778    return result
779
780
781def qtop_projection(qtop, target=0):
782    """Returns the projection to the topological charge sector defined by target.
783
784    Parameters
785    ----------
786    path : Obs
787        Topological charge.
788    target : int
789        Specifies the topological sector to be reweighted to (default 0)
790    """
791    if qtop.reweighted:
792        raise Exception('You can not use a reweighted observable for reweighting!')
793
794    proj_qtop = []
795    for n in qtop.deltas:
796        proj_qtop.append(np.array([1 if round(qtop.value + q) == target else 0 for q in qtop.deltas[n]]))
797
798    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
799    reto.is_merged = qtop.is_merged
800    return reto
801
802
803def read_qtop_sector(path, prefix, c, target=0, **kwargs):
804    """Constructs reweighting factors to a specified topological sector.
805
806    Parameters
807    ----------
808    path : str
809        path of the measurement files
810    prefix : str
811        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
812    c : double
813        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
814    target : int
815        Specifies the topological sector to be reweighted to (default 0)
816    dtr_cnfg : int
817        (optional) parameter that specifies the number of trajectories
818        between two configs.
819        if it is not set, the distance between two measurements
820        in the file is assumed to be the distance between two configurations.
821    steps : int
822        (optional) Distance between two configurations in units of trajectories /
823         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
824    version : str
825        version string of the openQCD (sfqcd) version used to create
826        the ensemble. Default is 2.0. May also be set to sfqcd.
827    L : int
828        spatial length of the lattice in L/a.
829        HAS to be set if version != sfqcd, since openQCD does not provide
830        this in the header
831    r_start : list
832        offset of the first ensemble, making it easier to match
833        later on with other Obs
834    r_stop : list
835        last configurations that need to be read (per replicum)
836    files : list
837        specify the exact files that need to be read
838        from path, practical if e.g. only one replicum is needed
839    names : list
840        Alternative labeling for replicas/ensembles.
841        Has to have the appropriate length
842    Zeuthen_flow : bool
843        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
844        for version=='sfqcd' If False, the Wilson flow is used.
845    """
846
847    if not isinstance(target, int):
848        raise Exception("'target' has to be an integer.")
849
850    kwargs['integer_charge'] = True
851    qtop = read_qtop(path, prefix, c, **kwargs)
852
853    return qtop_projection(qtop, target=target)