pyerrors.input.dobs

  1import datetime
  2import getpass
  3import gzip
  4import json
  5import socket
  6import warnings
  7from collections import defaultdict
  8
  9import lxml.etree as et
 10import numpy as np
 11
 12from .. import version as pyerrorsversion
 13from ..covobs import Covobs
 14from ..obs import Obs, _merge_idx
 15
 16
 17# Based on https://stackoverflow.com/a/10076823
 18def _etree_to_dict(t):
 19    """ Convert the content of an XML file to a python dict"""
 20    d = {t.tag: {} if t.attrib else None}
 21    children = list(t)
 22    if children:
 23        dd = defaultdict(list)
 24        for dc in map(_etree_to_dict, children):
 25            for k, v in dc.items():
 26                dd[k].append(v)
 27        d = {t.tag: {k: v[0] if len(v) == 1 else v
 28                     for k, v in dd.items()}}
 29    if t.attrib:
 30        d[t.tag].update(('@' + k, v)
 31                        for k, v in t.attrib.items())
 32    if t.text:
 33        text = t.text.strip()
 34        if children or t.attrib:
 35            if text:
 36                d[t.tag]['#data'] = [text]
 37        else:
 38            d[t.tag] = text
 39    return d
 40
 41
 42def _dict_to_xmlstring(d):
 43    if isinstance(d, dict):
 44        iters = ''
 45        for k in d:
 46            if k.startswith('#'):
 47                for la in d[k]:
 48                    iters += la
 49                iters = '<array>\n' + iters + '<{}array>\n'.format('/')
 50                return iters
 51            if isinstance(d[k], dict):
 52                iters += f'<{k}>\n' + _dict_to_xmlstring(d[k]) + '<{}{}>\n'.format('/', k)
 53            elif isinstance(d[k], str):
 54                if len(d[k]) > 100:
 55                    iters += f'<{k}>\n ' + d[k] + ' \n<{}{}>\n'.format('/', k)
 56                else:
 57                    iters += f'<{k}> ' + d[k] + ' <{}{}>\n'.format('/', k)
 58            elif isinstance(d[k], list):
 59                for i in range(len(d[k])):
 60                    iters += _dict_to_xmlstring(d[k][i])
 61            elif not d[k]:
 62                return '\n'
 63            else:
 64                raise Exception('Type', type(d[k]), 'not supported in export!')
 65    else:
 66        raise Exception('Type', type(d), 'not supported in export!')
 67    return iters
 68
 69
 70def _dict_to_xmlstring_spaces(d, space='  '):
 71    s = _dict_to_xmlstring(d)
 72    o = ''
 73    c = 0
 74    cm = False
 75    for li in s.split('\n'):
 76        if li.startswith('<{}'.format('/')):
 77            c -= 1
 78            cm = True
 79        for _i in range(c):
 80            o += space
 81        o += li + '\n'
 82        if li.startswith('<') and not cm:
 83            if '<{}'.format('/') not in li:
 84                c += 1
 85        cm = False
 86    return o
 87
 88
 89def create_pobs_string(obsl, name, spec='', origin='', symbol=None, enstag=None):
 90    """Export a list of Obs or structures containing Obs to an xml string
 91    according to the Zeuthen pobs format.
 92
 93    Tags are not written or recovered automatically. The separator | is removed from the replica names.
 94
 95    Parameters
 96    ----------
 97    obsl : list
 98        List of Obs that will be exported.
 99        The Obs inside a structure have to be defined on the same ensemble.
100    name : str
101        The name of the observable.
102    spec : str
103        Optional string that describes the contents of the file.
104    origin : str
105        Specify where the data has its origin.
106    symbol : list
107        A list of symbols that describe the observables to be written. May be empty.
108    enstag : str
109        Enstag that is written to pobs. If None, the ensemble name is used.
110
111    Returns
112    -------
113    xml_str : str
114        XML formatted string of the input data
115    """
116
117    if symbol is None:
118        symbol = []
119
120    od = {}
121    ename = obsl[0].e_names[0]
122    names = list(obsl[0].deltas.keys())
123    nr = len(names)
124    onames = [name.replace('|', '') for name in names]
125    for o in obsl:
126        if len(o.e_names) != 1:
127            raise Exception('You try to export dobs to obs!')
128        if o.e_names[0] != ename:
129            raise Exception('You try to export dobs to obs!')
130        if len(o.deltas.keys()) != nr:
131            raise Exception('Incompatible obses in list')
132    od['observables'] = {}
133    od['observables']['schema'] = {'name': 'lattobs', 'version': '1.0'}
134    od['observables']['origin'] = {
135        'who': getpass.getuser(),
136        'date': str(datetime.datetime.now())[:-7],
137        'host': socket.gethostname(),
138        'tool': {'name': 'pyerrors', 'version': pyerrorsversion.__version__}}
139    od['observables']['pobs'] = {}
140    pd = od['observables']['pobs']
141    pd['spec'] = spec
142    pd['origin'] = origin
143    pd['name'] = name
144    if enstag:
145        if not isinstance(enstag, str):
146            raise Exception('enstag has to be a string!')
147        pd['enstag'] = enstag
148    else:
149        pd['enstag'] = ename
150    pd['nr'] = f'{nr}'
151    pd['array'] = []
152    osymbol = 'cfg'
153    if not isinstance(symbol, list):
154        raise Exception('Symbol has to be a list!')
155    if not (len(symbol) == 0 or len(symbol) == len(obsl)):
156        raise Exception(f'Symbol has to be a list of lenght 0 or {len(obsl)}!')
157    for s in symbol:
158        osymbol += f' {s}'
159    for r in range(nr):
160        ad = {}
161        ad['id'] = onames[r]
162        Nconf = len(obsl[0].deltas[names[r]])
163        layout = f'{Nconf} i f{len(obsl)}'
164        ad['layout'] = layout
165        ad['symbol'] = osymbol
166        data = ''
167        for c in range(Nconf):
168            data += f'{obsl[0].idl[names[r]][c]} '
169            for o in obsl:
170                num = o.deltas[names[r]][c] + o.r_values[names[r]]
171                if num == 0:
172                    data += '0 '
173                else:
174                    data += f'{num:1.16e} '
175            data += '\n'
176        ad['#data'] = data
177        pd['array'].append(ad)
178
179    rs = '<?xml version="1.0" encoding="utf-8"?>\n' + _dict_to_xmlstring_spaces(od)
180    return rs
181
182
183def write_pobs(obsl, fname, name, spec='', origin='', symbol=None, enstag=None, gz=True):
184    """Export a list of Obs or structures containing Obs to a .xml.gz file
185    according to the Zeuthen pobs format.
186
187    Tags are not written or recovered automatically. The separator | is removed from the replica names.
188
189    Parameters
190    ----------
191    obsl : list
192        List of Obs that will be exported.
193        The Obs inside a structure have to be defined on the same ensemble.
194    fname : str
195        Filename of the output file.
196    name : str
197        The name of the observable.
198    spec : str
199        Optional string that describes the contents of the file.
200    origin : str
201        Specify where the data has its origin.
202    symbol : list
203        A list of symbols that describe the observables to be written. May be empty.
204    enstag : str
205        Enstag that is written to pobs. If None, the ensemble name is used.
206    gz : bool
207        If True, the output is a gzipped xml. If False, the output is an xml file.
208
209    Returns
210    -------
211    None
212    """
213    pobsstring = create_pobs_string(obsl, name, spec, origin, symbol, enstag)
214
215    if not fname.endswith('.xml') and not fname.endswith('.gz'):
216        fname += '.xml'
217
218    if gz:
219        if not fname.endswith('.gz'):
220            fname += '.gz'
221
222        fp = gzip.open(fname, 'wb')
223        fp.write(pobsstring.encode('utf-8'))
224    else:
225        fp = open(fname, 'w', encoding='utf-8')
226        fp.write(pobsstring)
227    fp.close()
228
229
230def _import_data(string):
231    return json.loads("[" + ",".join(string.replace(' +', ' ').split()) + "]")
232
233
234def _check(condition):
235    if not condition:
236        raise Exception("XML file format not supported")
237
238
239class _NoTagInDataError(Exception):
240    """Raised when tag is not in data"""
241    def __init__(self, tag):
242        self.tag = tag
243        super().__init__(f'Tag {self.tag} not in data!')
244
245
246def _find_tag(dat, tag):
247    for i in range(len(dat)):
248        if dat[i].tag == tag:
249            return i
250    raise _NoTagInDataError(tag)
251
252
253def _import_array(arr):
254    name = arr[_find_tag(arr, 'id')].text.strip()
255    index = _find_tag(arr, 'layout')
256    try:
257        sindex = _find_tag(arr, 'symbol')
258    except _NoTagInDataError:
259        sindex = 0
260    if sindex > index:
261        tmp = _import_data(arr[sindex].tail)
262    else:
263        tmp = _import_data(arr[index].tail)
264
265    li = arr[index].text.strip()
266    m = li.split()
267    if m[1] == "i" and m[2][0] == "f":
268        nc = int(m[0])
269        na = int(m[2].lstrip('f'))
270        _dat = []
271        mask = []
272        for a in range(na):
273            mask += [a]
274            _dat += [np.array(tmp[1 + a:: na + 1])]
275        _check(len(tmp[0:: na + 1]) == nc)
276        return [name, tmp[0:: na + 1], mask, _dat]
277    elif m[1][0] == 'f' and len(m) < 3:
278        sh = (int(m[0]), int(m[1].lstrip('f')))
279        return np.reshape(tmp, sh)
280    elif any(['f' in s for s in m]):
281        for si in range(len(m)):
282            if m[si] == 'f':
283                break
284        sh = [int(m[i]) for i in range(si)]
285        return np.reshape(tmp, sh)
286    else:
287        print(name, m)
288        _check(False)
289
290
291def _import_rdata(rd):
292    name, idx, _mask, deltas = _import_array(rd)
293    return deltas, name, idx
294
295
296def _import_cdata(cd):
297    _check(cd[0].tag == "id")
298    _check(cd[1][0].text.strip() == "cov")
299    cov = _import_array(cd[1])
300    grad = _import_array(cd[2])
301    return cd[0].text.strip(), cov, grad
302
303
304def read_pobs(fname, full_output=False, gz=True, separator_insertion=None):
305    """Import a list of Obs from an xml.gz file in the Zeuthen pobs format.
306
307    Tags are not written or recovered automatically.
308
309    Parameters
310    ----------
311    fname : str
312        Filename of the input file.
313    full_output : bool
314        If True, a dict containing auxiliary information and the data is returned.
315        If False, only the data is returned as list.
316    separatior_insertion: str or int
317        str: replace all occurences of "separator_insertion" within the replica names
318        by "|%s" % (separator_insertion) when constructing the names of the replica.
319        int: Insert the separator "|" at the position given by separator_insertion.
320        None (default): Replica names remain unchanged.
321
322    Returns
323    -------
324    res : list[Obs]
325        Imported data
326    or
327    res : dict
328        Imported data and meta-data
329    """
330
331    if not fname.endswith('.xml') and not fname.endswith('.gz'):
332        fname += '.xml'
333    if gz:
334        if not fname.endswith('.gz'):
335            fname += '.gz'
336        with gzip.open(fname, 'r') as fin:
337            content = fin.read()
338    else:
339        if fname.endswith('.gz'):
340            warnings.warn(f"Trying to read from {fname} without unzipping!", UserWarning, stacklevel=2)
341        with open(fname) as fin:
342            content = fin.read()
343
344    # parse xml file content
345    root = et.fromstring(content)
346
347    _check(root[2].tag == 'pobs')
348    pobs = root[2]
349
350    version = root[0][1].text.strip()
351
352    _check(root[1].tag == 'origin')
353    file_origin = _etree_to_dict(root[1])['origin']
354
355    deltas = []
356    names = []
357    idl = []
358    for i in range(5, len(pobs)):
359        delta, name, idx = _import_rdata(pobs[i])
360        deltas.append(delta)
361        if separator_insertion is None:
362            pass
363        elif isinstance(separator_insertion, int):
364            name = name[:separator_insertion] + '|' + name[separator_insertion:]
365        elif isinstance(separator_insertion, str):
366            name = name.replace(separator_insertion, f"|{separator_insertion}")
367        else:
368            raise Exception("separator_insertion has to be string or int, is ", type(separator_insertion))
369        names.append(name)
370        idl.append(idx)
371    res = [Obs([d[i] for d in deltas], names, idl=idl) for i in range(len(deltas[0]))]
372
373    descriptiond = {}
374    for i in range(4):
375        descriptiond[pobs[i].tag] = pobs[i].text.strip()
376
377    _check(pobs[4].tag == "nr")
378
379    _check(pobs[5].tag == 'array')
380    if pobs[5][1].tag == 'symbol':
381        symbol = pobs[5][1].text.strip()
382        descriptiond['symbol'] = symbol
383
384    if full_output:
385        retd = {}
386        tool = file_origin.get('tool', None)
387        if tool:
388            program = tool['name'] + ' ' + tool['version']
389        else:
390            program = ''
391        retd['program'] = program
392        retd['version'] = version
393        retd['who'] = file_origin['who']
394        retd['date'] = file_origin['date']
395        retd['host'] = file_origin['host']
396        retd['description'] = descriptiond
397        retd['obsdata'] = res
398        return retd
399    else:
400        return res
401
402
403# this is based on Mattia Bruno's implementation at https://github.com/mbruno46/pyobs/blob/master/pyobs/IO/xml.py
404def import_dobs_string(content, full_output=False, separator_insertion=True):
405    """Import a list of Obs from a string in the Zeuthen dobs format.
406
407    Tags are not written or recovered automatically.
408
409    Parameters
410    ----------
411    content : str
412        XML string containing the data
413    full_output : bool
414        If True, a dict containing auxiliary information and the data is returned.
415        If False, only the data is returned as list.
416    separatior_insertion: str, int or bool
417        str: replace all occurences of "separator_insertion" within the replica names
418        by "|%s" % (separator_insertion) when constructing the names of the replica.
419        int: Insert the separator "|" at the position given by separator_insertion.
420        True (default): separator "|" is inserted after len(ensname), assuming that the
421        ensemble name is a prefix to the replica name.
422        None or False: No separator is inserted.
423
424    Returns
425    -------
426    res : list[Obs]
427        Imported data
428    or
429    res : dict
430        Imported data and meta-data
431    """
432
433    root = et.fromstring(content)
434
435    _check(root.tag == 'OBSERVABLES')
436    _check(root[0].tag == 'SCHEMA')
437    version = root[0][1].text.strip()
438
439    _check(root[1].tag == 'origin')
440    file_origin = _etree_to_dict(root[1])['origin']
441
442    _check(root[2].tag == 'dobs')
443
444    dobs = root[2]
445
446    descriptiond = {}
447    for i in range(3):
448        descriptiond[dobs[i].tag] = dobs[i].text.strip()
449
450    _check(dobs[3].tag == 'array')
451
452    symbol = []
453    if dobs[3][1].tag == 'symbol':
454        symbol = dobs[3][1].text.strip()
455        descriptiond['symbol'] = symbol
456    mean = _import_array(dobs[3])[0]
457
458    _check(dobs[4].tag == "ne")
459    ne = int(dobs[4].text.strip())
460    _check(dobs[5].tag == "nc")
461
462    idld = {}
463    deltad = {}
464    covd = {}
465    gradd = {}
466    names = []
467    e_names = []
468    enstags = {}
469    for k in range(6, len(list(dobs))):
470        if dobs[k].tag == "edata":
471            _check(dobs[k][0].tag == "enstag")
472            ename = dobs[k][0].text.strip()
473            e_names.append(ename)
474            _check(dobs[k][1].tag == "nr")
475            R = int(dobs[k][1].text.strip())
476            for i in range(2, 2 + R):
477                deltas, rname, idx = _import_rdata(dobs[k][i])
478                if separator_insertion is None or False:
479                    pass
480                elif separator_insertion is True:
481                    if rname.startswith(ename):
482                        rname = rname[:len(ename)] + '|' + rname[len(ename):]
483                elif isinstance(separator_insertion, int):
484                    rname = rname[:separator_insertion] + '|' + rname[separator_insertion:]
485                elif isinstance(separator_insertion, str):
486                    rname = rname.replace(separator_insertion, f"|{separator_insertion}")
487                else:
488                    raise Exception("separator_insertion has to be string or int, is ", type(separator_insertion))
489                if '|' in rname:
490                    new_ename = rname[:rname.index('|')]
491                else:
492                    new_ename = ename
493                enstags[new_ename] = ename
494                idld[rname] = idx
495                deltad[rname] = deltas
496                names.append(rname)
497        elif dobs[k].tag == "cdata":
498            cname, cov, grad = _import_cdata(dobs[k])
499            covd[cname] = cov
500            if grad.shape[1] == 1:
501                gradd[cname] = [grad for i in range(len(mean))]
502            else:
503                gradd[cname] = grad.T
504        else:
505            _check(False)
506    names = list(set(names))
507
508    for name in names:
509        for i in range(len(deltad[name])):
510            tmp = np.zeros_like(deltad[name][i])
511            for j in range(len(deltad[name][i])):
512                if deltad[name][i][j] != 0.:
513                    tmp[j] = deltad[name][i][j] + mean[i]
514            deltad[name][i] = tmp
515
516    res = []
517    for i in range(len(mean)):
518        deltas = []
519        idl = []
520        obs_names = []
521        for name in names:
522            h = np.unique(deltad[name][i])
523            if len(h) == 1 and np.all(h == mean[i]):
524                continue
525            repdeltas = []
526            repidl = []
527            for j in range(len(deltad[name][i])):
528                if deltad[name][i][j] != 0.:
529                    repdeltas.append(deltad[name][i][j])
530                    repidl.append(idld[name][j])
531            if len(repdeltas) > 0:
532                obs_names.append(name)
533                deltas.append(repdeltas)
534                idl.append(repidl)
535
536        obsmeans = [np.average(deltas[j]) for j in range(len(deltas))]
537        res.append(Obs([np.array(deltas[j]) - obsmeans[j] for j in range(len(obsmeans))], obs_names, idl=idl, means=obsmeans))
538        res[-1]._value = mean[i]
539    _check(len(e_names) == ne)
540
541    cnames = list(covd.keys())
542    for i in range(len(res)):
543        new_covobs = {name: Covobs(0, covd[name], name, grad=gradd[name][i]) for name in cnames}
544        for name in cnames:
545            if np.all(new_covobs[name].grad == 0):
546                del new_covobs[name]
547        cnames_loc = list(new_covobs.keys())
548        for name in cnames_loc:
549            res[i].names.append(name)
550            res[i].shape[name] = 1
551            res[i].idl[name] = []
552        res[i]._covobs = new_covobs
553
554    if symbol:
555        for i in range(len(res)):
556            res[i].tag = symbol[i]
557            if res[i].tag == 'None':
558                res[i].tag = None
559    if full_output:
560        retd = {}
561        tool = file_origin.get('tool', None)
562        if tool:
563            program = tool['name'] + ' ' + tool['version']
564        else:
565            program = ''
566        retd['program'] = program
567        retd['version'] = version
568        retd['who'] = file_origin['who']
569        retd['date'] = file_origin['date']
570        retd['host'] = file_origin['host']
571        retd['description'] = descriptiond
572        retd['enstags'] = enstags
573        retd['obsdata'] = res
574        return retd
575    else:
576        return res
577
578
579def read_dobs(fname, full_output=False, gz=True, separator_insertion=True):
580    """Import a list of Obs from an xml.gz file in the Zeuthen dobs format.
581
582    Tags are not written or recovered automatically.
583
584    Parameters
585    ----------
586    fname : str
587        Filename of the input file.
588    full_output : bool
589        If True, a dict containing auxiliary information and the data is returned.
590        If False, only the data is returned as list.
591    gz : bool
592        If True, assumes that data is gzipped. If False, assumes XML file.
593    separatior_insertion: str, int or bool
594        str: replace all occurences of "separator_insertion" within the replica names
595        by "|%s" % (separator_insertion) when constructing the names of the replica.
596        int: Insert the separator "|" at the position given by separator_insertion.
597        True (default): separator "|" is inserted after len(ensname), assuming that the
598        ensemble name is a prefix to the replica name.
599        None or False: No separator is inserted.
600
601    Returns
602    -------
603    res : list[Obs]
604        Imported data
605    or
606    res : dict
607        Imported data and meta-data
608    """
609
610    if not fname.endswith('.xml') and not fname.endswith('.gz'):
611        fname += '.xml'
612    if gz:
613        if not fname.endswith('.gz'):
614            fname += '.gz'
615        with gzip.open(fname, 'r') as fin:
616            content = fin.read()
617    else:
618        if fname.endswith('.gz'):
619            warnings.warn(f"Trying to read from {fname} without unzipping!", UserWarning, stacklevel=2)
620        with open(fname) as fin:
621            content = fin.read()
622
623    return import_dobs_string(content, full_output, separator_insertion=separator_insertion)
624
625
626def _dobsdict_to_xmlstring(d):
627    if isinstance(d, dict):
628        iters = ''
629        for k in d:
630            if k.startswith('#value'):
631                for li in d[k]:
632                    iters += li
633                return iters + '\n'
634            elif k.startswith('#'):
635                for li in d[k]:
636                    iters += li
637                iters = '<array>\n' + iters + '<{}array>\n'.format('/')
638                return iters
639            if isinstance(d[k], dict):
640                iters += f'<{k}>\n' + _dobsdict_to_xmlstring(d[k]) + '<{}{}>\n'.format('/', k)
641            elif isinstance(d[k], str):
642                if len(d[k]) > 100:
643                    iters += f'<{k}>\n ' + d[k] + ' \n<{}{}>\n'.format('/', k)
644                else:
645                    iters += f'<{k}> ' + d[k] + ' <{}{}>\n'.format('/', k)
646            elif isinstance(d[k], list):
647                tmps = ''
648                if k in ['edata', 'cdata']:
649                    for i in range(len(d[k])):
650                        tmps += f'<{k}>\n' + _dobsdict_to_xmlstring(d[k][i]) + f'</{k}>\n'
651                else:
652                    for i in range(len(d[k])):
653                        tmps += _dobsdict_to_xmlstring(d[k][i])
654                iters += tmps
655            elif isinstance(d[k], (int, float)):
656                iters += f'<{k}> ' + str(d[k]) + ' <{}{}>\n'.format('/', k)
657            elif not d[k]:
658                return '\n'
659            else:
660                raise Exception('Type', type(d[k]), 'not supported in export!')
661    else:
662        raise Exception('Type', type(d), 'not supported in export!')
663    return iters
664
665
666def _dobsdict_to_xmlstring_spaces(d, space='  '):
667    s = _dobsdict_to_xmlstring(d)
668    o = ''
669    c = 0
670    cm = False
671    for li in s.split('\n'):
672        if li.startswith('<{}'.format('/')):
673            c -= 1
674            cm = True
675        for _i in range(c):
676            o += space
677        o += li + '\n'
678        if li.startswith('<') and not cm:
679            if '<{}'.format('/') not in li:
680                c += 1
681        cm = False
682    return o
683
684
685def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=None, who=None, enstags=None):
686    """Generate the string for the export of a list of Obs or structures containing Obs
687    to a .xml.gz file according to the Zeuthen dobs format.
688
689    Tags are not written or recovered automatically. The separator |is removed from the replica names.
690
691    Parameters
692    ----------
693    obsl : list
694        List of Obs that will be exported.
695        The Obs inside a structure do not have to be defined on the same set of configurations,
696        but the storage requirement is increased, if this is not the case.
697    name : str
698        The name of the observable.
699    spec : str
700        Optional string that describes the contents of the file.
701    origin : str
702        Specify where the data has its origin.
703    symbol : list
704        A list of symbols that describe the observables to be written. May be empty.
705    who : str
706        Provide the name of the person that exports the data.
707    enstags : dict
708        Provide alternative enstag for ensembles in the form enstags = {ename: enstag}
709        Otherwise, the ensemble name is used.
710
711    Returns
712    -------
713    xml_str : str
714        XML string generated from the data
715    """
716    if enstags is None:
717        enstags = {}
718    if symbol is None:
719        symbol = []
720    od = {}
721    r_names = []
722    for o in obsl:
723        r_names += [name for name in o.names if name.split('|')[0] in o.mc_names]
724    r_names = sorted(set(r_names))
725    mc_names = sorted(set([n.split('|')[0] for n in r_names]))
726    for tmpname in mc_names:
727        if tmpname not in enstags:
728            enstags[tmpname] = tmpname
729    ne = len(set(mc_names))
730    cov_names = []
731    for o in obsl:
732        cov_names += list(o.cov_names)
733    cov_names = sorted(set(cov_names))
734    nc = len(set(cov_names))
735    od['OBSERVABLES'] = {}
736    od['OBSERVABLES']['SCHEMA'] = {'NAME': 'lattobs', 'VERSION': '1.0'}
737    if who is None:
738        who = getpass.getuser()
739    od['OBSERVABLES']['origin'] = {
740        'who': who,
741        'date': str(datetime.datetime.now())[:-7],
742        'host': socket.gethostname(),
743        'tool': {'name': 'pyerrors', 'version': pyerrorsversion.__version__}}
744    od['OBSERVABLES']['dobs'] = {}
745    pd = od['OBSERVABLES']['dobs']
746    pd['spec'] = spec
747    pd['origin'] = origin
748    pd['name'] = name
749    pd['array'] = {}
750    pd['array']['id'] = 'val'
751    pd['array']['layout'] = f'1 f{len(obsl)}'
752    osymbol = ''
753    if symbol:
754        if not isinstance(symbol, list):
755            raise Exception('Symbol has to be a list!')
756        if not (len(symbol) == 0 or len(symbol) == len(obsl)):
757            raise Exception(f'Symbol has to be a list of lenght 0 or {len(obsl)}!')
758        osymbol = symbol[0]
759        for s in symbol[1:]:
760            osymbol += f' {s}'
761        pd['array']['symbol'] = osymbol
762
763    pd['array']['#values'] = ['  '.join([f'{o.value:1.16e}' for o in obsl])]
764    pd['ne'] = f'{ne}'
765    pd['nc'] = f'{nc}'
766    pd['edata'] = []
767    for name in mc_names:
768        ed = {}
769        ed['enstag'] = enstags[name]
770        onames = sorted([n for n in r_names if (n.startswith(name + '|') or n == name)])
771        nr = len(onames)
772        ed['nr'] = nr
773        ed[''] = []
774
775        for r in range(nr):
776            ad = {}
777            repname = onames[r]
778            ad['id'] = repname.replace('|', '')
779            idx = _merge_idx([o.idl.get(repname, []) for o in obsl])
780            Nconf = len(idx)
781            layout = f'{Nconf} i f{len(obsl)}'
782            ad['layout'] = layout
783            data = ''
784            counters = [0 for o in obsl]
785            offsets = [o.r_values[repname] - o.value if repname in o.r_values else 0 for o in obsl]
786            for ci in idx:
787                data += f'{ci} '
788                for oi in range(len(obsl)):
789                    o = obsl[oi]
790                    if repname in o.idl:
791                        if counters[oi] < 0:
792                            num = 0
793                            if num == 0:
794                                data += '0 '
795                            else:
796                                data += f'{num:1.16e} '
797                            continue
798                        if o.idl[repname][counters[oi]] == ci:
799                            num = o.deltas[repname][counters[oi]] + offsets[oi]
800                            if num == 0:
801                                data += '0 '
802                            else:
803                                data += f'{num:1.16e} '
804                            counters[oi] += 1
805                            if counters[oi] >= len(o.idl[repname]):
806                                counters[oi] = -1
807                        else:
808                            num = 0
809                            if num == 0:
810                                data += '0 '
811                            else:
812                                data += f'{num:1.16e} '
813                    else:
814                        data += '0 '
815                data += '\n'
816            ad['#data'] = data
817            ed[''].append(ad)
818        pd['edata'].append(ed)
819
820        allcov = {}
821        for o in obsl:
822            for cname in o.cov_names:
823                if cname in allcov:
824                    if not np.array_equal(allcov[cname], o.covobs[cname].cov):
825                        raise Exception(f'Inconsistent covariance matrices for {cname}!')
826                else:
827                    allcov[cname] = o.covobs[cname].cov
828        pd['cdata'] = []
829        for cname in cov_names:
830            cd = {}
831            cd['id'] = cname
832
833            covd = {'id': 'cov'}
834            if allcov[cname].shape == ():
835                ncov = 1
836                covd['layout'] = '1 1 f'
837                covd['#data'] = f'{allcov[cname]:1.14e}'
838            else:
839                shape = allcov[cname].shape
840                assert (shape[0] == shape[1])
841                ncov = shape[0]
842                covd['layout'] = f'{ncov} {ncov} f'
843                ds = ''
844                for i in range(ncov):
845                    for j in range(ncov):
846                        val = allcov[cname][i][j]
847                        if val == 0:
848                            ds += '0 '
849                        else:
850                            ds += f'{val:1.14e} '
851                    ds += '\n'
852                covd['#data'] = ds
853
854            gradd = {'id': 'grad'}
855            gradd['layout'] = f'{ncov} f{len(obsl)}'
856            ds = ''
857            for i in range(ncov):
858                for o in obsl:
859                    if cname in o.covobs:
860                        val = o.covobs[cname].grad[i].item()
861                        if val != 0:
862                            ds += f'{val:1.14e} '
863                        else:
864                            ds += '0 '
865                    else:
866                        ds += '0 '
867            gradd['#data'] = ds
868            cd['array'] = [covd, gradd]
869            pd['cdata'].append(cd)
870
871    rs = '<?xml version="1.0" encoding="utf-8"?>\n' + _dobsdict_to_xmlstring_spaces(od)
872
873    return rs
874
875
876def write_dobs(obsl, fname, name, spec='dobs v1.0', origin='', symbol=None, who=None, enstags=None, gz=True):
877    """Export a list of Obs or structures containing Obs to a .xml.gz file
878    according to the Zeuthen dobs format.
879
880    Tags are not written or recovered automatically. The separator | is removed from the replica names.
881
882    Parameters
883    ----------
884    obsl : list
885        List of Obs that will be exported.
886        The Obs inside a structure do not have to be defined on the same set of configurations,
887        but the storage requirement is increased, if this is not the case.
888    fname : str
889        Filename of the output file.
890    name : str
891        The name of the observable.
892    spec : str
893        Optional string that describes the contents of the file.
894    origin : str
895        Specify where the data has its origin.
896    symbol : list
897        A list of symbols that describe the observables to be written. May be empty.
898    who : str
899        Provide the name of the person that exports the data.
900    enstags : dict
901        Provide alternative enstag for ensembles in the form enstags = {ename: enstag}
902        Otherwise, the ensemble name is used.
903    gz : bool
904        If True, the output is a gzipped XML. If False, the output is a XML file.
905
906    Returns
907    -------
908    None
909    """
910    if enstags is None:
911        enstags = {}
912
913    dobsstring = create_dobs_string(obsl, name, spec, origin, symbol, who, enstags=enstags)
914
915    if not fname.endswith('.xml') and not fname.endswith('.gz'):
916        fname += '.xml'
917
918    if gz:
919        if not fname.endswith('.gz'):
920            fname += '.gz'
921
922        fp = gzip.open(fname, 'wb')
923        fp.write(dobsstring.encode('utf-8'))
924    else:
925        fp = open(fname, 'w', encoding='utf-8')
926        fp.write(dobsstring)
927    fp.close()