# pyerrors.input.openQCD

View Source

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

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

     14def read_rwms(path, prefix, version='2.0', names=None, **kwargs):
     15    """Read rwms format from given folder structure. Returns a list of length nrw
     16
     17    Parameters
     18    ----------
     19    path : str
     20        path that contains the data files
     21    prefix : str
     22        all files in path that start with prefix are considered as input files.
     23        May be used together postfix to consider only special file endings.
     24        Prefix is ignored, if the keyword 'files' is used.
     25    version : str
     26        version of openQCD, default 2.0
     27    names : list
     28        list of names that is assigned to the data according according
     29        to the order in the file list. Use careful, if you do not provide file names!
     30    r_start : list
     31        list which contains the first config to be read for each replicum
     32    r_stop : list
     33        list which contains the last config to be read for each replicum
     34    r_step : int
     35        integer that defines a fixed step size between two measurements (in units of configs)
     36        If not given, r_step=1 is assumed.
     37    postfix : str
     38        postfix of the file to read, e.g. '.ms1' for openQCD-files
     39    files : list
     40        list which contains the filenames to be read. No automatic detection of
     41        files performed if given.
     42    print_err : bool
     43        Print additional information that is useful for debugging.
     44
     45    Returns
     46    -------
     47    rwms : Obs
     48        Reweighting factors read
     49    """
     50    known_oqcd_versions = ['1.4', '1.6', '2.0']
     51    if not (version in known_oqcd_versions):
     52        raise Exception('Unknown openQCD version defined!')
     53    print("Working with openQCD version " + version)
     54    if 'postfix' in kwargs:
     55        postfix = kwargs.get('postfix')
     56    else:
     57        postfix = ''
     58
     59    if 'files' in kwargs:
     60        known_files = kwargs.get('files')
     61    else:
     62        known_files = []
     63
     64    ls = _find_files(path, prefix, postfix, 'dat', known_files=known_files)
     65
     66    replica = len(ls)
     67
     68    if 'r_start' in kwargs:
     69        r_start = kwargs.get('r_start')
     70        if len(r_start) != replica:
     71            raise Exception('r_start does not match number of replicas')
     72        r_start = [o if o else None for o in r_start]
     73    else:
     74        r_start = [None] * replica
     75
     76    if 'r_stop' in kwargs:
     77        r_stop = kwargs.get('r_stop')
     78        if len(r_stop) != replica:
     79            raise Exception('r_stop does not match number of replicas')
     80    else:
     81        r_stop = [None] * replica
     82
     83    if 'r_step' in kwargs:
     84        r_step = kwargs.get('r_step')
     85    else:
     86        r_step = 1
     87
     88    print('Read reweighting factors from', prefix[:-1], ',',
     89          replica, 'replica', end='')
     90
     91    if names is None:
     92        rep_names = []
     93        for entry in ls:
     94            truncated_entry = entry
     95            suffixes = [".dat", ".rwms", ".ms1"]
     96            for suffix in suffixes:
     97                if truncated_entry.endswith(suffix):
     98                    truncated_entry = truncated_entry[0:-len(suffix)]
     99            idx = truncated_entry.index('r')
    100            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
    101    else:
    102        rep_names = names
    103
    104    rep_names = sort_names(rep_names)
    105
    106    print_err = 0
    107    if 'print_err' in kwargs:
    108        print_err = 1
    109        print()
    110
    111    deltas = []
    112
    113    configlist = []
    114    r_start_index = []
    115    r_stop_index = []
    116
    117    for rep in range(replica):
    118        tmp_array = []
    119        with open(path + '/' + ls[rep], 'rb') as fp:
    120
    121            t = fp.read(4)  # number of reweighting factors
    122            if rep == 0:
    123                nrw = struct.unpack('i', t)[0]
    124                if version == '2.0':
    125                    nrw = int(nrw / 2)
    126                for k in range(nrw):
    127                    deltas.append([])
    128            else:
    129                if ((nrw != struct.unpack('i', t)[0] and (not version == '2.0')) or (nrw != struct.unpack('i', t)[0] / 2 and version == '2.0')):
    130                    raise Exception('Error: different number of reweighting factors for replicum', rep)
    131
    132            for k in range(nrw):
    133                tmp_array.append([])
    134
    135            # This block is necessary for openQCD1.6 and openQCD2.0 ms1 files
    136            nfct = []
    137            if version in ['1.6', '2.0']:
    138                for i in range(nrw):
    139                    t = fp.read(4)
    140                    nfct.append(struct.unpack('i', t)[0])
    141            else:
    142                for i in range(nrw):
    143                    nfct.append(1)
    144
    145            nsrc = []
    146            for i in range(nrw):
    147                t = fp.read(4)
    148                nsrc.append(struct.unpack('i', t)[0])
    149            if version == '2.0':
    150                if not struct.unpack('i', fp.read(4))[0] == 0:
    151                    raise Exception("You are using the input for openQCD version 2.0, this is not correct.")
    152
    153            configlist.append([])
    154            while True:
    155                t = fp.read(4)
    156                if len(t) < 4:
    157                    break
    158                config_no = struct.unpack('i', t)[0]
    159                configlist[-1].append(config_no)
    160                for i in range(nrw):
    161                    if (version == '2.0'):
    162                        tmpd = _read_array_openQCD2(fp)
    163                        tmpd = _read_array_openQCD2(fp)
    164                        tmp_rw = tmpd['arr']
    165                        tmp_nfct = 1.0
    166                        for j in range(tmpd['n'][0]):
    167                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw[j])))
    168                            if print_err:
    169                                print(config_no, i, j,
    170                                      np.mean(np.exp(-np.asarray(tmp_rw[j]))),
    171                                      np.std(np.exp(-np.asarray(tmp_rw[j]))))
    172                                print('Sources:',
    173                                      np.exp(-np.asarray(tmp_rw[j])))
    174                                print('Partial factor:', tmp_nfct)
    175                    elif version == '1.6' or version == '1.4':
    176                        tmp_nfct = 1.0
    177                        for j in range(nfct[i]):
    178                            t = fp.read(8 * nsrc[i])
    179                            t = fp.read(8 * nsrc[i])
    180                            tmp_rw = struct.unpack('d' * nsrc[i], t)
    181                            tmp_nfct *= np.mean(np.exp(-np.asarray(tmp_rw)))
    182                            if print_err:
    183                                print(config_no, i, j,
    184                                      np.mean(np.exp(-np.asarray(tmp_rw))),
    185                                      np.std(np.exp(-np.asarray(tmp_rw))))
    186                                print('Sources:', np.exp(-np.asarray(tmp_rw)))
    187                                print('Partial factor:', tmp_nfct)
    188                    tmp_array[i].append(tmp_nfct)
    189
    190            diffmeas = configlist[-1][-1] - configlist[-1][-2]
    191            configlist[-1] = [item // diffmeas for item in configlist[-1]]
    192            if configlist[-1][0] > 1 and diffmeas > 1:
    193                warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
    194                offset = configlist[-1][0] - 1
    195                configlist[-1] = [item - offset for item in configlist[-1]]
    196
    197            if r_start[rep] is None:
    198                r_start_index.append(0)
    199            else:
    200                try:
    201                    r_start_index.append(configlist[-1].index(r_start[rep]))
    202                except ValueError:
    203                    raise Exception('Config %d not in file with range [%d, %d]' % (
    204                        r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
    205
    206            if r_stop[rep] is None:
    207                r_stop_index.append(len(configlist[-1]) - 1)
    208            else:
    209                try:
    210                    r_stop_index.append(configlist[-1].index(r_stop[rep]))
    211                except ValueError:
    212                    raise Exception('Config %d not in file with range [%d, %d]' % (
    213                        r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
    214
    215            for k in range(nrw):
    216                deltas[k].append(tmp_array[k][r_start_index[rep]:r_stop_index[rep] + 1][::r_step])
    217
    218    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
    219        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
    220    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
    221    if np.any([step != 1 for step in stepsizes]):
    222        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
    223
    224    print(',', nrw, 'reweighting factors with', nsrc, 'sources')
    225    result = []
    226    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
    227
    228    for t in range(nrw):
    229        result.append(Obs(deltas[t], rep_names, idl=idl))
    230    return result

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

###### Parameters

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

###### Returns

- rwms (Obs): Reweighting factors read

def extract_t0( path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfix='ms', c=0.3, **kwargs): View Source

    426def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfix='ms', c=0.3, **kwargs):
    427    """Extract t0/a^2 from given .ms.dat files. Returns t0 as Obs.
    428
    429    It is assumed that all boundary effects have
    430    sufficiently decayed at x0=xmin.
    431    The data around the zero crossing of t^2<E> - c (where c=0.3 by default)
    432    is fitted with a linear function
    433    from which the exact root is extracted.
    434
    435    It is assumed that one measurement is performed for each config.
    436    If this is not the case, the resulting idl, as well as the handling
    437    of r_start, r_stop and r_step is wrong and the user has to correct
    438    this in the resulting observable.
    439
    440    Parameters
    441    ----------
    442    path : str
    443        Path to .ms.dat files
    444    prefix : str
    445        Ensemble prefix
    446    dtr_read : int
    447        Determines how many trajectories should be skipped
    448        when reading the ms.dat files.
    449        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
    450    xmin : int
    451        First timeslice where the boundary
    452        effects have sufficiently decayed.
    453    spatial_extent : int
    454        spatial extent of the lattice, required for normalization.
    455    fit_range : int
    456        Number of data points left and right of the zero
    457        crossing to be included in the linear fit. (Default: 5)
    458    postfix : str
    459        Postfix of measurement file (Default: ms)
    460    c: float
    461        Constant that defines the flow scale. Default 0.3 for t_0, choose 2./3 for t_1.
    462    r_start : list
    463        list which contains the first config to be read for each replicum.
    464    r_stop : list
    465        list which contains the last config to be read for each replicum.
    466    r_step : int
    467        integer that defines a fixed step size between two measurements (in units of configs)
    468        If not given, r_step=1 is assumed.
    469    plaquette : bool
    470        If true extract the plaquette estimate of t0 instead.
    471    names : list
    472        list of names that is assigned to the data according according
    473        to the order in the file list. Use careful, if you do not provide file names!
    474    files : list
    475        list which contains the filenames to be read. No automatic detection of
    476        files performed if given.
    477    plot_fit : bool
    478        If true, the fit for the extraction of t0 is shown together with the data.
    479    assume_thermalization : bool
    480        If True: If the first record divided by the distance between two measurements is larger than
    481        1, it is assumed that this is due to thermalization and the first measurement belongs
    482        to the first config (default).
    483        If False: The config numbers are assumed to be traj_number // difference
    484
    485    Returns
    486    -------
    487    t0 : Obs
    488        Extracted t0
    489    """
    490
    491    E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, postfix, **kwargs)
    492    t2E_dict = {}
    493    for t in sorted(E_dict.keys()):
    494        t2E_dict[t] = t ** 2 * E_dict[t] - c
    495
    496    return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit'))

Extract t0/a^2 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 - c (where c=0.3 by default) 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)
- postfix (str): Postfix of measurement file (Default: ms)
- c (float): Constant that defines the flow scale. Default 0.3 for t_0, choose 2./3 for t_1.
- 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 extract_w0( path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfix='ms', c=0.3, **kwargs): View Source

    499def extract_w0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfix='ms', c=0.3, **kwargs):
    500    """Extract w0/a from given .ms.dat files. Returns w0 as Obs.
    501
    502    It is assumed that all boundary effects have
    503    sufficiently decayed at x0=xmin.
    504    The data around the zero crossing of t d(t^2<E>)/dt -  (where c=0.3 by default)
    505    is fitted with a linear function
    506    from which the exact root is extracted.
    507
    508    It is assumed that one measurement is performed for each config.
    509    If this is not the case, the resulting idl, as well as the handling
    510    of r_start, r_stop and r_step is wrong and the user has to correct
    511    this in the resulting observable.
    512
    513    Parameters
    514    ----------
    515    path : str
    516        Path to .ms.dat files
    517    prefix : str
    518        Ensemble prefix
    519    dtr_read : int
    520        Determines how many trajectories should be skipped
    521        when reading the ms.dat files.
    522        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
    523    xmin : int
    524        First timeslice where the boundary
    525        effects have sufficiently decayed.
    526    spatial_extent : int
    527        spatial extent of the lattice, required for normalization.
    528    fit_range : int
    529        Number of data points left and right of the zero
    530        crossing to be included in the linear fit. (Default: 5)
    531    postfix : str
    532        Postfix of measurement file (Default: ms)
    533    c: float
    534        Constant that defines the flow scale. Default 0.3 for w_0, choose 2./3 for w_1.
    535    r_start : list
    536        list which contains the first config to be read for each replicum.
    537    r_stop : list
    538        list which contains the last config to be read for each replicum.
    539    r_step : int
    540        integer that defines a fixed step size between two measurements (in units of configs)
    541        If not given, r_step=1 is assumed.
    542    plaquette : bool
    543        If true extract the plaquette estimate of w0 instead.
    544    names : list
    545        list of names that is assigned to the data according according
    546        to the order in the file list. Use careful, if you do not provide file names!
    547    files : list
    548        list which contains the filenames to be read. No automatic detection of
    549        files performed if given.
    550    plot_fit : bool
    551        If true, the fit for the extraction of w0 is shown together with the data.
    552    assume_thermalization : bool
    553        If True: If the first record divided by the distance between two measurements is larger than
    554        1, it is assumed that this is due to thermalization and the first measurement belongs
    555        to the first config (default).
    556        If False: The config numbers are assumed to be traj_number // difference
    557
    558    Returns
    559    -------
    560    w0 : Obs
    561        Extracted w0
    562    """
    563
    564    E_dict = _extract_flowed_energy_density(path, prefix, dtr_read, xmin, spatial_extent, postfix, **kwargs)
    565
    566    ftimes = sorted(E_dict.keys())
    567
    568    t2E_dict = {}
    569    for t in ftimes:
    570        t2E_dict[t] = t ** 2 * E_dict[t]
    571
    572    tdtt2E_dict = {}
    573    tdtt2E_dict[ftimes[0]] = ftimes[0] * (t2E_dict[ftimes[1]] - t2E_dict[ftimes[0]]) / (ftimes[1] - ftimes[0]) - c
    574    for i in range(1, len(ftimes) - 1):
    575        tdtt2E_dict[ftimes[i]] = ftimes[i] * (t2E_dict[ftimes[i + 1]] - t2E_dict[ftimes[i - 1]]) / (ftimes[i + 1] - ftimes[i - 1]) - c
    576    tdtt2E_dict[ftimes[-1]] = ftimes[-1] * (t2E_dict[ftimes[-1]] - t2E_dict[ftimes[-2]]) / (ftimes[-1] - ftimes[-2]) - c
    577
    578    return np.sqrt(fit_t0(tdtt2E_dict, fit_range, plot_fit=kwargs.get('plot_fit'), observable='w0'))

Extract w0/a from given .ms.dat files. Returns w0 as Obs.

It is assumed that all boundary effects have sufficiently decayed at x0=xmin. The data around the zero crossing of t d(t^2)/dt - (where c=0.3 by default) 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)
- postfix (str): Postfix of measurement file (Default: ms)
- c (float): Constant that defines the flow scale. Default 0.3 for w_0, choose 2./3 for w_1.
- 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 w0 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 w0 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

