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

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

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

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

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