pyerrors.input.bdio

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

Extract generic MCMC data from a bdio file

read_ADerrors requires bdio to be compiled into a shared library. This can be achieved by adding the flag -fPIC to CC and changing the all target to

all: bdio.o $(LIBDIR) gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o cp $(BUILDDIR)/libbdio.so $(LIBDIR)/

Parameters
  • file_path -- path to the bdio file
  • bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so)
Returns
  • data (List[Obs]): Extracted data
def write_ADerrors(obs_list, file_path, bdio_path='./libbdio.so', **kwargs):
164def write_ADerrors(obs_list, file_path, bdio_path='./libbdio.so', **kwargs):
165    """ Write Obs to a bdio file according to ADerrors conventions
166
167    read_mesons requires bdio to be compiled into a shared library. This can be achieved by
168    adding the flag -fPIC to CC and changing the all target to
169
170    all:		bdio.o $(LIBDIR)
171                gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o
172                cp $(BUILDDIR)/libbdio.so $(LIBDIR)/
173
174    Parameters
175    ----------
176    file_path -- path to the bdio file
177    bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so)
178
179    Returns
180    -------
181    success : int
182        returns 0 is successful
183    """
184
185    for obs in obs_list:
186        if not hasattr(obs, 'e_names'):
187            raise Exception('Run the gamma method first for all obs.')
188
189    bdio = ctypes.cdll.LoadLibrary(bdio_path)
190
191    bdio_open = bdio.bdio_open
192    bdio_open.restype = ctypes.c_void_p
193
194    bdio_close = bdio.bdio_close
195    bdio_close.restype = ctypes.c_int
196    bdio_close.argtypes = [ctypes.c_void_p]
197
198    bdio_start_record = bdio.bdio_start_record
199    bdio_start_record.restype = ctypes.c_int
200    bdio_start_record.argtypes = [ctypes.c_size_t, ctypes.c_size_t, ctypes.c_void_p]
201
202    bdio_flush_record = bdio.bdio_flush_record
203    bdio_flush_record.restype = ctypes.c_int
204    bdio_flush_record.argytpes = [ctypes.c_void_p]
205
206    bdio_write_f64 = bdio.bdio_write_f64
207    bdio_write_f64.restype = ctypes.c_size_t
208    bdio_write_f64.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
209
210    bdio_write_int32 = bdio.bdio_write_int32
211    bdio_write_int32.restype = ctypes.c_size_t
212    bdio_write_int32.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
213
214    b_path = file_path.encode('utf-8')
215    write = 'w'
216    b_write = write.encode('utf-8')
217    form = 'pyerrors ADerror export'
218    b_form = form.encode('utf-8')
219
220    fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_write), b_form)
221
222    for obs in obs_list:
223        # mean = obs.value
224        neid = len(obs.e_names)
225        vrep = [[obs.shape[o] for o in sl] for sl in list(obs.e_content.values())]
226        vrep_write = [item for sublist in vrep for item in sublist]
227        ndata = [np.sum(o) for o in vrep]
228        nrep = [len(o) for o in vrep]
229        print('ndata', ndata)
230        print('nrep', nrep)
231        print('vrep', vrep)
232        keys = list(obs.e_content.keys())
233        ids = []
234        for key in keys:
235            try:  # Try to convert key to integer
236                ids.append(int(key))
237            except Exception:  # If not possible construct a hash
238                ids.append(int(hashlib.sha256(key.encode('utf-8')).hexdigest(), 16) % 10 ** 8)
239        print('ids', ids)
240        nt = []
241        for e, e_name in enumerate(obs.e_names):
242
243            r_length = []
244            for r_name in obs.e_content[e_name]:
245                r_length.append(len(obs.deltas[r_name]))
246
247            # e_N = np.sum(r_length)
248            nt.append(max(r_length) // 2)
249        print('nt', nt)
250        zero = neid * [0.0]
251        four = neid * [4.0]
252        print('zero', zero)
253        print('four', four)
254        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])
255
256        bdio_start_record(0x00, 8, fbdio)
257
258        def write_c_double(double):
259            pd_buf = ctypes.c_double(double)
260            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
261            bdio_write_f64(ppd_buf, ctypes.c_size_t(8), ctypes.c_void_p(fbdio))
262
263        def write_c_size_t(int32):
264            pd_buf = ctypes.c_size_t(int32)
265            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
266            bdio_write_int32(ppd_buf, ctypes.c_size_t(4), ctypes.c_void_p(fbdio))
267
268        write_c_double(obs.value)
269        write_c_size_t(neid)
270
271        for element in ndata:
272            write_c_size_t(element)
273        for element in nrep:
274            write_c_size_t(element)
275        for element in vrep_write:
276            write_c_size_t(element)
277        for element in ids:
278            write_c_size_t(element)
279        for element in nt:
280            write_c_size_t(element)
281
282        for element in zero:
283            write_c_double(element)
284        for element in four:
285            write_c_double(element)
286
287        for element in delta:
288            write_c_double(element)
289
290    bdio_close(fbdio)
291    return 0

