pyerrors.input.openQCD

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

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

Parameters
  • path (str): path that contains the data files
  • prefix (str): all files in path that start with prefix are considered as input files. May be used together postfix to consider only special file endings. Prefix is ignored, if the keyword 'files' is used.
  • version (str): version of openQCD, default 2.0
  • names (list): list of names that is assigned to the data according according to the order in the file list. Use careful, if you do not provide file names!
  • r_start (list): list which contains the first config to be read for each replicum
  • r_stop (list): list which contains the last config to be read for each replicum
  • r_step (int): integer that defines a fixed step size between two measurements (in units of configs) If not given, r_step=1 is assumed.
  • postfix (str): postfix of the file to read, e.g. '.ms1' for openQCD-files
  • files (list): list which contains the filenames to be read. No automatic detection of files performed if given.
  • print_err (bool): Print additional information that is useful for debugging.
Returns
  • rwms (Obs): Reweighting factors read
def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
240def extract_t0(path, prefix, dtr_read, xmin, spatial_extent, fit_range=5, **kwargs):
241    """Extract t0 from given .ms.dat files. Returns t0 as Obs.
242
243    It is assumed that all boundary effects have
244    sufficiently decayed at x0=xmin.
245    The data around the zero crossing of t^2<E> - 0.3
246    is fitted with a linear function
247    from which the exact root is extracted.
248
249    It is assumed that one measurement is performed for each config.
250    If this is not the case, the resulting idl, as well as the handling
251    of r_start, r_stop and r_step is wrong and the user has to correct
252    this in the resulting observable.
253
254    Parameters
255    ----------
256    path : str
257        Path to .ms.dat files
258    prefix : str
259        Ensemble prefix
260    dtr_read : int
261        Determines how many trajectories should be skipped
262        when reading the ms.dat files.
263        Corresponds to dtr_cnfg / dtr_ms in the openQCD input file.
264    xmin : int
265        First timeslice where the boundary
266        effects have sufficiently decayed.
267    spatial_extent : int
268        spatial extent of the lattice, required for normalization.
269    fit_range : int
270        Number of data points left and right of the zero
271        crossing to be included in the linear fit. (Default: 5)
272    r_start : list
273        list which contains the first config to be read for each replicum.
274    r_stop : list
275        list which contains the last config to be read for each replicum.
276    r_step : int
277        integer that defines a fixed step size between two measurements (in units of configs)
278        If not given, r_step=1 is assumed.
279    plaquette : bool
280        If true extract the plaquette estimate of t0 instead.
281    names : list
282        list of names that is assigned to the data according according
283        to the order in the file list. Use careful, if you do not provide file names!
284    files : list
285        list which contains the filenames to be read. No automatic detection of
286        files performed if given.
287    plot_fit : bool
288        If true, the fit for the extraction of t0 is shown together with the data.
289    assume_thermalization : bool
290        If True: If the first record divided by the distance between two measurements is larger than
291        1, it is assumed that this is due to thermalization and the first measurement belongs
292        to the first config (default).
293        If False: The config numbers are assumed to be traj_number // difference
294
295    Returns
296    -------
297    t0 : Obs
298        Extracted t0
299    """
300
301    ls = []
302    for (dirpath, dirnames, filenames) in os.walk(path):
303        ls.extend(filenames)
304        break
305
306    if not ls:
307        raise Exception('Error, directory not found')
308
309    if 'files' in kwargs:
310        ls = kwargs.get('files')
311    else:
312        for exc in ls:
313            if not fnmatch.fnmatch(exc, prefix + '*.ms.dat'):
314                ls = list(set(ls) - set([exc]))
315        if len(ls) > 1:
316            ls.sort(key=lambda x: int(re.findall(r'\d+', x[len(prefix):])[0]))
317    replica = len(ls)
318
319    if 'r_start' in kwargs:
320        r_start = kwargs.get('r_start')
321        if len(r_start) != replica:
322            raise Exception('r_start does not match number of replicas')
323        r_start = [o if o else None for o in r_start]
324    else:
325        r_start = [None] * replica
326
327    if 'r_stop' in kwargs:
328        r_stop = kwargs.get('r_stop')
329        if len(r_stop) != replica:
330            raise Exception('r_stop does not match number of replicas')
331    else:
332        r_stop = [None] * replica
333
334    if 'r_step' in kwargs:
335        r_step = kwargs.get('r_step')
336    else:
337        r_step = 1
338
339    print('Extract t0 from', prefix, ',', replica, 'replica')
340
341    if 'names' in kwargs:
342        rep_names = kwargs.get('names')
343    else:
344        rep_names = []
345        for entry in ls:
346            truncated_entry = entry.split('.')[0]
347            idx = truncated_entry.index('r')
348            rep_names.append(truncated_entry[:idx] + '|' + truncated_entry[idx:])
349
350    Ysum = []
351
352    configlist = []
353    r_start_index = []
354    r_stop_index = []
355
356    for rep in range(replica):
357
358        with open(path + '/' + ls[rep], 'rb') as fp:
359            t = fp.read(12)
360            header = struct.unpack('iii', t)
361            if rep == 0:
362                dn = header[0]
363                nn = header[1]
364                tmax = header[2]
365            elif dn != header[0] or nn != header[1] or tmax != header[2]:
366                raise Exception('Replica parameters do not match.')
367
368            t = fp.read(8)
369            if rep == 0:
370                eps = struct.unpack('d', t)[0]
371                print('Step size:', eps, ', Maximal t value:', dn * (nn) * eps)
372            elif eps != struct.unpack('d', t)[0]:
373                raise Exception('Values for eps do not match among replica.')
374
375            Ysl = []
376
377            configlist.append([])
378            while True:
379                t = fp.read(4)
380                if (len(t) < 4):
381                    break
382                nc = struct.unpack('i', t)[0]
383                configlist[-1].append(nc)
384
385                t = fp.read(8 * tmax * (nn + 1))
386                if kwargs.get('plaquette'):
387                    if nc % dtr_read == 0:
388                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
389                t = fp.read(8 * tmax * (nn + 1))
390                if not kwargs.get('plaquette'):
391                    if nc % dtr_read == 0:
392                        Ysl.append(struct.unpack('d' * tmax * (nn + 1), t))
393                t = fp.read(8 * tmax * (nn + 1))
394
395        Ysum.append([])
396        for i, item in enumerate(Ysl):
397            Ysum[-1].append([np.mean(item[current + xmin:
398                             current + tmax - xmin])
399                            for current in range(0, len(item), tmax)])
400
401        diffmeas = configlist[-1][-1] - configlist[-1][-2]
402        configlist[-1] = [item // diffmeas for item in configlist[-1]]
403        if kwargs.get('assume_thermalization', True) and configlist[-1][0] > 1:
404            warnings.warn('Assume thermalization and that the first measurement belongs to the first config.')
405            offset = configlist[-1][0] - 1
406            configlist[-1] = [item - offset for item in configlist[-1]]
407
408        if r_start[rep] is None:
409            r_start_index.append(0)
410        else:
411            try:
412                r_start_index.append(configlist[-1].index(r_start[rep]))
413            except ValueError:
414                raise Exception('Config %d not in file with range [%d, %d]' % (
415                    r_start[rep], configlist[-1][0], configlist[-1][-1])) from None
416
417        if r_stop[rep] is None:
418            r_stop_index.append(len(configlist[-1]) - 1)
419        else:
420            try:
421                r_stop_index.append(configlist[-1].index(r_stop[rep]))
422            except ValueError:
423                raise Exception('Config %d not in file with range [%d, %d]' % (
424                    r_stop[rep], configlist[-1][0], configlist[-1][-1])) from None
425
426    if np.any([len(np.unique(np.diff(cl))) != 1 for cl in configlist]):
427        raise Exception('Irregular spaced data in input file!', [len(np.unique(np.diff(cl))) for cl in configlist])
428    stepsizes = [list(np.unique(np.diff(cl)))[0] for cl in configlist]
429    if np.any([step != 1 for step in stepsizes]):
430        warnings.warn('Stepsize between configurations is greater than one!' + str(stepsizes), RuntimeWarning)
431
432    idl = [range(configlist[rep][r_start_index[rep]], configlist[rep][r_stop_index[rep]] + 1, r_step) for rep in range(replica)]
433    t2E_dict = {}
434    for n in range(nn + 1):
435        samples = []
436        for nrep, rep in enumerate(Ysum):
437            samples.append([])
438            for cnfg in rep:
439                samples[-1].append(cnfg[n])
440            samples[-1] = samples[-1][r_start_index[nrep]:r_stop_index[nrep] + 1][::r_step]
441        new_obs = Obs(samples, rep_names, idl=idl)
442        t2E_dict[n * dn * eps] = (n * dn * eps) ** 2 * new_obs / (spatial_extent ** 3) - 0.3
443
444    zero_crossing = np.argmax(np.array(
445        [o.value for o in t2E_dict.values()]) > 0.0)
446
447    x = list(t2E_dict.keys())[zero_crossing - fit_range:
448                              zero_crossing + fit_range]
449    y = list(t2E_dict.values())[zero_crossing - fit_range:
450                                zero_crossing + fit_range]
451    [o.gamma_method() for o in y]
452
453    fit_result = fit_lin(x, y)
454
455    if kwargs.get('plot_fit'):
456        plt.figure()
457        gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
458        ax0 = plt.subplot(gs[0])
459        xmore = list(t2E_dict.keys())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
460        ymore = list(t2E_dict.values())[zero_crossing - fit_range - 2: zero_crossing + fit_range + 2]
461        [o.gamma_method() for o in ymore]
462        ax0.errorbar(xmore, [yi.value for yi in ymore], yerr=[yi.dvalue for yi in ymore], fmt='x')
463        xplot = np.linspace(np.min(x), np.max(x))
464        yplot = [fit_result[0] + fit_result[1] * xi for xi in xplot]
465        [yi.gamma_method() for yi in yplot]
466        ax0.fill_between(xplot, y1=[yi.value - yi.dvalue for yi in yplot], y2=[yi.value + yi.dvalue for yi in yplot])
467        retval = (-fit_result[0] / fit_result[1])
468        retval.gamma_method()
469        ylim = ax0.get_ylim()
470        ax0.fill_betweenx(ylim, x1=retval.value - retval.dvalue, x2=retval.value + retval.dvalue, color='gray', alpha=0.4)
471        ax0.set_ylim(ylim)
472        ax0.set_ylabel(r'$t^2 \langle E(t) \rangle - 0.3 $')
473        xlim = ax0.get_xlim()
474
475        fit_res = [fit_result[0] + fit_result[1] * xi for xi in x]
476        residuals = (np.asarray([o.value for o in y]) - [o.value for o in fit_res]) / np.asarray([o.dvalue for o in y])
477        ax1 = plt.subplot(gs[1])
478        ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
479        ax1.tick_params(direction='out')
480        ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
481        ax1.axhline(y=0.0, ls='--', color='k')
482        ax1.fill_between(xlim, -1.0, 1.0, alpha=0.1, facecolor='k')
483        ax1.set_xlim(xlim)
484        ax1.set_ylabel('Residuals')
485        ax1.set_xlabel(r'$t/a^2$')
486
487        plt.draw()
488    return -fit_result[0] / fit_result[1]

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

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

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

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