pyerrors.input.openQCD

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

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

Parameters
  • path (str): path that contains the data files
  • prefix (str): all files in path that start with prefix are considered as input files. May be used together postfix to consider only special file endings. Prefix is ignored, if the keyword 'files' is used.
  • version (str): version of openQCD, default 2.0
  • names (list): list of names that is assigned to the data according according to the order in the file list. Use careful, if you do not provide file names!
  • r_start (list): list which contains the first config to be read for each replicum
  • r_stop (list): list which contains the last config to be read for each replicum
  • r_step (int): integer that defines a fixed step size between two measurements (in units of configs) If not given, r_step=1 is assumed.
  • postfix (str): postfix of the file to read, e.g. '.ms1' for openQCD-files
  • files (list): list which contains the filenames to be read. No automatic detection of files performed if given.
  • print_err (bool): Print additional information that is useful for debugging.
Returns
  • rwms (Obs): Reweighting factors read
def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
235def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
236    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
237
238    It is assumed that all boundary effects have
239    sufficiently decayed at x0=xmin.
240    The data around the zero crossing of t^2<E> - 0.3
241    is fitted with a linear function
242    from which the exact root is extracted.
243
244    It is assumed that one measurement is performed for each config.
245    If this is not the case, the resulting idl, as well as the handling
246    of r_start, r_stop and r_step is wrong and the user has to correct
247    this in the resulting observable.
248
249    Parameters
250    ----------
251    path : str
252        Path to .ms.dat files
253    prefix : str
254        Ensemble prefix
255    dtr_read : int
256        Determines how many trajectories should be skipped
257        when reading the ms.dat files.
258        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
259    xmin : int
260        First timeslice where the boundary
261        effects have sufficiently decayed.
262    spatial_extent : int
263        spatial extent of the lattice, required for normalization.
264    fit_range : int
265        Number of data points left and right of the zero
266        crossing to be included in the linear fit. (Default: 5)
267    r_start : list
268        list which contains the first config to be read for each replicum.
269    r_stop : list
270        list which contains the last config to be read for each replicum.
271    r_step : int
272        integer that defines a fixed step size between two measurements (in units of configs)
273        If not given, r_step=1 is assumed.
274    plaquette : bool
275        If true extract the plaquette estimate of t0 instead.
276    names : list
277        list of names that is assigned to the data according according
278        to the order in the file list. Use careful, if you do not provide file names!
279    files : list
280        list which contains the filenames to be read. No automatic detection of
281        files performed if given.
282    plot_fit : bool
283        If true, the fit for the extraction of t0 is shown together with the data.
284    assume_thermalization : bool
285        If True: If the first record divided by the distance between two measurements is larger than
286        1, it is assumed that this is due to thermalization and the first measurement belongs
287        to the first config (default).
288        If False: The config numbers are assumed to be traj_number // difference
289
290    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, 'ms', '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    zero_crossing = np.argmax(np.array(
431        [o.value for o in t2E_dict.values()]) > 0.0)
432
433    x = list(t2E_dict.keys())[zero_crossing - fit_range:
434                              zero_crossing + fit_range]
435    y = list(t2E_dict.values())[zero_crossing - fit_range:
436                                zero_crossing + fit_range]
437    [o.gamma_method() for o in y]
438
439    fit_result = fit_lin(x, y)
440
441    if kwargs.get('plot_fit'):
442        plt.figure()
443        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
444        ax0 = plt.subplot(gs[0])
445        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
446        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
447        [o.gamma_method() for o in ymore]
448        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
449        xplot = np.linspace(np.min(x), np.max(x))
450        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
451        [yi.gamma_method() for yi in yplot]
452        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
453        retval = (-fit_result[0] / fit_result[1])
454        retval.gamma_method()
455        ylim = ax0.get_ylim()
456        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
457        ax0.set_ylim(ylim)
458        ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
459        xlim = ax0.get_xlim()
460
461        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
462        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
463        ax1 = plt.subplot(gs[1])
464        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
465        ax1.tick_params(direction='out')
466        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
467        ax1.axhline(y=0.0, ls='--', color='k')
468        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
469        ax1.set_xlim(xlim)
470        ax1.set_ylabel('Residuals')
471        ax1.set_xlabel(r'$t/a^2$')
472
473        plt.draw()
474    return -fit_result[0] / fit_result[1]

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

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

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