Write Obs to a bdio file according to ADerrors conventions

read_mesons requires bdio to be compiled into a shared library. This can be achieved by adding the flag -fPIC to CC and changing the all target to

all: bdio.o $(LIBDIR) gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o cp $(BUILDDIR)/libbdio.so $(LIBDIR)/

Parameters
  • file_path -- path to the bdio file
  • bdio_path -- path to the shared bdio library libbdio.so (default ./libbdio.so)
Returns
  • success (int): returns 0 is successful
def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs):
302def read_mesons(file_path, bdio_path='./libbdio.so', **kwargs):
303    """ Extract mesons data from a bdio file and return it as a dictionary
304
305    The dictionary can be accessed with a tuple consisting of (type, source_position, kappa1, kappa2)
306
307    read_mesons requires bdio to be compiled into a shared library. This can be achieved by
308    adding the flag -fPIC to CC and changing the all target to
309
310    all:		bdio.o $(LIBDIR)
311                gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o
312                cp $(BUILDDIR)/libbdio.so $(LIBDIR)/
313
314    Parameters
315    ----------
316    file_path : str
317        path to the bdio file
318    bdio_path : str
319        path to the shared bdio library libbdio.so (default ./libbdio.so)
320    start : int
321        The first configuration to be read (default 1)
322    stop : int
323        The last configuration to be read (default None)
324    step : int
325        Fixed step size between two measurements (default 1)
326    alternative_ensemble_name : str
327        Manually overwrite ensemble name
328
329    Returns
330    -------
331    data : dict
332        Extracted meson data
333    """
334
335    start = kwargs.get('start', 1)
336    stop = kwargs.get('stop', None)
337    step = kwargs.get('step', 1)
338
339    bdio = ctypes.cdll.LoadLibrary(bdio_path)
340
341    bdio_open = bdio.bdio_open
342    bdio_open.restype = ctypes.c_void_p
343
344    bdio_close = bdio.bdio_close
345    bdio_close.restype = ctypes.c_int
346    bdio_close.argtypes = [ctypes.c_void_p]
347
348    bdio_seek_record = bdio.bdio_seek_record
349    bdio_seek_record.restype = ctypes.c_int
350    bdio_seek_record.argtypes = [ctypes.c_void_p]
351
352    bdio_get_rlen = bdio.bdio_get_rlen
353    bdio_get_rlen.restype = ctypes.c_int
354    bdio_get_rlen.argtypes = [ctypes.c_void_p]
355
356    bdio_get_ruinfo = bdio.bdio_get_ruinfo
357    bdio_get_ruinfo.restype = ctypes.c_int
358    bdio_get_ruinfo.argtypes = [ctypes.c_void_p]
359
360    bdio_read = bdio.bdio_read
361    bdio_read.restype = ctypes.c_size_t
362    bdio_read.argtypes = [ctypes.c_char_p, ctypes.c_size_t, ctypes.c_void_p]
363
364    bdio_read_f64 = bdio.bdio_read_f64
365    bdio_read_f64.restype = ctypes.c_size_t
366    bdio_read_f64.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
367
368    b_path = file_path.encode('utf-8')
369    read = 'r'
370    b_read = read.encode('utf-8')
371    form = 'Generic Correlator Format 1.0'
372    b_form = form.encode('utf-8')
373
374    ensemble_name = ''
375    volume = []  # lattice volume
376    boundary_conditions = []
377    corr_name = []  # Contains correlator names
378    corr_type = []  # Contains correlator data type (important for reading out numerical data)
379    corr_props = []  # Contanis propagator types (Component of corr_kappa)
380    d0 = 0  # tvals
381    d1 = 0  # nnoise
382    prop_kappa = []  # Contains propagator kappas (Component of corr_kappa)
383    prop_source = []  # Contains propagator source positions
384    # Check noise type for multiple replica?
385    corr_no = -1
386    data = []
387    idl = []
388
389    fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), ctypes.c_char_p(b_form))
390
391    print('Reading of bdio file started')
392    while True:
393        bdio_seek_record(fbdio)
394        ruinfo = bdio_get_ruinfo(fbdio)
395        if ruinfo < 0:
396            # EOF reached
397            break
398        rlen = bdio_get_rlen(fbdio)
399        if ruinfo == 5:
400            d_buf = ctypes.c_double * (2 + d0 * d1 * 2)
401            pd_buf = d_buf()
402            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
403            bdio_read_f64(ppd_buf, ctypes.c_size_t(rlen), ctypes.c_void_p(fbdio))
404            if corr_type[corr_no] == 'complex':
405                tmp_mean = np.mean(np.asarray(np.split(np.asarray(pd_buf[2 + 2 * d1:-2 * d1:2]), d0 - 2)), axis=1)
406            else:
407                tmp_mean = np.mean(np.asarray(np.split(np.asarray(pd_buf[2 + d1:-d0 * d1 - d1]), d0 - 2)), axis=1)
408
409            data[corr_no].append(tmp_mean)
410            corr_no += 1
411        else:
412            alt_buf = ctypes.create_string_buffer(1024)
413            palt_buf = ctypes.c_char_p(ctypes.addressof(alt_buf))
414            iread = bdio_read(palt_buf, ctypes.c_size_t(rlen), ctypes.c_void_p(fbdio))
415            if rlen != iread:
416                print('Error')
417            for i, item in enumerate(alt_buf):
418                if item == b'\x00':
419                    alt_buf[i] = b' '
420            tmp_string = (alt_buf[:].decode("utf-8")).rstrip()
421            if ruinfo == 0:
422                ensemble_name = _get_kwd(tmp_string, 'ENSEMBLE=')
423                volume.append(int(_get_kwd(tmp_string, 'L0=')))
424                volume.append(int(_get_kwd(tmp_string, 'L1=')))
425                volume.append(int(_get_kwd(tmp_string, 'L2=')))
426                volume.append(int(_get_kwd(tmp_string, 'L3=')))
427                boundary_conditions.append(_get_kwd(tmp_string, 'BC0='))
428                boundary_conditions.append(_get_kwd(tmp_string, 'BC1='))
429                boundary_conditions.append(_get_kwd(tmp_string, 'BC2='))
430                boundary_conditions.append(_get_kwd(tmp_string, 'BC3='))
431
432            if ruinfo == 1:
433                corr_name.append(_get_corr_name(tmp_string, 'CORR_NAME='))
434                corr_type.append(_get_kwd(tmp_string, 'DATATYPE='))
435                corr_props.append([_get_kwd(tmp_string, 'PROP0='), _get_kwd(tmp_string, 'PROP1=')])
436                if d0 == 0:
437                    d0 = int(_get_kwd(tmp_string, 'D0='))
438                else:
439                    if d0 != int(_get_kwd(tmp_string, 'D0=')):
440                        print('Error: Varying number of time values')
441                if d1 == 0:
442                    d1 = int(_get_kwd(tmp_string, 'D1='))
443                else:
444                    if d1 != int(_get_kwd(tmp_string, 'D1=')):
445                        print('Error: Varying number of random sources')
446            if ruinfo == 2:
447                prop_kappa.append(_get_kwd(tmp_string, 'KAPPA='))
448                prop_source.append(_get_kwd(tmp_string, 'x0='))
449            if ruinfo == 4:
450                cnfg_no = int(_get_kwd(tmp_string, 'CNFG_ID='))
451                if stop:
452                    if cnfg_no > kwargs.get('stop'):
453                        break
454                idl.append(cnfg_no)
455                print('\r%s %i' % ('Reading configuration', cnfg_no), end='\r')
456                if len(idl) == 1:
457                    no_corrs = len(corr_name)
458                    data = []
459                    for c in range(no_corrs):
460                        data.append([])
461
462                corr_no = 0
463
464    bdio_close(fbdio)
465
466    print('\nEnsemble: ', ensemble_name)
467    if 'alternative_ensemble_name' in kwargs:
468        ensemble_name = kwargs.get('alternative_ensemble_name')
469        print('Ensemble name overwritten to', ensemble_name)
470    print('Lattice volume: ', volume)
471    print('Boundary conditions: ', boundary_conditions)
472    print('Number of time values: ', d0)
473    print('Number of random sources: ', d1)
474    print('Number of corrs: ', len(corr_name))
475    print('Number of configurations: ', len(idl))
476
477    corr_kappa = []  # Contains kappa values for both propagators of given correlation function
478    corr_source = []
479    for item in corr_props:
480        corr_kappa.append([float(prop_kappa[int(item[0])]), float(prop_kappa[int(item[1])])])
481        if prop_source[int(item[0])] != prop_source[int(item[1])]:
482            raise Exception('Source position do not match for correlator' + str(item))
483        else:
484            corr_source.append(int(prop_source[int(item[0])]))
485
486    if stop is None:
487        stop = idl[-1]
488    idl_target = range(start, stop + 1, step)
489
490    if set(idl) != set(idl_target):
491        try:
492            indices = [idl.index(i) for i in idl_target]
493        except ValueError as err:
494            raise Exception('Configurations in file do no match target list!', err)
495    else:
496        indices = None
497
498    result = {}
499    for c in range(no_corrs):
500        tmp_corr = []
501        tmp_data = np.asarray(data[c])
502        for t in range(d0 - 2):
503            if indices:
504                deltas = [tmp_data[:, t][index] for index in indices]
505            else:
506                deltas = tmp_data[:, t]
507            tmp_corr.append(Obs([deltas], [ensemble_name], idl=[idl_target]))
508        result[(corr_name[c], corr_source[c]) + tuple(corr_kappa[c])] = tmp_corr
509
510    # Check that all data entries have the same number of configurations
511    if len(set([o[0].N for o in list(result.values())])) != 1:
512        raise Exception('Error: Not all correlators have the same number of configurations. bdio file is possibly corrupted.')
513
514    return result

