pyerrors.input.bdio

  0import ctypes
  1import hashlib
  2import autograd.numpy as np  # Thinly-wrapped numpy
  3from ..obs import Obs
  4
  5
  6def read_ADerrors(file_path, bdio_path='./libbdio.so', **kwargs):
  7    """ Extract generic MCMC data from a bdio file
  8
  9    read_ADerrors requires bdio to be compiled into a shared library. This can be achieved by
 10    adding the flag -fPIC to CC and changing the all target to
 11
 12    all:		bdio.o $(LIBDIR)
 13                gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o
 14                cp $(BUILDDIR)/libbdio.so $(LIBDIR)/
 15
 16    Parameters
 17    ----------
 18    file_path -- path to the bdio file
 19    bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so)
 20    """
 21    bdio = ctypes.cdll.LoadLibrary(bdio_path)
 22
 23    bdio_open = bdio.bdio_open
 24    bdio_open.restype = ctypes.c_void_p
 25
 26    bdio_close = bdio.bdio_close
 27    bdio_close.restype = ctypes.c_int
 28    bdio_close.argtypes = [ctypes.c_void_p]
 29
 30    bdio_seek_record = bdio.bdio_seek_record
 31    bdio_seek_record.restype = ctypes.c_int
 32    bdio_seek_record.argtypes = [ctypes.c_void_p]
 33
 34    bdio_get_rlen = bdio.bdio_get_rlen
 35    bdio_get_rlen.restype = ctypes.c_int
 36    bdio_get_rlen.argtypes = [ctypes.c_void_p]
 37
 38    bdio_get_ruinfo = bdio.bdio_get_ruinfo
 39    bdio_get_ruinfo.restype = ctypes.c_int
 40    bdio_get_ruinfo.argtypes = [ctypes.c_void_p]
 41
 42    bdio_read = bdio.bdio_read
 43    bdio_read.restype = ctypes.c_size_t
 44    bdio_read.argtypes = [ctypes.c_char_p, ctypes.c_size_t, ctypes.c_void_p]
 45
 46    bdio_read_f64 = bdio.bdio_read_f64
 47    bdio_read_f64.restype = ctypes.c_size_t
 48    bdio_read_f64.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
 49
 50    bdio_read_int32 = bdio.bdio_read_int32
 51    bdio_read_int32.restype = ctypes.c_size_t
 52    bdio_read_int32.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
 53
 54    b_path = file_path.encode('utf-8')
 55    read = 'r'
 56    b_read = read.encode('utf-8')
 57
 58    fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), None)
 59
 60    return_list = []
 61
 62    print('Reading of bdio file started')
 63    while 1 > 0:
 64        bdio_seek_record(fbdio)
 65        ruinfo = bdio_get_ruinfo(fbdio)
 66
 67        if ruinfo == 7:
 68            print('MD5sum found')  # For now we just ignore these entries and do not perform any checks on them
 69            continue
 70
 71        if ruinfo < 0:
 72            # EOF reached
 73            break
 74        bdio_get_rlen(fbdio)
 75
 76        def read_c_double():
 77            d_buf = ctypes.c_double
 78            pd_buf = d_buf()
 79            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
 80            bdio_read_f64(ppd_buf, ctypes.c_size_t(8), ctypes.c_void_p(fbdio))
 81            return pd_buf.value
 82
 83        mean = read_c_double()
 84        print('mean', mean)
 85
 86        def read_c_size_t():
 87            d_buf = ctypes.c_size_t
 88            pd_buf = d_buf()
 89            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
 90            bdio_read_int32(ppd_buf, ctypes.c_size_t(4), ctypes.c_void_p(fbdio))
 91            return pd_buf.value
 92
 93        neid = read_c_size_t()
 94        print('neid', neid)
 95
 96        ndata = []
 97        for index in range(neid):
 98            ndata.append(read_c_size_t())
 99        print('ndata', ndata)