Parameters
  • path (str): Path to .ms.dat files
  • prefix (str): Ensemble prefix
  • dtr_read (int): Determines how many trajectories should be skipped when reading the ms.dat files. Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
  • xmin (int): First timeslice where the boundary effects have sufficiently decayed.
  • spatial_extent (int): spatial extent of the lattice, required for normalization.
  • fit_range (int): Number of data points left and right of the zero crossing to be included in the linear fit. (Default: 5)
  • r_start (list): list which contains the first config to be read for each replicum.
  • r_stop (list): list which contains the last config to be read for each replicum.
  • r_step (int): integer that defines a fixed step size between two measurements (in units of configs) If not given, r_step=1 is assumed.
  • plaquette (bool): If true extract the plaquette estimate of t0 instead.
  • names (list): list of names that is assigned to the data according according to the order in the file list. Use careful, if you do not provide file names!
  • files (list): list which contains the filenames to be read. No automatic detection of files performed if given.
  • plot_fit (bool): If true, the fit for the extraction of t0 is shown together with the data.
  • assume_thermalization (bool): If True: If the first record divided by the distance between two measurements is larger than 1, it is assumed that this is due to thermalization and the first measurement belongs to the first config (default). If False: The config numbers are assumed to be traj_number // difference
Returns
  • t0 (Obs): Extracted t0
