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