# 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_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfix='ms', **kwargs):
     233    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
     234
     235    It is assumed that all boundary effects have
     236    sufficiently decayed at x0=xmin.
     237    The data around the zero crossing of t^2<E> - 0.3
     238    is fitted with a linear function
     239    from which the exact root is extracted.
     240
     241    It is assumed that one measurement is performed for each config.
     242    If this is not the case, the resulting idl, as well as the handling
     243    of r_start, r_stop and r_step is wrong and the user has to correct
     244    this in the resulting observable.
     245
     246    Parameters
     247    ----------
     248    path : str
     249        Path to .ms.dat files
     250    prefix : str
     251        Ensemble prefix
     252    dtr_read : int
     253        Determines how many trajectories should be skipped
     254        when reading the ms.dat files.
     255        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
     256    xmin : int
     257        First timeslice where the boundary
     258        effects have sufficiently decayed.
     259    spatial_extent : int
     260        spatial extent of the lattice, required for normalization.
     261    fit_range : int
     262        Number of data points left and right of the zero
     263        crossing to be included in the linear fit. (Default: 5)
     264    postfix : str
     265        Postfix of measurement file (Default: ms)
     266    r_start : list
     267        list which contains the first config to be read for each replicum.
     268    r_stop : list
     269        list which contains the last config to be read for each replicum.
     270    r_step : int
     271        integer that defines a fixed step size between two measurements (in units of configs)
     272        If not given, r_step=1 is assumed.
     273    plaquette : bool
     274        If true extract the plaquette estimate of t0 instead.
     275    names : list
     276        list of names that is assigned to the data according according
     277        to the order in the file list. Use careful, if you do not provide file names!
     278    files : list
     279        list which contains the filenames to be read. No automatic detection of
     280        files performed if given.
     281    plot_fit : bool
     282        If true, the fit for the extraction of t0 is shown together with the data.
     283    assume_thermalization : bool
     284        If True: If the first record divided by the distance between two measurements is larger than
     285        1, it is assumed that this is due to thermalization and the first measurement belongs
     286        to the first config (default).
     287        If False: The config numbers are assumed to be traj_number // difference
     288
     289    Returns
     290    -------
     291    t0 : Obs
     292        Extracted t0
     293    """
     294
     295    if 'files' in kwargs:
     296        known_files = kwargs.get('files')
     297    else:
     298        known_files = []
     299
     300    ls = _find_files(path, prefix, postfix, 'dat', known_files=known_files)
     301
     302    replica = len(ls)
     303
     304    if 'r_start' in kwargs:
     305        r_start = kwargs.get('r_start')
     306        if len(r_start) != replica:
     307            raise Exception('r_start does not match number of replicas')
     308        r_start = [o if o else None for o in r_start]
     309    else:
     310        r_start = [None] * replica
     311
     312    if 'r_stop' in kwargs:
     313        r_stop = kwargs.get('r_stop')
     314        if len(r_stop) != replica:
     315            raise Exception('r_stop does not match number of replicas')
     316    else:
     317        r_stop = [None] * replica
     318
     319    if 'r_step' in kwargs:
     320        r_step = kwargs.get('r_step')
     321    else:
     322        r_step = 1
     323
     324    print('Extract t0 from', prefix, ',', replica, 'replica')
     325
     326    if 'names' in kwargs:
     327        rep_names = kwargs.get('names')
     328    else:
     329        rep_names = []
     330        for entry in ls:
     331            truncated_entry = entry.split('.')[0]
     332            idx = truncated_entry.index('r')
     333            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
     334
     335    Ysum = []
     336
     337    configlist = []
     338    r_start_index = []
     339    r_stop_index = []
     340
     341    for rep in range(replica):
     342
     343        with open(path + '/' + ls[rep], 'rb') as fp:
     344            t = fp.read(12)
     345            header = struct.unpack('iii', t)
     346            if rep == 0:
     347                dn = header[0]
     348                nn = header[1]
     349                tmax = header[2]
     350            elif dn != header[0] or nn != header[1] or tmax != header[2]:
     351                raise Exception('Replica parameters do not match.')
     352
     353            t = fp.read(8)
     354            if rep == 0:
     355                eps = struct.unpack('d', t)[0]
     356                print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
     357            elif eps != struct.unpack('d', t)[0]:
     358                raise Exception('Values for eps do not match among replica.')
     359
     360            Ysl = []
     361
     362            configlist.append([])
     363            while True:
     364                t = fp.read(4)
     365                if (len(t) < 4):
     366                    break
     367                nc = struct.unpack('i', t)[0]
     368                configlist[-1].append(nc)
     369
     370                t = fp.read(8 * tmax * (nn + 1))
     371                if kwargs.get('plaquette'):
     372                    if nc % dtr_read == 0:
     373                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
     374                t = fp.read(8 * tmax * (nn + 1))
     375                if not kwargs.get('plaquette'):
     376                    if nc % dtr_read == 0:
     377                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
     378                t = fp.read(8 * tmax * (nn + 1))
     379
     380        Ysum.append([])
     381        for i, item in enumerate(Ysl):
     382            Ysum[-1].append([np.mean(item[current + xmin:
     383                             current + tmax - xmin])
     384                            for current in range(0, len(item), tmax)])
     385
     386        diffmeas = configlist[-1][-1] - configlist[-1][-2]
     387        configlist[-1] = [item // diffmeas for item in configlist[-1]]
     388        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
     389            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
     390            offset = configlist[-1][0] - 1
     391            configlist[-1] = [item - offset for item in configlist[-1]]
     392
     393        if r_start[rep] is None:
     394            r_start_index.append(0)
     395        else:
     396            try:
     397                r_start_index.append(configlist[-1].index(r_start[rep]))
     398            except ValueError:
     399                raise Exception('Config %d not in file with range [%d, %d]' % (
     400                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
     401
     402        if r_stop[rep] is None:
     403            r_stop_index.append(len(configlist[-1]) - 1)
     404        else:
     405            try:
     406                r_stop_index.append(configlist[-1].index(r_stop[rep]))
     407            except ValueError:
     408                raise Exception('Config %d not in file with range [%d, %d]' % (
     409                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
     410
     411    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
     412        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
     413    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
     414    if np.any([step != 1 for step in stepsizes]):
     415        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
     416
     417    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
     418    t2E_dict = {}
     419    for n in range(nn + 1):
     420        samples = []
     421        for nrep, rep in enumerate(Ysum):
     422            samples.append([])
     423            for cnfg in rep:
     424                samples[-1].append(cnfg[n])
     425            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
     426        new_obs = Obs(samples, rep_names, idl=idl)
     427        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
     428
     429    return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit'))
     430
     431
     432def _parse_array_openQCD2(d, n, size, wa, quadrupel=False):
     433    arr = []
     434    if d == 2:
     435        for i in range(n[0]):
     436            tmp = wa[i * n[1]:(i + 1) * n[1]]
     437            if quadrupel:
     438                tmp2 = []
     439                for j in range(0, len(tmp), 2):
     440                    tmp2.append(tmp[j])
     441                arr.append(tmp2)
     442            else:
     443                arr.append(np.asarray(tmp))
     444
     445    else:
     446        raise Exception('Only two-dimensional arrays supported!')
     447
     448    return arr
     449
     450
     451def _find_files(path, prefix, postfix, ext, known_files=[]):
     452    found = []
     453    files = []
     454
     455    if postfix != "":
     456        if postfix[-1] != ".":
     457            postfix = postfix + "."
     458        if postfix[0] != ".":
     459            postfix = "." + postfix
     460
     461    if ext[0] == ".":
     462        ext = ext[1:]
     463
     464    pattern = prefix + "*" + postfix + ext
     465
     466    for (dirpath, dirnames, filenames) in os.walk(path + "/"):
     467        found.extend(filenames)
     468        break
     469
     470    if known_files != []:
     471        for kf in known_files:
     472            if kf not in found:
     473                raise FileNotFoundError("Given file " + kf + " does not exist!")
     474
     475        return known_files
     476
     477    if not found:
     478        raise FileNotFoundError(f"Error, directory '{path}' not found")
     479
     480    for f in found:
     481        if fnmatch.fnmatch(f, pattern):
     482            files.append(f)
     483
     484    if files == []:
     485        raise Exception("No files found after pattern filter!")
     486
     487    files = sort_names(files)
     488    return files
     489
     490
     491def _read_array_openQCD2(fp):
     492    t = fp.read(4)
     493    d = struct.unpack('i', t)[0]
     494    t = fp.read(4 * d)
     495    n = struct.unpack('%di' % (d), t)
     496    t = fp.read(4)
     497    size = struct.unpack('i', t)[0]
     498    if size == 4:
     499        types = 'i'
     500    elif size == 8:
     501        types = 'd'
     502    elif size == 16:
     503        types = 'dd'
     504    else:
     505        raise Exception("Type for size '" + str(size) + "' not known.")
     506    m = n[0]
     507    for i in range(1, d):
     508        m *= n[i]
     509
     510    t = fp.read(m * size)
     511    tmp = struct.unpack('%d%s' % (m, types), t)
     512
     513    arr = _parse_array_openQCD2(d, n, size, tmp, quadrupel=True)
     514    return {'d': d, 'n': n, 'size': size, 'arr': arr}
     515
     516
     517def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
     518    """Read the topologial charge based on openQCD gradient flow measurements.
     519
     520    Parameters
     521    ----------
     522    path : str
     523        path of the measurement files
     524    prefix : str
     525        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
     526        Ignored if file names are passed explicitly via keyword files.
     527    c : double
     528        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
     529    dtr_cnfg : int
     530        (optional) parameter that specifies the number of measurements
     531        between two configs.
     532        If it is not set, the distance between two measurements
     533        in the file is assumed to be the distance between two configurations.
     534    steps : int
     535        (optional) Distance between two configurations in units of trajectories /
     536         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
     537    version : str
     538        Either openQCD or sfqcd, depending on the data.
     539    L : int
     540        spatial length of the lattice in L/a.
     541        HAS to be set if version != sfqcd, since openQCD does not provide
     542        this in the header
     543    r_start : list
     544        list which contains the first config to be read for each replicum.
     545    r_stop : list
     546        list which contains the last config to be read for each replicum.
     547    files : list
     548        specify the exact files that need to be read
     549        from path, practical if e.g. only one replicum is needed
     550    postfix : str
     551        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
     552    names : list
     553        Alternative labeling for replicas/ensembles.
     554        Has to have the appropriate length.
     555    Zeuthen_flow : bool
     556        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
     557        for version=='sfqcd' If False, the Wilson flow is used.
     558    integer_charge : bool
     559        If True, the charge is rounded towards the nearest integer on each config.
     560
     561    Returns
     562    -------
     563    result : Obs
     564        Read topological charge
     565    """
     566
     567    return _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version=version, obspos=0, **kwargs)
     568
     569
     570def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
     571    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
     572
     573    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.
     574
     575    Parameters
     576    ----------
     577    path : str
     578        path of the measurement files
     579    prefix : str
     580        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
     581        Ignored if file names are passed explicitly via keyword files.
     582    c : double
     583        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
     584    dtr_cnfg : int
     585        (optional) parameter that specifies the number of measurements
     586        between two configs.
     587        If it is not set, the distance between two measurements
     588        in the file is assumed to be the distance between two configurations.
     589    steps : int
     590        (optional) Distance between two configurations in units of trajectories /
     591         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
     592    r_start : list
     593        list which contains the first config to be read for each replicum.
     594    r_stop : list
     595        list which contains the last config to be read for each replicum.
     596    files : list
     597        specify the exact files that need to be read
     598        from path, practical if e.g. only one replicum is needed
     599    names : list
     600        Alternative labeling for replicas/ensembles.
     601        Has to have the appropriate length.
     602    postfix : str
     603        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
     604    Zeuthen_flow : bool
     605        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
     606    """
     607
     608    if c != 0.3:
     609        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
     610
     611    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)
     612    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)
     613    L = plaq.tag["L"]
     614    T = plaq.tag["T"]
     615
     616    if T != L:
     617        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
     618
     619    if Zeuthen_flow is not True:
     620        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
     621
     622    t = (c * L) ** 2 / 8
     623
     624    normdict = {4: 0.012341170468270,
     625                6: 0.010162691462430,
     626                8: 0.009031614807931,
     627                10: 0.008744966371393,
     628                12: 0.008650917856809,
     629                14: 8.611154391267955E-03,
     630                16: 0.008591758449508,
     631                20: 0.008575359627103,
     632                24: 0.008569387847540,
     633                28: 8.566803713382559E-03,
     634                32: 0.008565541650006,
     635                40: 8.564480684962046E-03,
     636                48: 8.564098025073460E-03,
     637                64: 8.563853943383087E-03}
     638
     639    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]
     640
     641
     642def _read_flow_obs(path, prefix, c, dtr_cnfg=1, version="openQCD", obspos=0, sum_t=True, **kwargs):
     643    """Read a flow observable based on openQCD gradient flow measurements.
     644
     645    Parameters
     646    ----------
     647    path : str
     648        path of the measurement files
     649    prefix : str
     650        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
     651        Ignored if file names are passed explicitly via keyword files.
     652    c : double
     653        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
     654    dtr_cnfg : int
     655        (optional) parameter that specifies the number of measurements
     656        between two configs.
     657        If it is not set, the distance between two measurements
     658        in the file is assumed to be the distance between two configurations.
     659    steps : int
     660        (optional) Distance between two configurations in units of trajectories /
     661         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
     662    version : str
     663        Either openQCD or sfqcd, depending on the data.
     664    obspos : int
     665        position of the obeservable in the measurement file. Only relevant for sfqcd files.
     666    sum_t : bool
     667        If true sum over all timeslices, if false only take the value at T/2.
     668    L : int
     669        spatial length of the lattice in L/a.
     670        HAS to be set if version != sfqcd, since openQCD does not provide
     671        this in the header
     672    r_start : list
     673        list which contains the first config to be read for each replicum.
     674    r_stop : list
     675        list which contains the last config to be read for each replicum.
     676    files : list
     677        specify the exact files that need to be read
     678        from path, practical if e.g. only one replicum is needed
     679    names : list
     680        Alternative labeling for replicas/ensembles.
     681        Has to have the appropriate length.
     682    postfix : str
     683        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
     684    Zeuthen_flow : bool
     685        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
     686        for version=='sfqcd' If False, the Wilson flow is used.
     687    integer_charge : bool
     688        If True, the charge is rounded towards the nearest integer on each config.
     689
     690    Returns
     691    -------
     692    result : Obs
     693        flow observable specified
     694    """
     695    known_versions = ["openQCD", "sfqcd"]
     696
     697    if version not in known_versions:
     698        raise Exception("Unknown openQCD version.")
     699    if "steps" in kwargs:
     700        steps = kwargs.get("steps")
     701    if version == "sfqcd":
     702        if "L" in kwargs:
     703            supposed_L = kwargs.get("L")
     704        else:
     705            supposed_L = None
     706        postfix = "gfms"
     707    else:
     708        if "L" not in kwargs:
     709            raise Exception("This version of openQCD needs you to provide the spatial length of the lattice as parameter 'L'.")
     710        else:
     711            L = kwargs.get("L")
     712        postfix = "ms"
     713
     714    if "postfix" in kwargs:
     715        postfix = kwargs.get("postfix")
     716
     717    if "files" in kwargs:
     718        known_files = kwargs.get("files")
     719    else:
     720        known_files = []
     721
     722    files = _find_files(path, prefix, postfix, "dat", known_files=known_files)
     723
     724    if 'r_start' in kwargs:
     725        r_start = kwargs.get('r_start')
     726        if len(r_start) != len(files):
     727            raise Exception('r_start does not match number of replicas')
     728        r_start = [o if o else None for o in r_start]
     729    else:
     730        r_start = [None] * len(files)
     731
     732    if 'r_stop' in kwargs:
     733        r_stop = kwargs.get('r_stop')
     734        if len(r_stop) != len(files):
     735            raise Exception('r_stop does not match number of replicas')
     736    else:
     737        r_stop = [None] * len(files)
     738    rep_names = []
     739
     740    zeuthen = kwargs.get('Zeuthen_flow', False)
     741    if zeuthen and version not in ['sfqcd']:
     742        raise Exception('Zeuthen flow can only be used for version==sfqcd')
     743
     744    r_start_index = []
     745    r_stop_index = []
     746    deltas = []
     747    configlist = []
     748    if not zeuthen:
     749        obspos += 8
     750    for rep, file in enumerate(files):
     751        with open(path + "/" + file, "rb") as fp:
     752
     753            Q = []
     754            traj_list = []
     755            if version in ['sfqcd']:
     756                t = fp.read(12)
     757                header = struct.unpack('<iii', t)
     758                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)
     759                ncs = header[1]  # number of different values for c in t_flow=1/8 c² L² -> measurements done for ncs c's
     760                tmax = header[2]  # lattice T/a
     761
     762                t = fp.read(12)
     763                Ls = struct.unpack('<iii', t)
     764                if (Ls[0] == Ls[1] and Ls[1] == Ls[2]):
     765                    L = Ls[0]
     766                    if not (supposed_L == L) and supposed_L:
     767                        raise Exception("It seems the length given in the header and by you contradict each other")
     768                else:
     769                    raise Exception("Found more than one spatial length in header!")
     770
     771                t = fp.read(16)
     772                header2 = struct.unpack('<dd', t)
     773                tol = header2[0]
     774                cmax = header2[1]  # highest value of c used
     775
     776                if c > cmax:
     777                    raise Exception('Flow has been determined between c=0 and c=%lf with tolerance %lf' % (cmax, tol))
     778
     779                if (zthfl == 2):
     780                    nfl = 2  # number of flows
     781                else:
     782                    nfl = 1
     783                iobs = 8 * nfl  # number of flow observables calculated
     784
     785                while True:
     786                    t = fp.read(4)
     787                    if (len(t) < 4):
     788                        break
     789                    traj_list.append(struct.unpack('i', t)[0])   # trajectory number when measurement was done
     790
     791                    for j in range(ncs + 1):
     792                        for i in range(iobs):
     793                            t = fp.read(8 * tmax)
     794                            if (i == obspos):  # determines the flow observable -> i=0 <-> Zeuthen flow
     795                                Q.append(struct.unpack('d' * tmax, t))
     796
     797            else:
     798                t = fp.read(12)
     799                header = struct.unpack('<iii', t)
     800                # step size in integration steps "dnms"
     801                dn = header[0]
     802                # number of measurements, so "ntot"/dn
     803                nn = header[1]
     804                # lattice T/a
     805                tmax = header[2]
     806
     807                t = fp.read(8)
     808                eps = struct.unpack('d', t)[0]
     809
     810                while True:
     811                    t = fp.read(4)
     812                    if (len(t) < 4):
     813                        break
     814                    traj_list.append(struct.unpack('i', t)[0])
     815                    # Wsl
     816                    t = fp.read(8 * tmax * (nn + 1))
     817                    # Ysl
     818                    t = fp.read(8 * tmax * (nn + 1))
     819                    # Qsl, which is asked for in this method
     820                    t = fp.read(8 * tmax * (nn + 1))
     821                    # unpack the array of Qtops,
     822                    # on each timeslice t=0,...,tmax-1 and the
     823                    # measurement number in = 0...nn (see README.qcd1)
     824                    tmpd = struct.unpack('d' * tmax * (nn + 1), t)
     825                    Q.append(tmpd)
     826
     827        if len(np.unique(np.diff(traj_list))) != 1:
     828            raise Exception("Irregularities in stepsize found")
     829        else:
     830            if 'steps' in kwargs:
     831                if steps != traj_list[1] - traj_list[0]:
     832                    raise Exception("steps and the found stepsize are not the same")
     833            else:
     834                steps = traj_list[1] - traj_list[0]
     835
     836        configlist.append([tr // steps // dtr_cnfg for tr in traj_list])
     837        if configlist[-1][0] > 1:
     838            offset = configlist[-1][0] - 1
     839            warnings.warn('Assume thermalization and that the first measurement belongs to the first config. Offset = %d configs (%d trajectories / cycles)' % (
     840                offset, offset * steps))
     841            configlist[-1] = [item - offset for item in configlist[-1]]
     842
     843        if r_start[rep] is None:
     844            r_start_index.append(0)
     845        else:
     846            try:
     847                r_start_index.append(configlist[-1].index(r_start[rep]))
     848            except ValueError:
     849                raise Exception('Config %d not in file with range [%d, %d]' % (
     850                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
     851
     852        if r_stop[rep] is None:
     853            r_stop_index.append(len(configlist[-1]) - 1)
     854        else:
     855            try:
     856                r_stop_index.append(configlist[-1].index(r_stop[rep]))
     857            except ValueError:
     858                raise Exception('Config %d not in file with range [%d, %d]' % (
     859                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
     860
     861        if version in ['sfqcd']:
     862            cstepsize = cmax / ncs
     863            index_aim = round(c / cstepsize)
     864        else:
     865            t_aim = (c * L) ** 2 / 8
     866            index_aim = round(t_aim / eps / dn)
     867
     868        Q_sum = []
     869        for i, item in enumerate(Q):
     870            if sum_t is True:
     871                Q_sum.append([sum(item[current:current + tmax])
     872                             for current in range(0, len(item), tmax)])
     873            else:
     874                Q_sum.append([item[int(tmax / 2)]])
     875        Q_top = []
     876        if version in ['sfqcd']:
     877            for i in range(len(Q_sum) // (ncs + 1)):
     878                Q_top.append(Q_sum[i * (ncs + 1) + index_aim][0])
     879        else:
     880            for i in range(len(Q) // dtr_cnfg):
     881                Q_top.append(Q_sum[dtr_cnfg * i][index_aim])
     882        if len(Q_top) != len(traj_list) // dtr_cnfg:
     883            raise Exception("qtops and traj_list dont have the same length")
     884
     885        if kwargs.get('integer_charge', False):
     886            Q_top = [round(q) for q in Q_top]
     887
     888        truncated_file = file[:-len(postfix)]
     889
     890        if "names" not in kwargs:
     891            try:
     892                idx = truncated_file.index('r')
     893            except Exception:
     894                if "names" not in kwargs:
     895                    raise Exception("Automatic recognition of replicum failed, please enter the key word 'names'.")
     896            ens_name = truncated_file[:idx]
     897            rep_names.append(ens_name + '|' + truncated_file[idx:].split(".")[0])
     898        else:
     899            names = kwargs.get("names")
     900            rep_names = names
     901
     902        deltas.append(Q_top)
     903
     904    rep_names = sort_names(rep_names)
     905
     906    idl = [range(int(configlist[rep][r_start_index[rep]]), int(configlist[rep][r_stop_index[rep]]) + 1, 1) for rep in range(len(deltas))]
     907    deltas = [deltas[nrep][r_start_index[nrep]:r_stop_index[nrep] + 1] for nrep in range(len(deltas))]
     908    result = Obs(deltas, rep_names, idl=idl)
     909    result.tag = {"T": tmax - 1,
     910                  "L": L}
     911    return result
     912
     913
     914def qtop_projection(qtop, target=0):
     915    """Returns the projection to the topological charge sector defined by target.
     916
     917    Parameters
     918    ----------
     919    path : Obs
     920        Topological charge.
     921    target : int
     922        Specifies the topological sector to be reweighted to (default 0)
     923
     924    Returns
     925    -------
     926    reto : Obs
     927        projection to the topological charge sector defined by target
     928    """
     929    if qtop.reweighted:
     930        raise Exception('You can not use a reweighted observable for reweighting!')
     931
     932    proj_qtop = []
     933    for n in qtop.deltas:
     934        proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
     935
     936    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
     937    return reto
     938
     939
     940def read_qtop_sector(path, prefix, c, target=0, **kwargs):
     941    """Constructs reweighting factors to a specified topological sector.
     942
     943    Parameters
     944    ----------
     945    path : str
     946        path of the measurement files
     947    prefix : str
     948        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
     949    c : double
     950        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
     951    target : int
     952        Specifies the topological sector to be reweighted to (default 0)
     953    dtr_cnfg : int
     954        (optional) parameter that specifies the number of trajectories
     955        between two configs.
     956        if it is not set, the distance between two measurements
     957        in the file is assumed to be the distance between two configurations.
     958    steps : int
     959        (optional) Distance between two configurations in units of trajectories /
     960         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
     961    version : str
     962        version string of the openQCD (sfqcd) version used to create
     963        the ensemble. Default is 2.0. May also be set to sfqcd.
     964    L : int
     965        spatial length of the lattice in L/a.
     966        HAS to be set if version != sfqcd, since openQCD does not provide
     967        this in the header
     968    r_start : list
     969        offset of the first ensemble, making it easier to match
     970        later on with other Obs
     971    r_stop : list
     972        last configurations that need to be read (per replicum)
     973    files : list
     974        specify the exact files that need to be read
     975        from path, practical if e.g. only one replicum is needed
     976    names : list
     977        Alternative labeling for replicas/ensembles.
     978        Has to have the appropriate length
     979    Zeuthen_flow : bool
     980        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
     981        for version=='sfqcd' If False, the Wilson flow is used.
     982
     983    Returns
     984    -------
     985    reto : Obs
     986        projection to the topological charge sector defined by target
     987    """
     988
     989    if not isinstance(target, int):
     990        raise Exception("'target' has to be an integer.")
     991
     992    kwargs['integer_charge'] = True
     993    qtop = read_qtop(path, prefix, c, **kwargs)
     994
     995    return qtop_projection(qtop, target=target)
     996
     997
     998def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
     999    """
    1000    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
    1001
    1002    Parameters
    1003    ----------
    1004    path : str
    1005        The directory to search for the files in.
    1006    prefix : str
    1007        The prefix to match the files against.
    1008    qc : str
    1009        The quark combination extension to match the files against.
    1010    corr : str
    1011        The correlator to extract data for.
    1012    sep : str, optional
    1013        The separator to use when parsing the replika names.
    1014    **kwargs
    1015        Additional keyword arguments. The following keyword arguments are recognized:
    1016
    1017        - names (List[str]): A list of names to use for the replicas.
    1018
    1019    Returns
    1020    -------
    1021    Corr
    1022        A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
    1023    or
    1024    CObs
    1025        A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
    1026
    1027
    1028    Raises
    1029    ------
    1030    FileNotFoundError
    1031        If no files matching the specified prefix and quark combination extension are found in the specified directory.
    1032    IOError
    1033        If there is an error reading a file.
    1034    struct.error
    1035        If there is an error unpacking binary data.
    1036    """
    1037
    1038    # found = []
    1039    files = []
    1040    names = []
    1041
    1042    # test if the input is correct
    1043    if qc not in ['dd', 'ud', 'du', 'uu']:
    1044        raise Exception("Unknown quark conbination!")
    1045
    1046    if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
    1047        raise Exception("Unknown correlator!")
    1048
    1049    if "files" in kwargs:
    1050        known_files = kwargs.get("files")
    1051    else:
    1052        known_files = []
    1053    files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
    1054
    1055    if "names" in kwargs:
    1056        names = kwargs.get("names")
    1057    else:
    1058        for f in files:
    1059            if not sep == "":
    1060                se = f.split(".")[0]
    1061                for s in f.split(".")[1:-2]:
    1062                    se += "." + s
    1063                names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
    1064            else:
    1065                names.append(prefix)
    1066
    1067    names = sorted(names)
    1068    files = sorted(files)
    1069
    1070    cnfgs = []
    1071    realsamples = []
    1072    imagsamples = []
    1073    repnum = 0
    1074    for file in files:
    1075        with open(path + "/" + file, "rb") as fp:
    1076
    1077            t = fp.read(8)
    1078            kappa = struct.unpack('d', t)[0]
    1079            t = fp.read(8)
    1080            csw = struct.unpack('d', t)[0]
    1081            t = fp.read(8)
    1082            dF = struct.unpack('d', t)[0]
    1083            t = fp.read(8)
    1084            zF = struct.unpack('d', t)[0]
    1085
    1086            t = fp.read(4)
    1087            tmax = struct.unpack('i', t)[0]
    1088            t = fp.read(4)
    1089            bnd = struct.unpack('i', t)[0]
    1090
    1091            placesBI = ["gS", "gP",
    1092                        "gA", "gV",
    1093                        "gVt", "lA",
    1094                        "lV", "lVt",
    1095                        "lT", "lTt"]
    1096            placesBB = ["g1", "l1"]
    1097
    1098            # the chunks have the following structure:
    1099            # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
    1100
    1101            chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
    1102            packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
    1103            cnfgs.append([])
    1104            realsamples.append([])
    1105            imagsamples.append([])
    1106            for t in range(tmax):
    1107                realsamples[repnum].append([])
    1108                imagsamples[repnum].append([])
    1109
    1110            while True:
    1111                cnfgt = fp.read(chunksize)
    1112                if not cnfgt:
    1113                    break
    1114                asascii = struct.unpack(packstr, cnfgt)
    1115                cnfg = asascii[0]
    1116                cnfgs[repnum].append(cnfg)
    1117
    1118                if corr not in placesBB:
    1119                    tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
    1120                else:
    1121                    tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
    1122
    1123                corrres = [[], []]
    1124                for i in range(len(tmpcorr)):
    1125                    corrres[i % 2].append(tmpcorr[i])
    1126                for t in range(int(len(tmpcorr) / 2)):
    1127                    realsamples[repnum][t].append(corrres[0][t])
    1128                for t in range(int(len(tmpcorr) / 2)):
    1129                    imagsamples[repnum][t].append(corrres[1][t])
    1130        repnum += 1
    1131
    1132    s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
    1133    for rep in range(1, repnum):
    1134        s += ", " + str(len(realsamples[rep][t]))
    1135    s += " samples"
    1136    print(s)
    1137    print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
    1138
    1139    # we have the data now... but we need to re format the whole thing and put it into Corr objects.
    1140
    1141    compObs = []
    1142
    1143    for t in range(int(len(tmpcorr) / 2)):
    1144        compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
    1145                            Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
    1146
    1147    if len(compObs) == 1:
    1148        return compObs[0]
    1149    else:
    1150        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', **kwargs): View Source

    233def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, postfix='ms', **kwargs):
    234    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
    235
    236    It is assumed that all boundary effects have
    237    sufficiently decayed at x0=xmin.
    238    The data around the zero crossing of t^2<E> - 0.3
    239    is fitted with a linear function
    240    from which the exact root is extracted.
    241
    242    It is assumed that one measurement is performed for each config.
    243    If this is not the case, the resulting idl, as well as the handling
    244    of r_start, r_stop and r_step is wrong and the user has to correct
    245    this in the resulting observable.
    246
    247    Parameters
    248    ----------
    249    path : str
    250        Path to .ms.dat files
    251    prefix : str
    252        Ensemble prefix
    253    dtr_read : int
    254        Determines how many trajectories should be skipped
    255        when reading the ms.dat files.
    256        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
    257    xmin : int
    258        First timeslice where the boundary
    259        effects have sufficiently decayed.
    260    spatial_extent : int
    261        spatial extent of the lattice, required for normalization.
    262    fit_range : int
    263        Number of data points left and right of the zero
    264        crossing to be included in the linear fit. (Default: 5)
    265    postfix : str
    266        Postfix of measurement file (Default: ms)
    267    r_start : list
    268        list which contains the first config to be read for each replicum.
    269    r_stop : list
    270        list which contains the last config to be read for each replicum.
    271    r_step : int
    272        integer that defines a fixed step size between two measurements (in units of configs)
    273        If not given, r_step=1 is assumed.
    274    plaquette : bool
    275        If true extract the plaquette estimate of t0 instead.
    276    names : list
    277        list of names that is assigned to the data according according
    278        to the order in the file list. Use careful, if you do not provide file names!
    279    files : list
    280        list which contains the filenames to be read. No automatic detection of
    281        files performed if given.
    282    plot_fit : bool
    283        If true, the fit for the extraction of t0 is shown together with the data.
    284    assume_thermalization : bool
    285        If True: If the first record divided by the distance between two measurements is larger than
    286        1, it is assumed that this is due to thermalization and the first measurement belongs
    287        to the first config (default).
    288        If False: The config numbers are assumed to be traj_number // difference
    289
    290    Returns
    291    -------
    292    t0 : Obs
    293        Extracted t0
    294    """
    295
    296    if 'files' in kwargs:
    297        known_files = kwargs.get('files')
    298    else:
    299        known_files = []
    300
    301    ls = _find_files(path, prefix, postfix, 'dat', known_files=known_files)
    302
    303    replica = len(ls)
    304
    305    if 'r_start' in kwargs:
    306        r_start = kwargs.get('r_start')
    307        if len(r_start) != replica:
    308            raise Exception('r_start does not match number of replicas')
    309        r_start = [o if o else None for o in r_start]
    310    else:
    311        r_start = [None] * replica
    312
    313    if 'r_stop' in kwargs:
    314        r_stop = kwargs.get('r_stop')
    315        if len(r_stop) != replica:
    316            raise Exception('r_stop does not match number of replicas')
    317    else:
    318        r_stop = [None] * replica
    319
    320    if 'r_step' in kwargs:
    321        r_step = kwargs.get('r_step')
    322    else:
    323        r_step = 1
    324
    325    print('Extract t0 from', prefix, ',', replica, 'replica')
    326
    327    if 'names' in kwargs:
    328        rep_names = kwargs.get('names')
    329    else:
    330        rep_names = []
    331        for entry in ls:
    332            truncated_entry = entry.split('.')[0]
    333            idx = truncated_entry.index('r')
    334            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
    335
    336    Ysum = []
    337
    338    configlist = []
    339    r_start_index = []
    340    r_stop_index = []
    341
    342    for rep in range(replica):
    343
    344        with open(path + '/' + ls[rep], 'rb') as fp:
    345            t = fp.read(12)
    346            header = struct.unpack('iii', t)
    347            if rep == 0:
    348                dn = header[0]
    349                nn = header[1]
    350                tmax = header[2]
    351            elif dn != header[0] or nn != header[1] or tmax != header[2]:
    352                raise Exception('Replica parameters do not match.')
    353
    354            t = fp.read(8)
    355            if rep == 0:
    356                eps = struct.unpack('d', t)[0]
    357                print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
    358            elif eps != struct.unpack('d', t)[0]:
    359                raise Exception('Values for eps do not match among replica.')
    360
    361            Ysl = []
    362
    363            configlist.append([])
    364            while True:
    365                t = fp.read(4)
    366                if (len(t) < 4):
    367                    break
    368                nc = struct.unpack('i', t)[0]
    369                configlist[-1].append(nc)
    370
    371                t = fp.read(8 * tmax * (nn + 1))
    372                if kwargs.get('plaquette'):
    373                    if nc % dtr_read == 0:
    374                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
    375                t = fp.read(8 * tmax * (nn + 1))
    376                if not kwargs.get('plaquette'):
    377                    if nc % dtr_read == 0:
    378                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
    379                t = fp.read(8 * tmax * (nn + 1))
    380
    381        Ysum.append([])
    382        for i, item in enumerate(Ysl):
    383            Ysum[-1].append([np.mean(item[current + xmin:
    384                             current + tmax - xmin])
    385                            for current in range(0, len(item), tmax)])
    386
    387        diffmeas = configlist[-1][-1] - configlist[-1][-2]
    388        configlist[-1] = [item // diffmeas for item in configlist[-1]]
    389        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
    390            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
    391            offset = configlist[-1][0] - 1
    392            configlist[-1] = [item - offset for item in configlist[-1]]
    393
    394        if r_start[rep] is None:
    395            r_start_index.append(0)
    396        else:
    397            try:
    398                r_start_index.append(configlist[-1].index(r_start[rep]))
    399            except ValueError:
    400                raise Exception('Config %d not in file with range [%d, %d]' % (
    401                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
    402
    403        if r_stop[rep] is None:
    404            r_stop_index.append(len(configlist[-1]) - 1)
    405        else:
    406            try:
    407                r_stop_index.append(configlist[-1].index(r_stop[rep]))
    408            except ValueError:
    409                raise Exception('Config %d not in file with range [%d, %d]' % (
    410                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
    411
    412    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
    413        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
    414    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
    415    if np.any([step != 1 for step in stepsizes]):
    416        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
    417
    418    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
    419    t2E_dict = {}
    420    for n in range(nn + 1):
    421        samples = []
    422        for nrep, rep in enumerate(Ysum):
    423            samples.append([])
    424            for cnfg in rep:
    425                samples[-1].append(cnfg[n])
    426            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
    427        new_obs = Obs(samples, rep_names, idl=idl)
    428        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
    429
    430    return fit_t0(t2E_dict, fit_range, plot_fit=kwargs.get('plot_fit'))

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

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

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

###### Parameters

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

###### Returns

- t0 (Obs): Extracted t0

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

    518def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
    519    """Read the topologial charge based on openQCD gradient flow measurements.
    520
    521    Parameters
    522    ----------
    523    path : str
    524        path of the measurement files
    525    prefix : str
    526        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
    527        Ignored if file names are passed explicitly via keyword files.
    528    c : double
    529        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
    530    dtr_cnfg : int
    531        (optional) parameter that specifies the number of measurements
    532        between two configs.
    533        If it is not set, the distance between two measurements
    534        in the file is assumed to be the distance between two configurations.
    535    steps : int
    536        (optional) Distance between two configurations in units of trajectories /
    537         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    538    version : str
    539        Either openQCD or sfqcd, depending on the data.
    540    L : int
    541        spatial length of the lattice in L/a.
    542        HAS to be set if version != sfqcd, since openQCD does not provide
    543        this in the header
    544    r_start : list
    545        list which contains the first config to be read for each replicum.
    546    r_stop : list
    547        list which contains the last config to be read for each replicum.
    548    files : list
    549        specify the exact files that need to be read
    550        from path, practical if e.g. only one replicum is needed
    551    postfix : str
    552        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
    553    names : list
    554        Alternative labeling for replicas/ensembles.
    555        Has to have the appropriate length.
    556    Zeuthen_flow : bool
    557        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
    558        for version=='sfqcd' If False, the Wilson flow is used.
    559    integer_charge : bool
    560        If True, the charge is rounded towards the nearest integer on each config.
    561
    562    Returns
    563    -------
    564    result : Obs
    565        Read topological charge
    566    """
    567
    568    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

    571def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
    572    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
    573
    574    Note: The current implementation only works for c=0.3 and T=L. The definition of the coupling in 1607.06423 requires projection to topological charge zero which is not done within this function but has to be performed in a separate step.
    575
    576    Parameters
    577    ----------
    578    path : str
    579        path of the measurement files
    580    prefix : str
    581        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
    582        Ignored if file names are passed explicitly via keyword files.
    583    c : double
    584        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
    585    dtr_cnfg : int
    586        (optional) parameter that specifies the number of measurements
    587        between two configs.
    588        If it is not set, the distance between two measurements
    589        in the file is assumed to be the distance between two configurations.
    590    steps : int
    591        (optional) Distance between two configurations in units of trajectories /
    592         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    593    r_start : list
    594        list which contains the first config to be read for each replicum.
    595    r_stop : list
    596        list which contains the last config to be read for each replicum.
    597    files : list
    598        specify the exact files that need to be read
    599        from path, practical if e.g. only one replicum is needed
    600    names : list
    601        Alternative labeling for replicas/ensembles.
    602        Has to have the appropriate length.
    603    postfix : str
    604        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
    605    Zeuthen_flow : bool
    606        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
    607    """
    608
    609    if c != 0.3:
    610        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
    611
    612    plaq = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=6, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
    613    C2x1 = _read_flow_obs(path, prefix, c, dtr_cnfg=dtr_cnfg, version="sfqcd", obspos=7, sum_t=False, Zeuthen_flow=Zeuthen_flow, integer_charge=False, **kwargs)
    614    L = plaq.tag["L"]
    615    T = plaq.tag["T"]
    616
    617    if T != L:
    618        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
    619
    620    if Zeuthen_flow is not True:
    621        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
    622
    623    t = (c * L) ** 2 / 8
    624
    625    normdict = {4: 0.012341170468270,
    626                6: 0.010162691462430,
    627                8: 0.009031614807931,
    628                10: 0.008744966371393,
    629                12: 0.008650917856809,
    630                14: 8.611154391267955E-03,
    631                16: 0.008591758449508,
    632                20: 0.008575359627103,
    633                24: 0.008569387847540,
    634                28: 8.566803713382559E-03,
    635                32: 0.008565541650006,
    636                40: 8.564480684962046E-03,
    637                48: 8.564098025073460E-03,
    638                64: 8.563853943383087E-03}
    639
    640    return t * t * (5 / 3 * plaq - 1 / 12 * C2x1) / normdict[L]

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

    915def qtop_projection(qtop, target=0):
    916    """Returns the projection to the topological charge sector defined by target.
    917
    918    Parameters
    919    ----------
    920    path : Obs
    921        Topological charge.
    922    target : int
    923        Specifies the topological sector to be reweighted to (default 0)
    924
    925    Returns
    926    -------
    927    reto : Obs
    928        projection to the topological charge sector defined by target
    929    """
    930    if qtop.reweighted:
    931        raise Exception('You can not use a reweighted observable for reweighting!')
    932
    933    proj_qtop = []
    934    for n in qtop.deltas:
    935        proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
    936
    937    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
    938    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

    941def read_qtop_sector(path, prefix, c, target=0, **kwargs):
    942    """Constructs reweighting factors to a specified topological sector.
    943
    944    Parameters
    945    ----------
    946    path : str
    947        path of the measurement files
    948    prefix : str
    949        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
    950    c : double
    951        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
    952    target : int
    953        Specifies the topological sector to be reweighted to (default 0)
    954    dtr_cnfg : int
    955        (optional) parameter that specifies the number of trajectories
    956        between two configs.
    957        if it is not set, the distance between two measurements
    958        in the file is assumed to be the distance between two configurations.
    959    steps : int
    960        (optional) Distance between two configurations in units of trajectories /
    961         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
    962    version : str
    963        version string of the openQCD (sfqcd) version used to create
    964        the ensemble. Default is 2.0. May also be set to sfqcd.
    965    L : int
    966        spatial length of the lattice in L/a.
    967        HAS to be set if version != sfqcd, since openQCD does not provide
    968        this in the header
    969    r_start : list
    970        offset of the first ensemble, making it easier to match
    971        later on with other Obs
    972    r_stop : list
    973        last configurations that need to be read (per replicum)
    974    files : list
    975        specify the exact files that need to be read
    976        from path, practical if e.g. only one replicum is needed
    977    names : list
    978        Alternative labeling for replicas/ensembles.
    979        Has to have the appropriate length
    980    Zeuthen_flow : bool
    981        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
    982        for version=='sfqcd' If False, the Wilson flow is used.
    983
    984    Returns
    985    -------
    986    reto : Obs
    987        projection to the topological charge sector defined by target
    988    """
    989
    990    if not isinstance(target, int):
    991        raise Exception("'target' has to be an integer.")
    992
    993    kwargs['integer_charge'] = True
    994    qtop = read_qtop(path, prefix, c, **kwargs)
    995
    996    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

     999def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
    1000    """
    1001    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
    1002
    1003    Parameters
    1004    ----------
    1005    path : str
    1006        The directory to search for the files in.
    1007    prefix : str
    1008        The prefix to match the files against.
    1009    qc : str
    1010        The quark combination extension to match the files against.
    1011    corr : str
    1012        The correlator to extract data for.
    1013    sep : str, optional
    1014        The separator to use when parsing the replika names.
    1015    **kwargs
    1016        Additional keyword arguments. The following keyword arguments are recognized:
    1017
    1018        - names (List[str]): A list of names to use for the replicas.
    1019
    1020    Returns
    1021    -------
    1022    Corr
    1023        A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
    1024    or
    1025    CObs
    1026        A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
    1027
    1028
    1029    Raises
    1030    ------
    1031    FileNotFoundError
    1032        If no files matching the specified prefix and quark combination extension are found in the specified directory.
    1033    IOError
    1034        If there is an error reading a file.
    1035    struct.error
    1036        If there is an error unpacking binary data.
    1037    """
    1038
    1039    # found = []
    1040    files = []
    1041    names = []
    1042
    1043    # test if the input is correct
    1044    if qc not in ['dd', 'ud', 'du', 'uu']:
    1045        raise Exception("Unknown quark conbination!")
    1046
    1047    if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
    1048        raise Exception("Unknown correlator!")
    1049
    1050    if "files" in kwargs:
    1051        known_files = kwargs.get("files")
    1052    else:
    1053        known_files = []
    1054    files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
    1055
    1056    if "names" in kwargs:
    1057        names = kwargs.get("names")
    1058    else:
    1059        for f in files:
    1060            if not sep == "":
    1061                se = f.split(".")[0]
    1062                for s in f.split(".")[1:-2]:
    1063                    se += "." + s
    1064                names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
    1065            else:
    1066                names.append(prefix)
    1067
    1068    names = sorted(names)
    1069    files = sorted(files)
    1070
    1071    cnfgs = []
    1072    realsamples = []
    1073    imagsamples = []
    1074    repnum = 0
    1075    for file in files:
    1076        with open(path + "/" + file, "rb") as fp:
    1077
    1078            t = fp.read(8)
    1079            kappa = struct.unpack('d', t)[0]
    1080            t = fp.read(8)
    1081            csw = struct.unpack('d', t)[0]
    1082            t = fp.read(8)
    1083            dF = struct.unpack('d', t)[0]
    1084            t = fp.read(8)
    1085            zF = struct.unpack('d', t)[0]
    1086
    1087            t = fp.read(4)
    1088            tmax = struct.unpack('i', t)[0]
    1089            t = fp.read(4)
    1090            bnd = struct.unpack('i', t)[0]
    1091
    1092            placesBI = ["gS", "gP",
    1093                        "gA", "gV",
    1094                        "gVt", "lA",
    1095                        "lV", "lVt",
    1096                        "lT", "lTt"]
    1097            placesBB = ["g1", "l1"]
    1098
    1099            # the chunks have the following structure:
    1100            # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
    1101
    1102            chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
    1103            packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
    1104            cnfgs.append([])
    1105            realsamples.append([])
    1106            imagsamples.append([])
    1107            for t in range(tmax):
    1108                realsamples[repnum].append([])
    1109                imagsamples[repnum].append([])
    1110
    1111            while True:
    1112                cnfgt = fp.read(chunksize)
    1113                if not cnfgt:
    1114                    break
    1115                asascii = struct.unpack(packstr, cnfgt)
    1116                cnfg = asascii[0]
    1117                cnfgs[repnum].append(cnfg)
    1118
    1119                if corr not in placesBB:
    1120                    tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
    1121                else:
    1122                    tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
    1123
    1124                corrres = [[], []]
    1125                for i in range(len(tmpcorr)):
    1126                    corrres[i % 2].append(tmpcorr[i])
    1127                for t in range(int(len(tmpcorr) / 2)):
    1128                    realsamples[repnum][t].append(corrres[0][t])
    1129                for t in range(int(len(tmpcorr) / 2)):
    1130                    imagsamples[repnum][t].append(corrres[1][t])
    1131        repnum += 1
    1132
    1133    s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
    1134    for rep in range(1, repnum):
    1135        s += ", " + str(len(realsamples[rep][t]))
    1136    s += " samples"
    1137    print(s)
    1138    print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
    1139
    1140    # we have the data now... but we need to re format the whole thing and put it into Corr objects.
    1141
    1142    compObs = []
    1143
    1144    for t in range(int(len(tmpcorr) / 2)):
    1145        compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
    1146                            Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
    1147
    1148    if len(compObs) == 1:
    1149        return compObs[0]
    1150    else:
    1151        return Corr(compObs)

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

###### Parameters

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

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

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

###### Returns

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

###### Raises

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