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

    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a Corr object containing the data.

    + +
    Parameters
    + + + +
    Returns
    + + + +
    Raises
    + + +
    + +