100
101        nrep = []
102        for index in range(neid):
103            nrep.append(read_c_size_t())
104        print('nrep', nrep)
105
106        vrep = []
107        for index in range(neid):
108            vrep.append([])
109            for jndex in range(nrep[index]):
110                vrep[-1].append(read_c_size_t())
111        print('vrep', vrep)
112
113        ids = []
114        for index in range(neid):
115            ids.append(read_c_size_t())
116        print('ids', ids)
117
118        nt = []
119        for index in range(neid):
120            nt.append(read_c_size_t())
121        print('nt', nt)
122
123        zero = []
124        for index in range(neid):
125            zero.append(read_c_double())
126        print('zero', zero)
127
128        four = []
129        for index in range(neid):
130            four.append(read_c_double())
131        print('four', four)
132
133        d_buf = ctypes.c_double * np.sum(ndata)
134        pd_buf = d_buf()
135        ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
136        bdio_read_f64(ppd_buf, ctypes.c_size_t(8 * np.sum(ndata)), ctypes.c_void_p(fbdio))
137        delta = pd_buf[:]
138
139        samples = np.split(np.asarray(delta) + mean, np.cumsum([a for su in vrep for a in su])[:-1])
140        no_reps = [len(o) for o in vrep]
141        assert len(ids) == len(no_reps)
142        tmp_names = []
143        ens_length = max([len(str(o)) for o in ids])
144        for loc_id, reps in zip(ids, no_reps):
145            for index in range(reps):
146                missing_chars = ens_length - len(str(loc_id))
147                tmp_names.append(str(loc_id) + ' ' * missing_chars + '|r' + '{0:03d}'.format(index))
148
149        return_list.append(Obs(samples, tmp_names))
150
151    bdio_close(fbdio)
152    print()
153    print(len(return_list), 'observable(s) extracted.')
154    return return_list
155
156
157def write_ADerrors(obs_list, file_path, bdio_path='./libbdio.so', **kwargs):
158    """ Write Obs to a bdio file according to ADerrors conventions
159
160    read_mesons requires bdio to be compiled into a shared library. This can be achieved by
161    adding the flag -fPIC to CC and changing the all target to
162
163    all:		bdio.o $(LIBDIR)
164                gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o
165                cp $(BUILDDIR)/libbdio.so $(LIBDIR)/
166
167    Parameters
168    ----------
169    file_path -- path to the bdio file
170    bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so)
171    """
172
173    for obs in obs_list:
174        if not hasattr(obs, 'e_names'):
175            raise Exception('Run the gamma method first for all obs.')
176
177    bdio = ctypes.cdll.LoadLibrary(bdio_path)
178
179    bdio_open = bdio.bdio_open
180    bdio_open.restype = ctypes.c_void_p
181
182    bdio_close = bdio.bdio_close
183    bdio_close.restype = ctypes.c_int
184    bdio_close.argtypes = [ctypes.c_void_p]
185
186    bdio_start_record = bdio.bdio_start_record
187    bdio_start_record.restype = ctypes.c_int
188    bdio_start_record.argtypes = [ctypes.c_size_t, ctypes.c_size_t, ctypes.c_void_p]
189
190    bdio_flush_record = bdio.bdio_flush_record
191    bdio_flush_record.restype = ctypes.c_int
192    bdio_flush_record.argytpes = [ctypes.c_void_p]
193
194    bdio_write_f64 = bdio.bdio_write_f64
195    bdio_write_f64.restype = ctypes.c_size_t
196    bdio_write_f64.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
197
198    bdio_write_int32 = bdio.bdio_write_int32
199    bdio_write_int32.restype = ctypes.c_size_t
200    bdio_write_int32.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
201
202    b_path = file_path.encode('utf-8')
203    write = 'w'
204    b_write = write.encode('utf-8')
205    form = 'pyerrors ADerror export'
206    b_form = form.encode('utf-8')
207
208    fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_write), b_form)
209
210    for obs in obs_list:
211        # mean = obs.value
212        neid = len(obs.e_names)
213        vrep = [[obs.shape[o] for o in sl] for sl in list(obs.e_content.values())]
214        vrep_write = [item for sublist in vrep for item in sublist]
215        ndata = [np.sum(o) for o in vrep]
216        nrep = [len(o) for o in vrep]
217        print('ndata', ndata)
218        print('nrep', nrep)
219        print('vrep', vrep)
220        keys = list(obs.e_content.keys())
221        ids = []
222        for key in keys:
223            try:  # Try to convert key to integer
224                ids.append(int(key))
225            except Exception:  # If not possible construct a hash
226                ids.append(int(hashlib.sha256(key.encode('utf-8')).hexdigest(), 16) % 10 ** 8)
227        print('ids', ids)
228        nt = []
229        for e, e_name in enumerate(obs.e_names):
230
231            r_length = []
232            for r_name in obs.e_content[e_name]:
233                r_length.append(len(obs.deltas[r_name]))
234
235            # e_N = np.sum(r_length)
236            nt.append(max(r_length) // 2)
237        print('nt', nt)
238        zero = neid * [0.0]
239        four = neid * [4.0]
240        print('zero', zero)
241        print('four', four)
242        delta = np.concatenate([item for sublist in [[obs.deltas[o] for o in sl] for sl in list(obs.e_content.values())] for item in sublist])
243
244        bdio_start_record(0x00, 8, fbdio)
245
246        def write_c_double(double):
247            pd_buf = ctypes.c_double(double)
248            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
249            bdio_write_f64(ppd_buf, ctypes.c_size_t(8), ctypes.c_void_p(fbdio))
250
251        def write_c_size_t(int32):
252            pd_buf = ctypes.c_size_t(int32)
253            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
254            bdio_write_int32(ppd_buf, ctypes.c_size_t(4), ctypes.c_void_p(fbdio))
255
256        write_c_double(obs.value)
257        write_c_size_t(neid)
258
259        for element in ndata:
260            write_c_size_t(element)
261        for element in nrep:
262            write_c_size_t(element)
263        for element in vrep_write:
264            write_c_size_t(element)
265        for element in ids:
266            write_c_size_t(element)
267        for element in nt:
268            write_c_size_t(element)
269
270        for element in zero:
271            write_c_double(element)
272        for element in four:
273            write_c_double(element)
274
275        for element in delta:
276            write_c_double(element)
277
278    bdio_close(fbdio)
279    return 0
280
281
282def _get_kwd(string, key):
283    return (string.split(key, 1)[1]).split(" ", 1)[0]
284
285
286def _get_corr_name(string, key):
287    return (string.split(key, 1)[1]).split(' NDIM=', 1)[0]
288
289
290def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs):
291    """ Extract mesons data from a bdio file and return it as a dictionary
292
293    The dictionary can be accessed with a tuple consisting of (type, source_position, kappa1, kappa2)
294
295    read_mesons requires bdio to be compiled into a shared library. This can be achieved by
296    adding the flag -fPIC to CC and changing the all target to
297
298    all:		bdio.o $(LIBDIR)
299                gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o
300                cp $(BUILDDIR)/libbdio.so $(LIBDIR)/
301
302    Parameters
303    ----------
304    file_path : str
305        path to the bdio file
306    bdio_path : str
307        path to the shared bdio library libbdio.so (default ./libbdio.so)
308    start : int
309        The first configuration to be read (default 1)
310    stop : int
311        The last configuration to be read (default None)
312    step : int
313        Fixed step size between two measurements (default 1)
314    alternative_ensemble_name : str
315        Manually overwrite ensemble name
316    """
317
318    start = kwargs.get('start', 1)
319    stop = kwargs.get('stop', None)
320    step = kwargs.get('step', 1)
321
322    bdio = ctypes.cdll.LoadLibrary(bdio_path)
323
324    bdio_open = bdio.bdio_open
325    bdio_open.restype = ctypes.c_void_p
326
327    bdio_close = bdio.bdio_close
328    bdio_close.restype = ctypes.c_int
329    bdio_close.argtypes = [ctypes.c_void_p]
330
331    bdio_seek_record = bdio.bdio_seek_record
332    bdio_seek_record.restype = ctypes.c_int
333    bdio_seek_record.argtypes = [ctypes.c_void_p]
334
335    bdio_get_rlen = bdio.bdio_get_rlen
336    bdio_get_rlen.restype = ctypes.c_int
337    bdio_get_rlen.argtypes = [ctypes.c_void_p]
338
339    bdio_get_ruinfo = bdio.bdio_get_ruinfo
340    bdio_get_ruinfo.restype = ctypes.c_int
341    bdio_get_ruinfo.argtypes = [ctypes.c_void_p]
342
343    bdio_read = bdio.bdio_read
344    bdio_read.restype = ctypes.c_size_t
345    bdio_read.argtypes = [ctypes.c_char_p, ctypes.c_size_t, ctypes.c_void_p]
346
347    bdio_read_f64 = bdio.bdio_read_f64
348    bdio_read_f64.restype = ctypes.c_size_t
349    bdio_read_f64.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
350
351    b_path = file_path.encode('utf-8')
352    read = 'r'
353    b_read = read.encode('utf-8')
354    form = 'Generic Correlator Format 1.0'
355    b_form = form.encode('utf-8')
356
357    ensemble_name = ''
358    volume = []  # lattice volume
359    boundary_conditions = []
360    corr_name = []  # Contains correlator names
361    corr_type = []  # Contains correlator data type (important for reading out numerical data)
362    corr_props = []  # Contanis propagator types (Component of corr_kappa)
363    d0 = 0  # tvals
364    d1 = 0  # nnoise
365    prop_kappa = []  # Contains propagator kappas (Component of corr_kappa)
366    prop_source = []  # Contains propagator source positions
367    # Check noise type for multiple replica?
368    corr_no = -1
369    data = []
370    idl = []
371
372    fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), ctypes.c_char_p(b_form))
373
374    print('Reading of bdio file started')
375    while 1 > 0:
376        bdio_seek_record(fbdio)
377        ruinfo = bdio_get_ruinfo(fbdio)
378        if ruinfo < 0:
379            # EOF reached
380            break
381        rlen = bdio_get_rlen(fbdio)
382        if ruinfo == 5:
383            d_buf = ctypes.c_double * (2 + d0 * d1 * 2)
384            pd_buf = d_buf()
385            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
386            iread = bdio_read_f64(ppd_buf, ctypes.c_size_t(rlen), ctypes.c_void_p(fbdio))
387            if corr_type[corr_no] == 'complex':
388                tmp_mean = np.mean(np.asarray(np.split(np.asarray(pd_buf[2 + 2 * d1:-2 * d1:2]), d0 - 2)), axis=1)
389            else:
390                tmp_mean = np.mean(np.asarray(np.split(np.asarray(pd_buf[2 + d1:-d0 * d1 - d1]), d0 - 2)), axis=1)
391
392            data[corr_no].append(tmp_mean)
393            corr_no += 1
394        else:
395            alt_buf = ctypes.create_string_buffer(1024)
396            palt_buf = ctypes.c_char_p(ctypes.addressof(alt_buf))
397            iread = bdio_read(palt_buf, ctypes.c_size_t(rlen), ctypes.c_void_p(fbdio))
398            if rlen != iread:
399                print('Error')
400            for i, item in enumerate(alt_buf):
401                if item == b'\x00':
402                    alt_buf[i] = b' '
403            tmp_string = (alt_buf[:].decode("utf-8")).rstrip()
404            if ruinfo == 0:
405                ensemble_name = _get_kwd(tmp_string, 'ENSEMBLE=')
406                volume.append(int(_get_kwd(tmp_string, 'L0=')))
407                volume.append(int(_get_kwd(tmp_string, 'L1=')))
408                volume.append(int(_get_kwd(tmp_string, 'L2=')))
409                volume.append(int(_get_kwd(tmp_string, 'L3=')))
410                boundary_conditions.append(_get_kwd(tmp_string, 'BC0='))
411                boundary_conditions.append(_get_kwd(tmp_string, 'BC1='))
412                boundary_conditions.append(_get_kwd(tmp_string, 'BC2='))
413                boundary_conditions.append(_get_kwd(tmp_string, 'BC3='))
414
415            if ruinfo == 1:
416                corr_name.append(_get_corr_name(tmp_string, 'CORR_NAME='))
417                corr_type.append(_get_kwd(tmp_string, 'DATATYPE='))
418                corr_props.append([_get_kwd(tmp_string, 'PROP0='), _get_kwd(tmp_string, 'PROP1=')])
419                if d0 == 0:
420                    d0 = int(_get_kwd(tmp_string, 'D0='))
421                else:
422                    if d0 != int(_get_kwd(tmp_string, 'D0=')):
423                        print('Error: Varying number of time values')
424                if d1 == 0:
425                    d1 = int(_get_kwd(tmp_string, 'D1='))
426                else:
427                    if d1 != int(_get_kwd(tmp_string, 'D1=')):
428                        print('Error: Varying number of random sources')
429            if ruinfo == 2:
430                prop_kappa.append(_get_kwd(tmp_string, 'KAPPA='))
431                prop_source.append(_get_kwd(tmp_string, 'x0='))
432            if ruinfo == 4:
433                cnfg_no = int(_get_kwd(tmp_string, 'CNFG_ID='))
434                if stop:
435                    if cnfg_no > kwargs.get('stop'):
436                        break
437                idl.append(cnfg_no)
438                print('\r%s %i' % ('Reading configuration', cnfg_no), end='\r')
439                if len(idl) == 1:
440                    no_corrs = len(corr_name)
441                    data = []
442                    for c in range(no_corrs):
443                        data.append([])
444
445                corr_no = 0
446
447    bdio_close(fbdio)
448
449    print('\nEnsemble: ', ensemble_name)
450    if 'alternative_ensemble_name' in kwargs:
451        ensemble_name = kwargs.get('alternative_ensemble_name')
452        print('Ensemble name overwritten to', ensemble_name)
453    print('Lattice volume: ', volume)
454    print('Boundary conditions: ', boundary_conditions)
455    print('Number of time values: ', d0)
456    print('Number of random sources: ', d1)
457    print('Number of corrs: ', len(corr_name))
458    print('Number of configurations: ', len(idl))
459
460    corr_kappa = []  # Contains kappa values for both propagators of given correlation function
461    corr_source = []
462    for item in corr_props:
463        corr_kappa.append([float(prop_kappa[int(item[0])]), float(prop_kappa[int(item[1])])])
464        if prop_source[int(item[0])] != prop_source[int(item[1])]:
465            raise Exception('Source position do not match for correlator' + str(item))
466        else:
467            corr_source.append(int(prop_source[int(item[0])]))
468
469    if stop is None:
470        stop = idl[-1]
471    idl_target = range(start, stop + 1, step)
472
473    if set(idl) != set(idl_target):
474        try:
475            indices = [idl.index(i) for i in idl_target]
476        except ValueError as err:
477            raise Exception('Configurations in file do no match target list!', err)
478    else:
479        indices = None
480
481    result = {}
482    for c in range(no_corrs):
483        tmp_corr = []
484        tmp_data = np.asarray(data[c])
485        for t in range(d0 - 2):
486            if indices:
487                deltas = [tmp_data[:, t][index] for index in indices]
488            else:
489                deltas = tmp_data[:, t]
490            tmp_corr.append(Obs([deltas], [ensemble_name], idl=[idl_target]))
491        result[(corr_name[c], corr_source[c]) + tuple(corr_kappa[c])] = tmp_corr
492
493    # Check that all data entries have the same number of configurations
494    if len(set([o[0].N for o in list(result.values())])) != 1:
495        raise Exception('Error: Not all correlators have the same number of configurations. bdio file is possibly corrupted.')
496
497    return result
498
499
500def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs):
501    """ Extract dSdm data from a bdio file and return it as a dictionary
502
503    The dictionary can be accessed with a tuple consisting of (type, kappa)
504
505    read_dSdm requires bdio to be compiled into a shared library. This can be achieved by
506    adding the flag -fPIC to CC and changing the all target to
507
508    all:		bdio.o $(LIBDIR)
509                gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o
510                cp $(BUILDDIR)/libbdio.so $(LIBDIR)/
511
512    Parameters
513    ----------
514    file_path : str
515        path to the bdio file
516    bdio_path : str
517        path to the shared bdio library libbdio.so (default ./libbdio.so)
518    start : int
519        The first configuration to be read (default 1)
520    stop : int
521        The last configuration to be read (default None)
522    step : int
523        Fixed step size between two measurements (default 1)
524    alternative_ensemble_name : str
525        Manually overwrite ensemble name
526    """
527
528    start = kwargs.get('start', 1)
529    stop = kwargs.get('stop', None)
530    step = kwargs.get('step', 1)
531
532    bdio = ctypes.cdll.LoadLibrary(bdio_path)
533
534    bdio_open = bdio.bdio_open
535    bdio_open.restype = ctypes.c_void_p
536
537    bdio_close = bdio.bdio_close
538    bdio_close.restype = ctypes.c_int
539    bdio_close.argtypes = [ctypes.c_void_p]
540
541    bdio_seek_record = bdio.bdio_seek_record
542    bdio_seek_record.restype = ctypes.c_int
543    bdio_seek_record.argtypes = [ctypes.c_void_p]
544
545    bdio_get_rlen = bdio.bdio_get_rlen
546    bdio_get_rlen.restype = ctypes.c_int
547    bdio_get_rlen.argtypes = [ctypes.c_void_p]
548
549    bdio_get_ruinfo = bdio.bdio_get_ruinfo
550    bdio_get_ruinfo.restype = ctypes.c_int
551    bdio_get_ruinfo.argtypes = [ctypes.c_void_p]
552
553    bdio_read = bdio.bdio_read
554    bdio_read.restype = ctypes.c_size_t
555    bdio_read.argtypes = [ctypes.c_char_p, ctypes.c_size_t, ctypes.c_void_p]
556
557    bdio_read_f64 = bdio.bdio_read_f64
558    bdio_read_f64.restype = ctypes.c_size_t
559    bdio_read_f64.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
560
561    b_path = file_path.encode('utf-8')
562    read = 'r'
563    b_read = read.encode('utf-8')
564    form = 'Generic Correlator Format 1.0'
565    b_form = form.encode('utf-8')
566
567    ensemble_name = ''
568    volume = []  # lattice volume
569    boundary_conditions = []
570    corr_name = []  # Contains correlator names
571    corr_type = []  # Contains correlator data type (important for reading out numerical data)
572    corr_props = []  # Contains propagator types (Component of corr_kappa)
573    d0 = 0  # tvals
574    # d1 = 0  # nnoise
575    prop_kappa = []  # Contains propagator kappas (Component of corr_kappa)
576    # Check noise type for multiple replica?
577    corr_no = -1
578    data = []
579    idl = []
580
581    fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), ctypes.c_char_p(b_form))
582
583    print('Reading of bdio file started')
584    while 1 > 0:
585        bdio_seek_record(fbdio)
586        ruinfo = bdio_get_ruinfo(fbdio)
587        if ruinfo < 0:
588            # EOF reached
589            break
590        rlen = bdio_get_rlen(fbdio)
591        if ruinfo == 5:
592            d_buf = ctypes.c_double * (2 + d0)
593            pd_buf = d_buf()
594            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
595            iread = bdio_read_f64(ppd_buf, ctypes.c_size_t(rlen), ctypes.c_void_p(fbdio))
596            tmp_mean = np.mean(np.asarray(pd_buf[2:]))
597
598            data[corr_no].append(tmp_mean)
599            corr_no += 1
600        else:
601            alt_buf = ctypes.create_string_buffer(1024)
602            palt_buf = ctypes.c_char_p(ctypes.addressof(alt_buf))
603            iread = bdio_read(palt_buf, ctypes.c_size_t(rlen), ctypes.c_void_p(fbdio))
604            if rlen != iread:
605                print('Error')
606            for i, item in enumerate(alt_buf):
607                if item == b'\x00':
608                    alt_buf[i] = b' '
609            tmp_string = (alt_buf[:].decode("utf-8")).rstrip()
610            if ruinfo == 0:
611                creator = _get_kwd(tmp_string, 'CREATOR=')
612                ensemble_name = _get_kwd(tmp_string, 'ENSEMBLE=')
613                volume.append(int(_get_kwd(tmp_string, 'L0=')))
614                volume.append(int(_get_kwd(tmp_string, 'L1=')))
615                volume.append(int(_get_kwd(tmp_string, 'L2=')))
616                volume.append(int(_get_kwd(tmp_string, 'L3=')))
617                boundary_conditions.append(_get_kwd(tmp_string, 'BC0='))
618                boundary_conditions.append(_get_kwd(tmp_string, 'BC1='))
619                boundary_conditions.append(_get_kwd(tmp_string, 'BC2='))
620                boundary_conditions.append(_get_kwd(tmp_string, 'BC3='))
621
622            if ruinfo == 1:
623                corr_name.append(_get_corr_name(tmp_string, 'CORR_NAME='))
624                corr_type.append(_get_kwd(tmp_string, 'DATATYPE='))
625                corr_props.append(_get_kwd(tmp_string, 'PROP0='))
626                if d0 == 0:
627                    d0 = int(_get_kwd(tmp_string, 'D0='))
628                else:
629                    if d0 != int(_get_kwd(tmp_string, 'D0=')):
630                        print('Error: Varying number of time values')
631            if ruinfo == 2:
632                prop_kappa.append(_get_kwd(tmp_string, 'KAPPA='))
633            if ruinfo == 4:
634                cnfg_no = int(_get_kwd(tmp_string, 'CNFG_ID='))
635                if stop:
636                    if cnfg_no > kwargs.get('stop'):
637                        break
638                idl.append(cnfg_no)
639                print('\r%s %i' % ('Reading configuration', cnfg_no), end='\r')
640                if len(idl) == 1:
641                    no_corrs = len(corr_name)
642                    data = []
643                    for c in range(no_corrs):
644                        data.append([])
645
646                corr_no = 0
647    bdio_close(fbdio)
648
649    print('\nCreator: ', creator)
650    print('Ensemble: ', ensemble_name)
651    print('Lattice volume: ', volume)
652    print('Boundary conditions: ', boundary_conditions)
653    print('Number of random sources: ', d0)
654    print('Number of corrs: ', len(corr_name))
655    print('Number of configurations: ', cnfg_no + 1)
656
657    corr_kappa = []  # Contains kappa values for both propagators of given correlation function
658    for item in corr_props:
659        corr_kappa.append(float(prop_kappa[int(item)]))
660
661    if stop is None:
662        stop = idl[-1]
663    idl_target = range(start, stop + 1, step)
664    try:
665        indices = [idl.index(i) for i in idl_target]
666    except ValueError as err:
667        raise Exception('Configurations in file do no match target list!', err)
668
669    result = {}
670    for c in range(no_corrs):
671        deltas = [np.asarray(data[c])[index] for index in indices]
672        result[(corr_name[c], str(corr_kappa[c]))] = Obs([deltas], [ensemble_name], idl=[idl_target])
673
674    # Check that all data entries have the same number of configurations
675    if len(set([o.N for o in list(result.values())])) != 1:
676        raise Exception('Error: Not all correlators have the same number of configurations. bdio file is possibly corrupted.')
677
678    return result