# pyerrors.input.openQCD API documentation

# pyerrors.input.openQCD

View Source

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

def read_rwms(path, prefix, version='2.0', names=None, **kwargs): View Source

     85def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
     86    """Read rwms format from given folder structure. Returns a list of length nrw
     87
     88    Parameters
     89    ----------
     90    path : str
     91        path that contains the data files
     92    prefix : str
     93        all files in path that start with prefix are considered as input files.
     94        May be used together postfix to consider only special file endings.
     95        Prefix is ignored, if the keyword 'files' is used.
     96    version : str
     97        version of openQCD, default 2.0
     98    names : list
     99        list of names that is assigned to the data according according
    100        to the order in the file list. Use careful, if you do not provide file names!
    101    r_start : list
    102        list which contains the first config to be read for each replicum
    103    r_stop : list
    104        list which contains the last config to be read for each replicum
    105    r_step : int
    106        integer that defines a fixed step size between two measurements (in units of configs)
    107        If not given, r_step=1 is assumed.
    108    postfix : str
    109        postfix of the file to read, e.g. '.ms1' for openQCD-files
    110    files : list
    111        list which contains the filenames to be read. No automatic detection of
    112        files performed if given.
    113    print_err : bool
    114        Print additional information that is useful for debugging.
    115
    116    Returns
    117    -------
    118    rwms : Obs
    119        Reweighting factors read
    120    """
    121    known_oqcd_versions = ['1.4', '1.6', '2.0']
    122    if not (version in known_oqcd_versions):
    123        raise Exception('Unknown openQCD version defined!')
    124    print("Working with openQCD version " + version)
    125    if 'postfix' in kwargs:
    126        postfix = kwargs.get('postfix')
    127    else:
    128        postfix = ''
    129
    130    if 'files' in kwargs:
    131        known_files = kwargs.get('files')
    132    else:
    133        known_files = []
    134
    135    ls = _find_files(path, prefix, postfix, 'dat', known_files=known_files)
    136
    137    replica = len(ls)
    138
    139    if 'r_start' in kwargs:
    140        r_start = kwargs.get('r_start')
    141        if len(r_start) != replica:
    142            raise Exception('r_start does not match number of replicas')
    143        r_start = [o if o else None for o in r_start]
    144    else:
    145        r_start = [None] * replica
    146
    147    if 'r_stop' in kwargs:
    148        r_stop = kwargs.get('r_stop')
    149        if len(r_stop) != replica:
    150            raise Exception('r_stop does not match number of replicas')
    151    else:
    152        r_stop = [None] * replica
    153
    154    if 'r_step' in kwargs:
    155        r_step = kwargs.get('r_step')
    156    else:
    157        r_step = 1
    158
    159    print('Read reweighting factors from', prefix[:-1], ',',
    160          replica, 'replica', end='')
    161
    162    if names is None:
    163        rep_names = []
    164        for entry in ls:
    165            truncated_entry = entry
    166            suffixes = [".dat", ".rwms", ".ms1"]
    167            for suffix in suffixes:
    168                if truncated_entry.endswith(suffix):
    169                    truncated_entry = truncated_entry[0:-len(suffix)]
    170            idx = truncated_entry.index('r')
    171            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
    172    else:
    173        rep_names = names
    174
    175    rep_names = _sort_names(rep_names)
    176
    177    print_err = 0
    178    if 'print_err' in kwargs:
    179        print_err = 1
    180        print()
    181
    182    deltas = []
    183
    184    configlist = []
    185    r_start_index = []
    186    r_stop_index = []
    187
    188    for rep in range(replica):
    189        tmp_array = []
    190        with open(path + '/' + ls[rep], 'rb') as fp:
    191
    192            t = fp.read(4)  # number of reweighting factors
    193            if rep == 0:
    194                nrw = struct.unpack('i', t)[0]
    195                if version == '2.0':
    196                    nrw = int(nrw / 2)
    197                for k in range(nrw):
    198                    deltas.append([])
    199            else:
    200                if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
    201                    raise Exception('Error: different number of reweighting factors for replicum', rep)
    202
    203            for k in range(nrw):
    204                tmp_array.append([])
    205
    206            # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
    207            nfct = []
    208            if version in ['1.6', '2.0']:
    209                for i in range(nrw):
    210                    t = fp.read(4)
    211                    nfct.append(struct.unpack('i', t)[0])
    212            else:
    213                for i in range(nrw):
    214                    nfct.append(1)
    215
    216            nsrc = []
    217            for i in range(nrw):
    218                t = fp.read(4)
    219                nsrc.append(struct.unpack('i', t)[0])
    220            if version == '2.0':
    221                if not struct.unpack('i', fp.read(4))[0] == 0:
    222                    raise Exception("You are using the input for openQCD version 2.0, this is not correct.")
    223
    224            configlist.append([])
    225            while True:
    226                t = fp.read(4)
    227                if len(t) < 4:
    228                    break
    229                config_no = struct.unpack('i', t)[0]
    230                configlist[-1].append(config_no)
    231                for i in range(nrw):
    232                    if (version == '2.0'):
    233                        tmpd = _read_array_openQCD2(fp)
    234                        tmpd = _read_array_openQCD2(fp)
    235                        tmp_rw = tmpd['arr']
    236                        tmp_nfct = 1.0
    237                        for j in range(tmpd['n'][0]):
    238                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
    239                            if print_err:
    240                                print(config_no, i, j,
    241                                      np.mean(np.exp(-np.asarray(tmp_rw[j]))),
    242                                      np.std(np.exp(-np.asarray(tmp_rw[j]))))
    243                                print('Sources:',
    244                                      np.exp(-np.asarray(tmp_rw[j])))
    245                                print('Partial factor:', tmp_nfct)
    246                    elif version == '1.6' or version == '1.4':
    247                        tmp_nfct = 1.0
    248                        for j in range(nfct[i]):
    249                            t = fp.read(8 * nsrc[i])
    250                            t = fp.read(8 * nsrc[i])
    251                            tmp_rw = struct.unpack('d' * nsrc[i], t)
    252                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
    253                            if print_err:
    254                                print(config_no, i, j,
    255                                      np.mean(np.exp(-np.asarray(tmp_rw))),
    256                                      np.std(np.exp(-np.asarray(tmp_rw))))
    257                                print('Sources:', np.exp(-np.asarray(tmp_rw)))
    258                                print('Partial factor:', tmp_nfct)
    259                    tmp_array[i].append(tmp_nfct)
    260
    261            diffmeas = configlist[-1][-1] - configlist[-1][-2]
    262            configlist[-1] = [item // diffmeas for item in configlist[-1]]
    263            if configlist[-1][0] > 1 and diffmeas > 1:
    264                warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
    265                offset = configlist[-1][0] - 1
    266                configlist[-1] = [item - offset for item in configlist[-1]]
    267
    268            if r_start[rep] is None:
    269                r_start_index.append(0)
    270            else:
    271                try:
    272                    r_start_index.append(configlist[-1].index(r_start[rep]))
    273                except ValueError:
    274                    raise Exception('Config %d not in file with range [%d, %d]' % (
    275                        r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
    276
    277            if r_stop[rep] is None:
    278                r_stop_index.append(len(configlist[-1]) - 1)
    279            else:
    280                try:
    281                    r_stop_index.append(configlist[-1].index(r_stop[rep]))
    282                except ValueError:
    283                    raise Exception('Config %d not in file with range [%d, %d]' % (
    284                        r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
    285
    286            for k in range(nrw):
    287                deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
    288
    289    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
    290        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
    291    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
    292    if np.any([step != 1 for step in stepsizes]):
    293        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
    294
    295    print(',', nrw, 'reweighting factors with', nsrc, 'sources')
    296    result = []
    297    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
    298
    299    for t in range(nrw):
    300        result.append(Obs(deltas[t], rep_names, idl=idl))
    301    return result

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

###### Parameters

- path (str): path that contains the data files
- prefix (str): all files in path that start with prefix are considered as input files. May be used together postfix to consider only special file endings. Prefix is ignored, if the keyword 'files' is used.
- version (str): version of openQCD, default 2.0
- names (list): list of names that is assigned to the data according according to the order in the file list. Use careful, if you do not provide file names!
- r_start (list): list which contains the first config to be read for each replicum
- r_stop (list): list which contains the last config to be read for each replicum
- r_step (int): integer that defines a fixed step size between two measurements (in units of configs) If not given, r_step=1 is assumed.
- postfix (str): postfix of the file to read, e.g. '.ms1' for openQCD-files
- files (list): list which contains the filenames to be read. No automatic detection of files performed if given.
- print_err (bool): Print additional information that is useful for debugging.

###### Returns

- rwms (Obs): Reweighting factors read

def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs): View Source

    304def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
    305    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
    306
    307    It is assumed that all boundary effects have
    308    sufficiently decayed at x0=xmin.
    309    The data around the zero crossing of t^2<E> - 0.3
    310    is fitted with a linear function
    311    from which the exact root is extracted.
    312
    313    It is assumed that one measurement is performed for each config.
    314    If this is not the case, the resulting idl, as well as the handling
    315    of r_start, r_stop and r_step is wrong and the user has to correct
    316    this in the resulting observable.
    317
    318    Parameters
    319    ----------
    320    path : str
    321        Path to .ms.dat files
    322    prefix : str
    323        Ensemble prefix
    324    dtr_read : int
    325        Determines how many trajectories should be skipped
    326        when reading the ms.dat files.
    327        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
    328    xmin : int
    329        First timeslice where the boundary
    330        effects have sufficiently decayed.
    331    spatial_extent : int
    332        spatial extent of the lattice, required for normalization.
    333    fit_range : int
    334        Number of data points left and right of the zero
    335        crossing to be included in the linear fit. (Default: 5)
    336    r_start : list
    337        list which contains the first config to be read for each replicum.
    338    r_stop : list
    339        list which contains the last config to be read for each replicum.
    340    r_step : int
    341        integer that defines a fixed step size between two measurements (in units of configs)
    342        If not given, r_step=1 is assumed.
    343    plaquette : bool
    344        If true extract the plaquette estimate of t0 instead.
    345    names : list
    346        list of names that is assigned to the data according according
    347        to the order in the file list. Use careful, if you do not provide file names!
    348    files : list
    349        list which contains the filenames to be read. No automatic detection of
    350        files performed if given.
    351    plot_fit : bool
    352        If true, the fit for the extraction of t0 is shown together with the data.
    353    assume_thermalization : bool
    354        If True: If the first record divided by the distance between two measurements is larger than
    355        1, it is assumed that this is due to thermalization and the first measurement belongs
    356        to the first config (default).
    357        If False: The config numbers are assumed to be traj_number // difference
    358
    359    Returns
    360    -------
    361    t0 : Obs
    362        Extracted t0
    363    """
    364
    365    if 'files' in kwargs:
    366        known_files = kwargs.get('files')
    367    else:
    368        known_files = []
    369
    370    ls = _find_files(path, prefix, 'ms', 'dat', known_files=known_files)
    371
    372    replica = len(ls)
    373
    374    if 'r_start' in kwargs:
    375        r_start = kwargs.get('r_start')
    376        if len(r_start) != replica:
    377            raise Exception('r_start does not match number of replicas')
    378        r_start = [o if o else None for o in r_start]
    379    else:
    380        r_start = [None] * replica
    381
    382    if 'r_stop' in kwargs:
    383        r_stop = kwargs.get('r_stop')
    384        if len(r_stop) != replica:
    385            raise Exception('r_stop does not match number of replicas')
    386    else:
    387        r_stop = [None] * replica
    388
    389    if 'r_step' in kwargs:
    390        r_step = kwargs.get('r_step')
    391    else:
    392        r_step = 1
    393
    394    print('Extract t0 from', prefix, ',', replica, 'replica')
    395
    396    if 'names' in kwargs:
    397        rep_names = kwargs.get('names')
    398    else:
    399        rep_names = []
    400        for entry in ls:
    401            truncated_entry = entry.split('.')[0]
    402            idx = truncated_entry.index('r')
    403            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
    404
    405    Ysum = []
    406
    407    configlist = []
    408    r_start_index = []
    409    r_stop_index = []
    410
    411    for rep in range(replica):
    412
    413        with open(path + '/' + ls[rep], 'rb') as fp:
    414            t = fp.read(12)
    415            header = struct.unpack('iii', t)
    416            if rep == 0:
    417                dn = header[0]
    418                nn = header[1]
    419                tmax = header[2]
    420            elif dn != header[0] or nn != header[1] or tmax != header[2]:
    421                raise Exception('Replica parameters do not match.')
    422
    423            t = fp.read(8)
    424            if rep == 0:
    425                eps = struct.unpack('d', t)[0]
    426                print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
    427            elif eps != struct.unpack('d', t)[0]:
    428                raise Exception('Values for eps do not match among replica.')
    429
    430            Ysl = []
    431
    432            configlist.append([])
    433            while True:
    434                t = fp.read(4)
    435                if (len(t) < 4):
    436                    break
    437                nc = struct.unpack('i', t)[0]
    438                configlist[-1].append(nc)
    439
    440                t = fp.read(8 * tmax * (nn + 1))
    441                if kwargs.get('plaquette'):
    442                    if nc % dtr_read == 0:
    443                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
    444                t = fp.read(8 * tmax * (nn + 1))
    445                if not kwargs.get('plaquette'):
    446                    if nc % dtr_read == 0:
    447                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
    448                t = fp.read(8 * tmax * (nn + 1))
    449
    450        Ysum.append([])
    451        for i, item in enumerate(Ysl):
    452            Ysum[-1].append([np.mean(item[current + xmin:
    453                             current + tmax - xmin])
    454                            for current in range(0, len(item), tmax)])
    455
    456        diffmeas = configlist[-1][-1] - configlist[-1][-2]
    457        configlist[-1] = [item // diffmeas for item in configlist[-1]]
    458        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
    459            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
    460            offset = configlist[-1][0] - 1
    461            configlist[-1] = [item - offset for item in configlist[-1]]
    462
    463        if r_start[rep] is None:
    464            r_start_index.append(0)
    465        else:
    466            try:
    467                r_start_index.append(configlist[-1].index(r_start[rep]))
    468            except ValueError:
    469                raise Exception('Config %d not in file with range [%d, %d]' % (
    470                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
    471
    472        if r_stop[rep] is None:
    473            r_stop_index.append(len(configlist[-1]) - 1)
    474        else:
    475            try:
    476                r_stop_index.append(configlist[-1].index(r_stop[rep]))
    477            except ValueError:
    478                raise Exception('Config %d not in file with range [%d, %d]' % (
    479                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
    480
    481    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
    482        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
    483    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
    484    if np.any([step != 1 for step in stepsizes]):
    485        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
    486
    487    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
    488    t2E_dict = {}
    489    for n in range(nn + 1):
    490        samples = []
    491        for nrep, rep in enumerate(Ysum):
    492            samples.append([])
    493            for cnfg in rep:
    494                samples[-1].append(cnfg[n])
    495            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
    496        new_obs = Obs(samples, rep_names, idl=idl)
    497        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
    498
    499    zero_crossing = np.argmax(np.array(
    500        [o.value for o in t2E_dict.values()]) > 0.0)
    501
    502    x = list(t2E_dict.keys())[zero_crossing - fit_range:
    503                              zero_crossing + fit_range]
    504    y = list(t2E_dict.values())[zero_crossing - fit_range:
    505                                zero_crossing + fit_range]
    506    [o.gamma_method() for o in y]
    507
    508    fit_result = fit_lin(x, y)
    509
    510    if kwargs.get('plot_fit'):
    511        plt.figure()
    512        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
    513        ax0 = plt.subplot(gs[0])
    514        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
    515        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
    516        [o.gamma_method() for o in ymore]
    517        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
    518        xplot = np.linspace(np.min(x), np.max(x))
    519        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
    520        [yi.gamma_method() for yi in yplot]
    521        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
    522        retval = (-fit_result[0] / fit_result[1])
    523        retval.gamma_method()
    524        ylim = ax0.get_ylim()
    525        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
    526        ax0.set_ylim(ylim)
    527        ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
    528        xlim = ax0.get_xlim()
    529
    530        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
    531        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
    532        ax1 = plt.subplot(gs[1])
    533        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
    534        ax1.tick_params(direction='out')
    535        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
    536        ax1.axhline(y=0.0, ls='--', color='k')
    537        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
    538        ax1.set_xlim(xlim)
    539        ax1.set_ylabel('Residuals')
    540        ax1.set_xlabel(r'$t/a^2$')
    541
    542        plt.draw()
    543    return -fit_result[0] / fit_result[1]

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

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

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

###### Parameters

- path (str): Path to .ms.dat files
- prefix (str): Ensemble prefix
- dtr_read (int): Determines how many trajectories should be skipped when reading the ms.dat files. Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
- xmin (int): First timeslice where the boundary effects have sufficiently decayed.
- spatial_extent (int): spatial extent of the lattice, required for normalization.
- fit_range (int): Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5)
- r_start (list): list which contains the first config to be read for each replicum.
- r_stop (list): list which contains the last config to be read for each replicum.
- r_step (int): integer that defines a fixed step size between two measurements (in units of configs) If not given, r_step=1 is assumed.
- plaquette (bool): If true extract the plaquette estimate of t0 instead.
- names (list): list of names that is assigned to the data according according to the order in the file list. Use careful, if you do not provide file names!
- files (list): list which contains the filenames to be read. No automatic detection of files performed if given.
- plot_fit (bool): If true, the fit for the extraction of t0 is shown together with the data.
- assume_thermalization (bool): If True: If the first record divided by the distance between two measurements is larger than 1, it is assumed that this is due to thermalization and the first measurement belongs to the first config (default). If False: The config numbers are assumed to be traj_number // difference

###### Returns

- t0 (Obs): Extracted t0

def read_qtop(path, prefix, c, dtr_cnfg=1, version='openQCD', **kwargs): View Source

    591def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
    592    """Read the topologial charge based on openQCD gradient flow measurements.
    593
    594    Parameters
    595    ----------
    596    path : str
    597        path of the measurement files
    598    prefix : str
    599        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
    600        Ignored if file names are passed explicitly via keyword files.
    601    c : double
    602        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
    603    dtr_cnfg : int
    604        (optional) parameter that specifies the number of measurements
    605        between two configs.
    606        If it is not set, the distance between two measurements
    607        in the file is assumed to be the distance between two configurations.
    608    steps : int
    609        (optional) Distance between two configurations in units of trajectories /
    610         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    611    version : str
    612        Either openQCD or sfqcd, depending on the data.
    613    L : int
    614        spatial length of the lattice in L/a.
    615        HAS to be set if version != sfqcd, since openQCD does not provide
    616        this in the header
    617    r_start : list
    618        list which contains the first config to be read for each replicum.
    619    r_stop : list
    620        list which contains the last config to be read for each replicum.
    621    files : list
    622        specify the exact files that need to be read
    623        from path, practical if e.g. only one replicum is needed
    624    postfix : str
    625        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
    626    names : list
    627        Alternative labeling for replicas/ensembles.
    628        Has to have the appropriate length.
    629    Zeuthen_flow : bool
    630        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
    631        for version=='sfqcd' If False, the Wilson flow is used.
    632    integer_charge : bool
    633        If True, the charge is rounded towards the nearest integer on each config.
    634
    635    Returns
    636    -------
    637    result : Obs
    638        Read topological charge
    639    """
    640
    641    return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)

Read the topologial charge based on openQCD gradient flow measurements.

###### Parameters

- path (str): path of the measurement files
- prefix (str): prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. Ignored if file names are passed explicitly via keyword files.
- c (double): Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
- dtr_cnfg (int): (optional) parameter that specifies the number of measurements between two configs. If it is not set, the distance between two measurements in the file is assumed to be the distance between two configurations.
- steps (int): (optional) Distance between two configurations in units of trajectories / cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- version (str): Either openQCD or sfqcd, depending on the data.
- L (int): spatial length of the lattice in L/a. HAS to be set if version != sfqcd, since openQCD does not provide this in the header
- r_start (list): list which contains the first config to be read for each replicum.
- r_stop (list): list which contains the last config to be read for each replicum.
- files (list): specify the exact files that need to be read from path, practical if e.g. only one replicum is needed
- postfix (str): postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
- names (list): Alternative labeling for replicas/ensembles. Has to have the appropriate length.
- Zeuthen_flow (bool): (optional) If True, the Zeuthen flow is used for Qtop. Only possible for version=='sfqcd' If False, the Wilson flow is used.
- integer_charge (bool): If True, the charge is rounded towards the nearest integer on each config.

###### Returns

- result (Obs): Read topological charge

def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs): View Source

    644def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
    645    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
    646
    647    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.
    648
    649    Parameters
    650    ----------
    651    path : str
    652        path of the measurement files
    653    prefix : str
    654        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
    655        Ignored if file names are passed explicitly via keyword files.
    656    c : double
    657        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
    658    dtr_cnfg : int
    659        (optional) parameter that specifies the number of measurements
    660        between two configs.
    661        If it is not set, the distance between two measurements
    662        in the file is assumed to be the distance between two configurations.
    663    steps : int
    664        (optional) Distance between two configurations in units of trajectories /
    665         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    666    r_start : list
    667        list which contains the first config to be read for each replicum.
    668    r_stop : list
    669        list which contains the last config to be read for each replicum.
    670    files : list
    671        specify the exact files that need to be read
    672        from path, practical if e.g. only one replicum is needed
    673    names : list
    674        Alternative labeling for replicas/ensembles.
    675        Has to have the appropriate length.
    676    postfix : str
    677        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
    678    Zeuthen_flow : bool
    679        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
    680    """
    681
    682    if c != 0.3:
    683        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
    684
    685    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)
    686    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)
    687    L = plaq.tag["L"]
    688    T = plaq.tag["T"]
    689
    690    if T != L:
    691        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
    692
    693    if Zeuthen_flow is not True:
    694        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
    695
    696    t = (c * L) ** 2 / 8
    697
    698    normdict = {4: 0.012341170468270,
    699                6: 0.010162691462430,
    700                8: 0.009031614807931,
    701                10: 0.008744966371393,
    702                12: 0.008650917856809,
    703                14: 8.611154391267955E-03,
    704                16: 0.008591758449508,
    705                20: 0.008575359627103,
    706                24: 0.008569387847540,
    707                28: 8.566803713382559E-03,
    708                32: 0.008565541650006,
    709                40: 8.564480684962046E-03,
    710                48: 8.564098025073460E-03,
    711                64: 8.563853943383087E-03}
    712
    713    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]

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

Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step.

###### Parameters

- path (str): path of the measurement files
- prefix (str): prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat. Ignored if file names are passed explicitly via keyword files.
- c (double): Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
- dtr_cnfg (int): (optional) parameter that specifies the number of measurements between two configs. If it is not set, the distance between two measurements in the file is assumed to be the distance between two configurations.
- steps (int): (optional) Distance between two configurations in units of trajectories / cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- r_start (list): list which contains the first config to be read for each replicum.
- r_stop (list): list which contains the last config to be read for each replicum.
- files (list): specify the exact files that need to be read from path, practical if e.g. only one replicum is needed
- names (list): Alternative labeling for replicas/ensembles. Has to have the appropriate length.
- postfix (str): postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
- Zeuthen_flow (bool): (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.

def qtop_projection(qtop, target=0): View Source

     988def qtop_projection(qtop, target=0):
     989    """Returns the projection to the topological charge sector defined by target.
     990
     991    Parameters
     992    ----------
     993    path : Obs
     994        Topological charge.
     995    target : int
     996        Specifies the topological sector to be reweighted to (default 0)
     997
     998    Returns
     999    -------
    1000    reto : Obs
    1001        projection to the topological charge sector defined by target
    1002    """
    1003    if qtop.reweighted:
    1004        raise Exception('You can not use a reweighted observable for reweighting!')
    1005
    1006    proj_qtop = []
    1007    for n in qtop.deltas:
    1008        proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
    1009
    1010    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
    1011    return reto

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

###### Parameters

- path (Obs): Topological charge.
- target (int): Specifies the topological sector to be reweighted to (default 0)

###### Returns

- reto (Obs): projection to the topological charge sector defined by target

def read_qtop_sector(path, prefix, c, target=0, **kwargs): View Source

    1014def read_qtop_sector(path, prefix, c, target=0, **kwargs):
    1015    """Constructs reweighting factors to a specified topological sector.
    1016
    1017    Parameters
    1018    ----------
    1019    path : str
    1020        path of the measurement files
    1021    prefix : str
    1022        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
    1023    c : double
    1024        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
    1025    target : int
    1026        Specifies the topological sector to be reweighted to (default 0)
    1027    dtr_cnfg : int
    1028        (optional) parameter that specifies the number of trajectories
    1029        between two configs.
    1030        if it is not set, the distance between two measurements
    1031        in the file is assumed to be the distance between two configurations.
    1032    steps : int
    1033        (optional) Distance between two configurations in units of trajectories /
    1034         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    1035    version : str
    1036        version string of the openQCD (sfqcd) version used to create
    1037        the ensemble. Default is 2.0. May also be set to sfqcd.
    1038    L : int
    1039        spatial length of the lattice in L/a.
    1040        HAS to be set if version != sfqcd, since openQCD does not provide
    1041        this in the header
    1042    r_start : list
    1043        offset of the first ensemble, making it easier to match
    1044        later on with other Obs
    1045    r_stop : list
    1046        last configurations that need to be read (per replicum)
    1047    files : list
    1048        specify the exact files that need to be read
    1049        from path, practical if e.g. only one replicum is needed
    1050    names : list
    1051        Alternative labeling for replicas/ensembles.
    1052        Has to have the appropriate length
    1053    Zeuthen_flow : bool
    1054        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
    1055        for version=='sfqcd' If False, the Wilson flow is used.
    1056
    1057    Returns
    1058    -------
    1059    reto : Obs
    1060        projection to the topological charge sector defined by target
    1061    """
    1062
    1063    if not isinstance(target, int):
    1064        raise Exception("'target' has to be an integer.")
    1065
    1066    kwargs['integer_charge'] = True
    1067    qtop = read_qtop(path, prefix, c, **kwargs)
    1068
    1069    return qtop_projection(qtop, target=target)

Constructs reweighting factors to a specified topological sector.

###### Parameters

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

###### Returns

- reto (Obs): projection to the topological charge sector defined by target

def read_ms5_xsf(path, prefix, qc, corr, sep='r', **kwargs): View Source

    1072def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
    1073    """
    1074    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
    1075
    1076    Parameters
    1077    ----------
    1078    path : str
    1079        The directory to search for the files in.
    1080    prefix : str
    1081        The prefix to match the files against.
    1082    qc : str
    1083        The quark combination extension to match the files against.
    1084    corr : str
    1085        The correlator to extract data for.
    1086    sep : str, optional
    1087        The separator to use when parsing the replika names.
    1088    **kwargs
    1089        Additional keyword arguments. The following keyword arguments are recognized:
    1090
    1091        - names (List[str]): A list of names to use for the replicas.
    1092
    1093    Returns
    1094    -------
    1095    Corr
    1096        A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
    1097    or
    1098    CObs
    1099        A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
    1100
    1101
    1102    Raises
    1103    ------
    1104    FileNotFoundError
    1105        If no files matching the specified prefix and quark combination extension are found in the specified directory.
    1106    IOError
    1107        If there is an error reading a file.
    1108    struct.error
    1109        If there is an error unpacking binary data.
    1110    """
    1111
    1112    # found = []
    1113    files = []
    1114    names = []
    1115
    1116    # test if the input is correct
    1117    if qc not in ['dd', 'ud', 'du', 'uu']:
    1118        raise Exception("Unknown quark conbination!")
    1119
    1120    if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
    1121        raise Exception("Unknown correlator!")
    1122
    1123    if "files" in kwargs:
    1124        known_files = kwargs.get("files")
    1125    else:
    1126        known_files = []
    1127    files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
    1128
    1129    if "names" in kwargs:
    1130        names = kwargs.get("names")
    1131    else:
    1132        for f in files:
    1133            if not sep == "":
    1134                se = f.split(".")[0]
    1135                for s in f.split(".")[1:-2]:
    1136                    se += "." + s
    1137                names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
    1138            else:
    1139                names.append(prefix)
    1140
    1141    names = sorted(names)
    1142    files = sorted(files)
    1143
    1144    cnfgs = []
    1145    realsamples = []
    1146    imagsamples = []
    1147    repnum = 0
    1148    for file in files:
    1149        with open(path + "/" + file, "rb") as fp:
    1150
    1151            t = fp.read(8)
    1152            kappa = struct.unpack('d', t)[0]
    1153            t = fp.read(8)
    1154            csw = struct.unpack('d', t)[0]
    1155            t = fp.read(8)
    1156            dF = struct.unpack('d', t)[0]
    1157            t = fp.read(8)
    1158            zF = struct.unpack('d', t)[0]
    1159
    1160            t = fp.read(4)
    1161            tmax = struct.unpack('i', t)[0]
    1162            t = fp.read(4)
    1163            bnd = struct.unpack('i', t)[0]
    1164
    1165            placesBI = ["gS", "gP",
    1166                        "gA", "gV",
    1167                        "gVt", "lA",
    1168                        "lV", "lVt",
    1169                        "lT", "lTt"]
    1170            placesBB = ["g1", "l1"]
    1171
    1172            # the chunks have the following structure:
    1173            # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
    1174
    1175            chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
    1176            packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
    1177            cnfgs.append([])
    1178            realsamples.append([])
    1179            imagsamples.append([])
    1180            for t in range(tmax):
    1181                realsamples[repnum].append([])
    1182                imagsamples[repnum].append([])
    1183
    1184            while True:
    1185                cnfgt = fp.read(chunksize)
    1186                if not cnfgt:
    1187                    break
    1188                asascii = struct.unpack(packstr, cnfgt)
    1189                cnfg = asascii[0]
    1190                cnfgs[repnum].append(cnfg)
    1191
    1192                if corr not in placesBB:
    1193                    tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
    1194                else:
    1195                    tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
    1196
    1197                corrres = [[], []]
    1198                for i in range(len(tmpcorr)):
    1199                    corrres[i % 2].append(tmpcorr[i])
    1200                for t in range(int(len(tmpcorr) / 2)):
    1201                    realsamples[repnum][t].append(corrres[0][t])
    1202                for t in range(int(len(tmpcorr) / 2)):
    1203                    imagsamples[repnum][t].append(corrres[1][t])
    1204        repnum += 1
    1205
    1206    s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
    1207    for rep in range(1, repnum):
    1208        s += ", " + str(len(realsamples[rep][t]))
    1209    s += " samples"
    1210    print(s)
    1211    print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
    1212
    1213    # we have the data now... but we need to re format the whole thing and put it into Corr objects.
    1214
    1215    compObs = []
    1216
    1217    for t in range(int(len(tmpcorr) / 2)):
    1218        compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
    1219                            Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
    1220
    1221    if len(compObs) == 1:
    1222        return compObs[0]
    1223    else:
    1224        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

- path (str): The directory to search for the files in.
- prefix (str): The prefix to match the files against.
- qc (str): The quark combination extension to match the files against.
- corr (str): The correlator to extract data for.
- sep (str, optional): The separator to use when parsing the replika names.

- **kwargs: Additional keyword arguments. The following keyword arguments are recognized:

  - names (List[str]): A list of names to use for the replicas.

###### Returns

- Corr: A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
- or
- CObs: A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.

###### Raises

- FileNotFoundError: If no files matching the specified prefix and quark combination extension are found in the specified directory.
- IOError: If there is an error reading a file.
- struct.error: If there is an error unpacking binary data.