- w0 (Obs): Extracted w0

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

    666def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
    667    """Read the topologial charge based on openQCD gradient flow measurements.
    668
    669    Parameters
    670    ----------
    671    path : str
    672        path of the measurement files
    673    prefix : str
    674        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
    675        Ignored if file names are passed explicitly via keyword files.
    676    c : double
    677        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
    678    dtr_cnfg : int
    679        (optional) parameter that specifies the number of measurements
    680        between two configs.
    681        If it is not set, the distance between two measurements
    682        in the file is assumed to be the distance between two configurations.
    683    steps : int
    684        (optional) Distance between two configurations in units of trajectories /
    685         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    686    version : str
    687        Either openQCD or sfqcd, depending on the data.
    688    L : int
    689        spatial length of the lattice in L/a.
    690        HAS to be set if version != sfqcd, since openQCD does not provide
    691        this in the header
    692    r_start : list
    693        list which contains the first config to be read for each replicum.
    694    r_stop : list
    695        list which contains the last config to be read for each replicum.
    696    files : list
    697        specify the exact files that need to be read
    698        from path, practical if e.g. only one replicum is needed
    699    postfix : str
    700        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
    701    names : list
    702        Alternative labeling for replicas/ensembles.
    703        Has to have the appropriate length.
    704    Zeuthen_flow : bool
    705        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
    706        for version=='sfqcd' If False, the Wilson flow is used.
    707    integer_charge : bool
    708        If True, the charge is rounded towards the nearest integer on each config.
    709
    710    Returns
    711    -------
    712    result : Obs
    713        Read topological charge
    714    """
    715
    716    return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)