Extract mesons data from a bdio file and return it as a dictionary

The dictionary can be accessed with a tuple consisting of (type, source_position, kappa1, kappa2)

read_mesons requires bdio to be compiled into a shared library. This can be achieved by adding the flag -fPIC to CC and changing the all target to

all: bdio.o $(LIBDIR) gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o cp $(BUILDDIR)/libbdio.so $(LIBDIR)/

Parameters
  • file_path (str): path to the bdio file
  • bdio_path (str): path to the shared bdio library libbdio.so (default ./libbdio.so)
  • start (int): The first configuration to be read (default 1)
  • stop (int): The last configuration to be read (default None)
  • step (int): Fixed step size between two measurements (default 1)
  • alternative_ensemble_name (str): Manually overwrite ensemble name
Returns
  • data (dict): Extracted meson data
def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs):
517def read_dSdm(file_path, bdio_path='./libbdio.so', **kwargs):
518    """ Extract dSdm data from a bdio file and return it as a dictionary
519
520    The dictionary can be accessed with a tuple consisting of (type, kappa)
521
522    read_dSdm requires bdio to be compiled into a shared library. This can be achieved by
523    adding the flag -fPIC to CC and changing the all target to
524
525    all:		bdio.o $(LIBDIR)
526                gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o
527                cp $(BUILDDIR)/libbdio.so $(LIBDIR)/
528
529    Parameters
530    ----------
531    file_path : str
532        path to the bdio file
533    bdio_path : str
534        path to the shared bdio library libbdio.so (default ./libbdio.so)
535    start : int
536        The first configuration to be read (default 1)
537    stop : int
538        The last configuration to be read (default None)
539    step : int
540        Fixed step size between two measurements (default 1)
541    alternative_ensemble_name : str
542        Manually overwrite ensemble name
543    """
544
545    start = kwargs.get('start', 1)
546    stop = kwargs.get('stop', None)
547    step = kwargs.get('step', 1)
548
549    bdio = ctypes.cdll.LoadLibrary(bdio_path)
550
551    bdio_open = bdio.bdio_open
552    bdio_open.restype = ctypes.c_void_p
553
554    bdio_close = bdio.bdio_close
555    bdio_close.restype = ctypes.c_int
556    bdio_close.argtypes = [ctypes.c_void_p]
557
558    bdio_seek_record = bdio.bdio_seek_record
559    bdio_seek_record.restype = ctypes.c_int
560    bdio_seek_record.argtypes = [ctypes.c_void_p]
561
562    bdio_get_rlen = bdio.bdio_get_rlen
563    bdio_get_rlen.restype = ctypes.c_int
564    bdio_get_rlen.argtypes = [ctypes.c_void_p]
565
566    bdio_get_ruinfo = bdio.bdio_get_ruinfo
567    bdio_get_ruinfo.restype = ctypes.c_int
568    bdio_get_ruinfo.argtypes = [ctypes.c_void_p]
569
570    bdio_read = bdio.bdio_read
571    bdio_read.restype = ctypes.c_size_t
572    bdio_read.argtypes = [ctypes.c_char_p, ctypes.c_size_t, ctypes.c_void_p]
573
574    bdio_read_f64 = bdio.bdio_read_f64
575    bdio_read_f64.restype = ctypes.c_size_t
576    bdio_read_f64.argtypes = [ctypes.c_void_p, ctypes.c_size_t, ctypes.c_void_p]
577
578    b_path = file_path.encode('utf-8')
579    read = 'r'
580    b_read = read.encode('utf-8')
581    form = 'Generic Correlator Format 1.0'
582    b_form = form.encode('utf-8')
583
584    ensemble_name = ''
585    volume = []  # lattice volume
586    boundary_conditions = []
587    corr_name = []  # Contains correlator names
588    corr_type = []  # Contains correlator data type (important for reading out numerical data)
589    corr_props = []  # Contains propagator types (Component of corr_kappa)
590    d0 = 0  # tvals
591    # d1 = 0  # nnoise
592    prop_kappa = []  # Contains propagator kappas (Component of corr_kappa)
593    # Check noise type for multiple replica?
594    corr_no = -1
595    data = []
596    idl = []
597
598    fbdio = bdio_open(ctypes.c_char_p(b_path), ctypes.c_char_p(b_read), ctypes.c_char_p(b_form))
599
600    print('Reading of bdio file started')
601    while True:
602        bdio_seek_record(fbdio)
603        ruinfo = bdio_get_ruinfo(fbdio)
604        if ruinfo < 0:
605            # EOF reached
606            break
607        rlen = bdio_get_rlen(fbdio)
608        if ruinfo == 5:
609            d_buf = ctypes.c_double * (2 + d0)
610            pd_buf = d_buf()
611            ppd_buf = ctypes.c_void_p(ctypes.addressof(pd_buf))
612            bdio_read_f64(ppd_buf, ctypes.c_size_t(rlen), ctypes.c_void_p(fbdio))
613            tmp_mean = np.mean(np.asarray(pd_buf[2:]))
614
615            data[corr_no].append(tmp_mean)
616            corr_no += 1
617        else:
618            alt_buf = ctypes.create_string_buffer(1024)
619            palt_buf = ctypes.c_char_p(ctypes.addressof(alt_buf))
620            iread = bdio_read(palt_buf, ctypes.c_size_t(rlen), ctypes.c_void_p(fbdio))
621            if rlen != iread:
622                print('Error')
623            for i, item in enumerate(alt_buf):
624                if item == b'\x00':
625                    alt_buf[i] = b' '
626            tmp_string = (alt_buf[:].decode("utf-8")).rstrip()
627            if ruinfo == 0:
628                creator = _get_kwd(tmp_string, 'CREATOR=')
629                ensemble_name = _get_kwd(tmp_string, 'ENSEMBLE=')
630                volume.append(int(_get_kwd(tmp_string, 'L0=')))
631                volume.append(int(_get_kwd(tmp_string, 'L1=')))
632                volume.append(int(_get_kwd(tmp_string, 'L2=')))
633                volume.append(int(_get_kwd(tmp_string, 'L3=')))
634                boundary_conditions.append(_get_kwd(tmp_string, 'BC0='))
635                boundary_conditions.append(_get_kwd(tmp_string, 'BC1='))
636                boundary_conditions.append(_get_kwd(tmp_string, 'BC2='))
637                boundary_conditions.append(_get_kwd(tmp_string, 'BC3='))
638
639            if ruinfo == 1:
640                corr_name.append(_get_corr_name(tmp_string, 'CORR_NAME='))
641                corr_type.append(_get_kwd(tmp_string, 'DATATYPE='))
642                corr_props.append(_get_kwd(tmp_string, 'PROP0='))
643                if d0 == 0:
644                    d0 = int(_get_kwd(tmp_string, 'D0='))
645                else:
646                    if d0 != int(_get_kwd(tmp_string, 'D0=')):
647                        print('Error: Varying number of time values')
648            if ruinfo == 2:
649                prop_kappa.append(_get_kwd(tmp_string, 'KAPPA='))
650            if ruinfo == 4:
651                cnfg_no = int(_get_kwd(tmp_string, 'CNFG_ID='))
652                if stop:
653                    if cnfg_no > kwargs.get('stop'):
654                        break
655                idl.append(cnfg_no)
656                print('\r%s %i' % ('Reading configuration', cnfg_no), end='\r')
657                if len(idl) == 1:
658                    no_corrs = len(corr_name)
659                    data = []
660                    for c in range(no_corrs):
661                        data.append([])
662
663                corr_no = 0
664    bdio_close(fbdio)
665
666    print('\nCreator: ', creator)
667    print('Ensemble: ', ensemble_name)
668    print('Lattice volume: ', volume)
669    print('Boundary conditions: ', boundary_conditions)
670    print('Number of random sources: ', d0)
671    print('Number of corrs: ', len(corr_name))
672    print('Number of configurations: ', cnfg_no + 1)
673
674    corr_kappa = []  # Contains kappa values for both propagators of given correlation function
675    for item in corr_props:
676        corr_kappa.append(float(prop_kappa[int(item)]))
677
678    if stop is None:
679        stop = idl[-1]
680    idl_target = range(start, stop + 1, step)
681    try:
682        indices = [idl.index(i) for i in idl_target]
683    except ValueError as err:
684        raise Exception('Configurations in file do no match target list!', err)
685
686    result = {}
687    for c in range(no_corrs):
688        deltas = [np.asarray(data[c])[index] for index in indices]
689        result[(corr_name[c], str(corr_kappa[c]))] = Obs([deltas], [ensemble_name], idl=[idl_target])
690
691    # Check that all data entries have the same number of configurations
692    if len(set([o.N for o in list(result.values())])) != 1:
693        raise Exception('Error: Not all correlators have the same number of configurations. bdio file is possibly corrupted.')
694
695    return result

Extract dSdm data from a bdio file and return it as a dictionary

The dictionary can be accessed with a tuple consisting of (type, kappa)

read_dSdm requires bdio to be compiled into a shared library. This can be achieved by adding the flag -fPIC to CC and changing the all target to

all: bdio.o $(LIBDIR) gcc -shared -Wl,-soname,libbdio.so -o $(BUILDDIR)/libbdio.so $(BUILDDIR)/bdio.o cp $(BUILDDIR)/libbdio.so $(LIBDIR)/

Parameters
  • file_path (str): path to the bdio file
  • bdio_path (str): path to the shared bdio library libbdio.so (default ./libbdio.so)
  • start (int): The first configuration to be read (default 1)
  • stop (int): The last configuration to be read (default None)
  • step (int): Fixed step size between two measurements (default 1)
  • alternative_ensemble_name (str): Manually overwrite ensemble name