pyerrors.input.openQCD

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