Read the topologial charge based on openQCD gradient flow measurements.

###### Parameters

- path (str): path of the measurement files
- prefix (str): prefix of the measurement files, e.g. _id0_r0.ms.dat. Ignored if file names are passed explicitly via keyword files.
- c (double): Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
- dtr_cnfg (int): (optional) parameter that specifies the number of measurements between two configs. If it is not set, the distance between two measurements in the file is assumed to be the distance between two configurations.
- steps (int): (optional) Distance between two configurations in units of trajectories / cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
- version (str): Either openQCD or sfqcd, depending on the data.
- L (int): spatial length of the lattice in L/a. HAS to be set if version != sfqcd, since openQCD does not provide this in the header
- r_start (list): list which contains the first config to be read for each replicum.
- r_stop (list): list which contains the last config to be read for each replicum.
- files (list): specify the exact files that need to be read from path, practical if e.g. only one replicum is needed
- 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

    719def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
    720    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
    721
    722    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.
    723
    724    Parameters
    725    ----------
    726    path : str
    727        path of the measurement files
    728    prefix : str
    729        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
    730        Ignored if file names are passed explicitly via keyword files.
    731    c : double
    732        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
    733    dtr_cnfg : int
    734        (optional) parameter that specifies the number of measurements
    735        between two configs.
    736        If it is not set, the distance between two measurements
    737        in the file is assumed to be the distance between two configurations.
    738    steps : int
    739        (optional) Distance between two configurations in units of trajectories /
    740         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    741    r_start : list
    742        list which contains the first config to be read for each replicum.
    743    r_stop : list
    744        list which contains the last config to be read for each replicum.
    745    files : list
    746        specify the exact files that need to be read
    747        from path, practical if e.g. only one replicum is needed
    748    names : list
    749        Alternative labeling for replicas/ensembles.
    750        Has to have the appropriate length.
    751    postfix : str
    752        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
    753    Zeuthen_flow : bool
    754        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
    755    """
    756
    757    if c != 0.3:
    758        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
    759
    760    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)
    761    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)
    762    L = plaq.tag["L"]
    763    T = plaq.tag["T"]
    764
    765    if T != L:
    766        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
    767
    768    if Zeuthen_flow is not True:
    769        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
    770
    771    t = (c * L) ** 2 / 8
    772
    773    normdict = {4: 0.012341170468270,
    774                6: 0.010162691462430,
    775                8: 0.009031614807931,
    776                10: 0.008744966371393,
    777                12: 0.008650917856809,
    778                14: 8.611154391267955E-03,
    779                16: 0.008591758449508,
    780                20: 0.008575359627103,
    781                24: 0.008569387847540,
    782                28: 8.566803713382559E-03,
    783                32: 0.008565541650006,
    784                40: 8.564480684962046E-03,
    785                48: 8.564098025073460E-03,
    786                64: 8.563853943383087E-03}
    787
    788    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]

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

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

###### Parameters

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

    1063def qtop_projection(qtop, target=0):
    1064    """Returns the projection to the topological charge sector defined by target.
    1065
    1066    Parameters
    1067    ----------
    1068    path : Obs
    1069        Topological charge.
    1070    target : int
    1071        Specifies the topological sector to be reweighted to (default 0)
    1072
    1073    Returns
    1074    -------
    1075    reto : Obs
    1076        projection to the topological charge sector defined by target
    1077    """
    1078    if qtop.reweighted:
    1079        raise Exception('You can not use a reweighted observable for reweighting!')
    1080
    1081    proj_qtop = []
    1082    for n in qtop.deltas:
    1083        proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
    1084
    1085    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
    1086    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

    1089def read_qtop_sector(path, prefix, c, target=0, **kwargs):
    1090    """Constructs reweighting factors to a specified topological sector.
    1091
    1092    Parameters
    1093    ----------
    1094    path : str
    1095        path of the measurement files
    1096    prefix : str
    1097        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
    1098    c : double
    1099        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
    1100    target : int
    1101        Specifies the topological sector to be reweighted to (default 0)
    1102    dtr_cnfg : int
    1103        (optional) parameter that specifies the number of trajectories
    1104        between two configs.
    1105        if it is not set, the distance between two measurements
    1106        in the file is assumed to be the distance between two configurations.
    1107    steps : int
    1108        (optional) Distance between two configurations in units of trajectories /
    1109         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    1110    version : str
    1111        version string of the openQCD (sfqcd) version used to create
    1112        the ensemble. Default is 2.0. May also be set to sfqcd.
    1113    L : int
    1114        spatial length of the lattice in L/a.
    1115        HAS to be set if version != sfqcd, since openQCD does not provide
    1116        this in the header
    1117    r_start : list
    1118        offset of the first ensemble, making it easier to match
    1119        later on with other Obs
    1120    r_stop : list
    1121        last configurations that need to be read (per replicum)
    1122    files : list
    1123        specify the exact files that need to be read
    1124        from path, practical if e.g. only one replicum is needed
    1125    names : list
    1126        Alternative labeling for replicas/ensembles.
    1127        Has to have the appropriate length
    1128    Zeuthen_flow : bool
    1129        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
    1130        for version=='sfqcd' If False, the Wilson flow is used.
    1131
    1132    Returns
    1133    -------
    1134    reto : Obs
    1135        projection to the topological charge sector defined by target
    1136    """
    1137
    1138    if not isinstance(target, int):
    1139        raise Exception("'target' has to be an integer.")
    1140
    1141    kwargs['integer_charge'] = True
    1142    qtop = read_qtop(path, prefix, c, **kwargs)
    1143
    1144    return qtop_projection(qtop, target=target)

Constructs reweighting factors to a specified topological sector.

###### Parameters

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

###### 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

    1147def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
    1148    """
    1149    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
    1150
    1151    Parameters
    1152    ----------
    1153    path : str
    1154        The directory to search for the files in.
    1155    prefix : str
    1156        The prefix to match the files against.
    1157    qc : str
    1158        The quark combination extension to match the files against.
    1159    corr : str
    1160        The correlator to extract data for.
    1161    sep : str, optional
    1162        The separator to use when parsing the replika names.
    1163    **kwargs
    1164        Additional keyword arguments. The following keyword arguments are recognized:
    1165
    1166        - names (List[str]): A list of names to use for the replicas.
    1167        - files (List[str]): A list of files to read data from.
    1168        - idl (List[List[int]]): A list of idls per replicum, resticting data to the idls given.
    1169
    1170    Returns
    1171    -------
    1172    Corr
    1173        A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
    1174    or
    1175    CObs
    1176        A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
    1177
    1178
    1179    Raises
    1180    ------
    1181    FileNotFoundError
    1182        If no files matching the specified prefix and quark combination extension are found in the specified directory.
    1183    IOError
    1184        If there is an error reading a file.
    1185    struct.error
    1186        If there is an error unpacking binary data.
    1187    """
    1188
    1189    # found = []
    1190    files = []
    1191    names = []
    1192
    1193    # test if the input is correct
    1194    if qc not in ['dd', 'ud', 'du', 'uu']:
    1195        raise Exception("Unknown quark conbination!")
    1196
    1197    if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
    1198        raise Exception("Unknown correlator!")
    1199
    1200    if "files" in kwargs:
    1201        known_files = kwargs.get("files")
    1202    else:
    1203        known_files = []
    1204    files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
    1205
    1206    if "names" in kwargs:
    1207        names = kwargs.get("names")
    1208    else:
    1209        for f in files:
    1210            if not sep == "":
    1211                se = f.split(".")[0]
    1212                for s in f.split(".")[1:-2]:
    1213                    se += "." + s
    1214                names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
    1215            else:
    1216                names.append(prefix)
    1217    if 'idl' in kwargs:
    1218        expected_idl = kwargs.get('idl')
    1219    names = sorted(names)
    1220    files = sorted(files)
    1221
    1222    cnfgs = []
    1223    realsamples = []
    1224    imagsamples = []
    1225    repnum = 0
    1226    for file in files:
    1227        with open(path + "/" + file, "rb") as fp:
    1228
    1229            t = fp.read(8)
    1230            kappa = struct.unpack('d', t)[0]
    1231            t = fp.read(8)
    1232            csw = struct.unpack('d', t)[0]
    1233            t = fp.read(8)
    1234            dF = struct.unpack('d', t)[0]
    1235            t = fp.read(8)
    1236            zF = struct.unpack('d', t)[0]
    1237
    1238            t = fp.read(4)
    1239            tmax = struct.unpack('i', t)[0]
    1240            t = fp.read(4)
    1241            bnd = struct.unpack('i', t)[0]
    1242
    1243            placesBI = ["gS", "gP",
    1244                        "gA", "gV",
    1245                        "gVt", "lA",
    1246                        "lV", "lVt",
    1247                        "lT", "lTt"]
    1248            placesBB = ["g1", "l1"]
    1249
    1250            # the chunks have the following structure:
    1251            # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
    1252
    1253            chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
    1254            packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
    1255            cnfgs.append([])
    1256            realsamples.append([])
    1257            imagsamples.append([])
    1258            for t in range(tmax):
    1259                realsamples[repnum].append([])
    1260                imagsamples[repnum].append([])
    1261            if 'idl' in kwargs:
    1262                left_idl = set(expected_idl[repnum])
    1263            while True:
    1264                cnfgt = fp.read(chunksize)
    1265                if not cnfgt:
    1266                    break
    1267                asascii = struct.unpack(packstr, cnfgt)
    1268                cnfg = asascii[0]
    1269                idl_wanted = True
    1270                if 'idl' in kwargs:
    1271                    idl_wanted = (cnfg in expected_idl[repnum])
    1272                    left_idl = left_idl - set([cnfg])
    1273                if idl_wanted:
    1274                    cnfgs[repnum].append(cnfg)
    1275
    1276                    if corr not in placesBB:
    1277                        tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
    1278                    else:
    1279                        tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
    1280
    1281                    corrres = [[], []]
    1282                    for i in range(len(tmpcorr)):
    1283                        corrres[i % 2].append(tmpcorr[i])
    1284                    for t in range(int(len(tmpcorr) / 2)):
    1285                        realsamples[repnum][t].append(corrres[0][t])
    1286                    for t in range(int(len(tmpcorr) / 2)):
    1287                        imagsamples[repnum][t].append(corrres[1][t])
    1288            if 'idl' in kwargs:
    1289                left_idl = list(left_idl)
    1290                if expected_idl[repnum] == left_idl:
    1291                    raise ValueError("None of the idls searched for were found in replikum of file " + file)
    1292                elif len(left_idl) > 0:
    1293                    warnings.warn('Could not find idls ' + str(left_idl) + ' in replikum of file ' + file, UserWarning)
    1294        repnum += 1
    1295    s = "Read correlator " + corr + " from " + str(repnum) + " replika with idls" + str(realsamples[0][t])
    1296    for rep in range(1, repnum):
    1297        s += ", " + str(realsamples[rep][t])
    1298    print(s)
    1299    print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
    1300
    1301    # we have the data now... but we need to re format the whole thing and put it into Corr objects.
    1302
    1303    compObs = []
    1304
    1305    for t in range(int(len(tmpcorr) / 2)):
    1306        compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
    1307                            Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
    1308
    1309    if len(compObs) == 1:
    1310        return compObs[0]
    1311    else:
    1312        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.
  - files (List[str]): A list of files to read data from.
  - idl (List[List[int]]): A list of idls per replicum, resticting data to the idls given.

###### 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.