pyerrors.input.openQCD

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