def read_qtop(path, prefix, c, dtr_cnfg=1, version='openQCD', **kwargs):
562def read_qtop(path, prefix, c, dtr_cnfg=1, version="openQCD", **kwargs):
563    """Read the topologial charge based on openQCD gradient flow measurements.
564
565    Parameters
566    ----------
567    path : str
568        path of the measurement files
569    prefix : str
570        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
571        Ignored if file names are passed explicitly via keyword files.
572    c : double
573        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
574    dtr_cnfg : int
575        (optional) parameter that specifies the number of measurements
576        between two configs.
577        If it is not set, the distance between two measurements
578        in the file is assumed to be the distance between two configurations.
579    steps : int
580        (optional) Distance between two configurations in units of trajectories /
581         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
582    version : str
583        Either openQCD or sfqcd, depending on the data.
584    L : int
585        spatial length of the lattice in L/a.
586        HAS to be set if version != sfqcd, since openQCD does not provide
587        this in the header
588    r_start : list
589        list which contains the first config to be read for each replicum.
590    r_stop : list
591        list which contains the last config to be read for each replicum.
592    files : list
593        specify the exact files that need to be read
594        from path, practical if e.g. only one replicum is needed
595    postfix : str
596        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
597    names : list
598        Alternative labeling for replicas/ensembles.
599        Has to have the appropriate length.
600    Zeuthen_flow : bool
601        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
602        for version=='sfqcd' If False, the Wilson flow is used.
603    integer_charge : bool
604        If True, the charge is rounded towards the nearest integer on each config.
605
606    Returns
607    -------
608    result : Obs
609        Read topological charge
610    """
611
612    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):
615def read_gf_coupling(path, prefix, c, dtr_cnfg=1, Zeuthen_flow=True, **kwargs):
616    """Read the gradient flow coupling based on sfqcd gradient flow measurements. See 1607.06423 for details.
617
618    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.
619
620    Parameters
621    ----------
622    path : str
623        path of the measurement files
624    prefix : str
625        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat.
626        Ignored if file names are passed explicitly via keyword files.
627    c : double
628        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L.
629    dtr_cnfg : int
630        (optional) parameter that specifies the number of measurements
631        between two configs.
632        If it is not set, the distance between two measurements
633        in the file is assumed to be the distance between two configurations.
634    steps : int
635        (optional) Distance between two configurations in units of trajectories /
636         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
637    r_start : list
638        list which contains the first config to be read for each replicum.
639    r_stop : list
640        list which contains the last config to be read for each replicum.
641    files : list
642        specify the exact files that need to be read
643        from path, practical if e.g. only one replicum is needed
644    names : list
645        Alternative labeling for replicas/ensembles.
646        Has to have the appropriate length.
647    postfix : str
648        postfix of the file to read, e.g. '.gfms.dat' for openQCD-files
649    Zeuthen_flow : bool
650        (optional) If True, the Zeuthen flow is used for the coupling. If False, the Wilson flow is used.
651    """
652
653    if c != 0.3:
654        raise Exception("The required lattice norm is only implemented for c=0.3 at the moment.")
655
656    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)
657    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)
658    L = plaq.tag["L"]
659    T = plaq.tag["T"]
660
661    if T != L:
662        raise Exception("The required lattice norm is only implemented for T=L at the moment.")
663
664    if Zeuthen_flow is not True:
665        raise Exception("The required lattice norm is only implemented for the Zeuthen flow at the moment.")
666
667    t = (c * L) ** 2 / 8
668
669    normdict = {4: 0.012341170468270,
670                6: 0.010162691462430,
671                8: 0.009031614807931,
672                10: 0.008744966371393,
673                12: 0.008650917856809,
674                14: 8.611154391267955E-03,
675                16: 0.008591758449508,
676                20: 0.008575359627103,
677                24: 0.008569387847540,
678                28: 8.566803713382559E-03,
679                32: 0.008565541650006,
680                40: 8.564480684962046E-03,
681                48: 8.564098025073460E-03,
682                64: 8.563853943383087E-03}
683
684    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):
959def qtop_projection(qtop, target=0):
960    """Returns the projection to the topological charge sector defined by target.
961
962    Parameters
963    ----------
964    path : Obs
965        Topological charge.
966    target : int
967        Specifies the topological sector to be reweighted to (default 0)
968
969    Returns
970    -------
971    reto : Obs
972        projection to the topological charge sector defined by target
973    """
974    if qtop.reweighted:
975        raise Exception('You can not use a reweighted observable for reweighting!')
976
977    proj_qtop = []
978    for n in qtop.deltas:
979        proj_qtop.append(np.array([1 if round(qtop.r_values[n] + q) == target else 0 for q in qtop.deltas[n]]))
980
981    reto = Obs(proj_qtop, qtop.names, idl=[qtop.idl[name] for name in qtop.names])
982    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):
 985def read_qtop_sector(path, prefix, c, target=0, **kwargs):
 986    """Constructs reweighting factors to a specified topological sector.
 987
 988    Parameters
 989    ----------
 990    path : str
 991        path of the measurement files
 992    prefix : str
 993        prefix of the measurement files, e.g. <prefix>_id0_r0.ms.dat
 994    c : double
 995        Smearing radius in units of the lattice extent, c = sqrt(8 t0) / L
 996    target : int
 997        Specifies the topological sector to be reweighted to (default 0)
 998    dtr_cnfg : int
 999        (optional) parameter that specifies the number of trajectories
1000        between two configs.
1001        if it is not set, the distance between two measurements
1002        in the file is assumed to be the distance between two configurations.
1003    steps : int
1004        (optional) Distance between two configurations in units of trajectories /
1005         cycles. Assumed to be the distance between two measurements * dtr_cnfg if not given
1006    version : str
1007        version string of the openQCD (sfqcd) version used to create
1008        the ensemble. Default is 2.0. May also be set to sfqcd.
1009    L : int
1010        spatial length of the lattice in L/a.
1011        HAS to be set if version != sfqcd, since openQCD does not provide
1012        this in the header
1013    r_start : list
1014        offset of the first ensemble, making it easier to match
1015        later on with other Obs
1016    r_stop : list
1017        last configurations that need to be read (per replicum)
1018    files : list
1019        specify the exact files that need to be read
1020        from path, practical if e.g. only one replicum is needed
1021    names : list
1022        Alternative labeling for replicas/ensembles.
1023        Has to have the appropriate length
1024    Zeuthen_flow : bool
1025        (optional) If True, the Zeuthen flow is used for Qtop. Only possible
1026        for version=='sfqcd' If False, the Wilson flow is used.
1027
1028    Returns
1029    -------
1030    reto : Obs
1031        projection to the topological charge sector defined by target
1032    """
1033
1034    if not isinstance(target, int):
1035        raise Exception("'target' has to be an integer.")
1036
1037    kwargs['integer_charge'] = True
1038    qtop = read_qtop(path, prefix, c, **kwargs)
1039
1040    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):
1043def read_ms5_xsf(path, prefix, qc, corr, sep="r", **kwargs):
1044    """
1045    Read data from files in the specified directory with the specified prefix and quark combination extension, and return a `Corr` object containing the data.
1046
1047    Parameters
1048    ----------
1049    path : str
1050        The directory to search for the files in.
1051    prefix : str
1052        The prefix to match the files against.
1053    qc : str
1054        The quark combination extension to match the files against.
1055    corr : str
1056        The correlator to extract data for.
1057    sep : str, optional
1058        The separator to use when parsing the replika names.
1059    **kwargs
1060        Additional keyword arguments. The following keyword arguments are recognized:
1061
1062        - names (List[str]): A list of names to use for the replicas.
1063
1064    Returns
1065    -------
1066    Corr
1067        A complex valued `Corr` object containing the data read from the files. In case of boudary to bulk correlators.
1068    or
1069    CObs
1070        A complex valued `CObs` object containing the data read from the files. In case of boudary to boundary correlators.
1071
1072
1073    Raises
1074    ------
1075    FileNotFoundError
1076        If no files matching the specified prefix and quark combination extension are found in the specified directory.
1077    IOError
1078        If there is an error reading a file.
1079    struct.error
1080        If there is an error unpacking binary data.
1081    """
1082
1083    # found = []
1084    files = []
1085    names = []
1086
1087    # test if the input is correct
1088    if qc not in ['dd', 'ud', 'du', 'uu']:
1089        raise Exception("Unknown quark conbination!")
1090
1091    if corr not in ["gS", "gP", "gA", "gV", "gVt", "lA", "lV", "lVt", "lT", "lTt", "g1", "l1"]:
1092        raise Exception("Unknown correlator!")
1093
1094    if "files" in kwargs:
1095        known_files = kwargs.get("files")
1096    else:
1097        known_files = []
1098    files = _find_files(path, prefix, "ms5_xsf_" + qc, "dat", known_files=known_files)
1099
1100    if "names" in kwargs:
1101        names = kwargs.get("names")
1102    else:
1103        for f in files:
1104            if not sep == "":
1105                se = f.split(".")[0]
1106                for s in f.split(".")[1:-2]:
1107                    se += "." + s
1108                names.append(se.split(sep)[0] + "|r" + se.split(sep)[1])
1109            else:
1110                names.append(prefix)
1111
1112    names = sorted(names)
1113    files = sorted(files)
1114
1115    cnfgs = []
1116    realsamples = []
1117    imagsamples = []
1118    repnum = 0
1119    for file in files:
1120        with open(path + "/" + file, "rb") as fp:
1121
1122            t = fp.read(8)
1123            kappa = struct.unpack('d', t)[0]
1124            t = fp.read(8)
1125            csw = struct.unpack('d', t)[0]
1126            t = fp.read(8)
1127            dF = struct.unpack('d', t)[0]
1128            t = fp.read(8)
1129            zF = struct.unpack('d', t)[0]
1130
1131            t = fp.read(4)
1132            tmax = struct.unpack('i', t)[0]
1133            t = fp.read(4)
1134            bnd = struct.unpack('i', t)[0]
1135
1136            placesBI = ["gS", "gP",
1137                        "gA", "gV",
1138                        "gVt", "lA",
1139                        "lV", "lVt",
1140                        "lT", "lTt"]
1141            placesBB = ["g1", "l1"]
1142
1143            # the chunks have the following structure:
1144            # confignumber, 10x timedependent complex correlators as doubles, 2x timeindependent complex correlators as doubles
1145
1146            chunksize = 4 + (8 * 2 * tmax * 10) + (8 * 2 * 2)
1147            packstr = '=i' + ('d' * 2 * tmax * 10) + ('d' * 2 * 2)
1148            cnfgs.append([])
1149            realsamples.append([])
1150            imagsamples.append([])
1151            for t in range(tmax):
1152                realsamples[repnum].append([])
1153                imagsamples[repnum].append([])
1154
1155            while True:
1156                cnfgt = fp.read(chunksize)
1157                if not cnfgt:
1158                    break
1159                asascii = struct.unpack(packstr, cnfgt)
1160                cnfg = asascii[0]
1161                cnfgs[repnum].append(cnfg)
1162
1163                if corr not in placesBB:
1164                    tmpcorr = asascii[1 + 2 * tmax * placesBI.index(corr):1 + 2 * tmax * placesBI.index(corr) + 2 * tmax]
1165                else:
1166                    tmpcorr = asascii[1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr):1 + 2 * tmax * len(placesBI) + 2 * placesBB.index(corr) + 2]
1167
1168                corrres = [[], []]
1169                for i in range(len(tmpcorr)):
1170                    corrres[i % 2].append(tmpcorr[i])
1171                for t in range(int(len(tmpcorr) / 2)):
1172                    realsamples[repnum][t].append(corrres[0][t])
1173                for t in range(int(len(tmpcorr) / 2)):
1174                    imagsamples[repnum][t].append(corrres[1][t])
1175        repnum += 1
1176
1177    s = "Read correlator " + corr + " from " + str(repnum) + " replika with " + str(len(realsamples[0][t]))
1178    for rep in range(1, repnum):
1179        s += ", " + str(len(realsamples[rep][t]))
1180    s += " samples"
1181    print(s)
1182    print("Asserted run parameters:\n T:", tmax, "kappa:", kappa, "csw:", csw, "dF:", dF, "zF:", zF, "bnd:", bnd)
1183
1184    # we have the data now... but we need to re format the whole thing and put it into Corr objects.
1185
1186    compObs = []
1187
1188    for t in range(int(len(tmpcorr) / 2)):
1189        compObs.append(CObs(Obs([realsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs),
1190                            Obs([imagsamples[rep][t] for rep in range(repnum)], names=names, idl=cnfgs)))
1191
1192    if len(compObs) == 1:
1193        return compObs[0]
1194    else:
1195        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.