Merge pull request #112 from s-kuberski/feature/dobs

Feature/dobs
This commit is contained in:
Fabian Joswig 2022-06-19 12:50:40 +01:00 committed by GitHub
commit 3f1fbd1869
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 36 additions and 34 deletions

View file

@ -379,10 +379,6 @@ def read_pobs(fname, full_output=False, gz=True, separator_insertion=None):
return res return res
# Reading (and writing) dobs is not yet working properly:
# we have to loop over root[2:] because each entry is a dobs
# But maybe this is just a problem with Ben's implementation
# this is based on Mattia Bruno's implementation at https://github.com/mbruno46/pyobs/blob/master/pyobs/IO/xml.py # this is based on Mattia Bruno's implementation at https://github.com/mbruno46/pyobs/blob/master/pyobs/IO/xml.py
def import_dobs_string(content, noempty=False, full_output=False, separator_insertion=True): def import_dobs_string(content, noempty=False, full_output=False, separator_insertion=True):
"""Import a list of Obs from a string in the Zeuthen dobs format. """Import a list of Obs from a string in the Zeuthen dobs format.
@ -501,7 +497,7 @@ def import_dobs_string(content, noempty=False, full_output=False, separator_inse
obs_names.append(name) obs_names.append(name)
idl.append(idld[name]) idl.append(idld[name])
res.append(Obs(deltas, obs_names, idl=idl)) res.append(Obs(deltas, obs_names, idl=idl))
print(mean, 'vs', res) res[-1]._value = mean[i]
_check(len(e_names) == ne) _check(len(e_names) == ne)
cnames = list(covd.keys()) cnames = list(covd.keys())
@ -579,17 +575,13 @@ def read_dobs(fname, noempty=False, full_output=False, gz=True, separator_insert
if not fname.endswith('.gz'): if not fname.endswith('.gz'):
fname += '.gz' fname += '.gz'
with gzip.open(fname, 'r') as fin: with gzip.open(fname, 'r') as fin:
content = fin.read().decode('utf-8') content = fin.read()
else: else:
if fname.endswith('.gz'): if fname.endswith('.gz'):
warnings.warn("Trying to read from %s without unzipping!" % fname, UserWarning) warnings.warn("Trying to read from %s without unzipping!" % fname, UserWarning)
with open(fname, 'r', encoding='utf-8') as fin: with open(fname, 'r') as fin:
content = fin.read() content = fin.read()
# open and read gzipped xml file
infile = gzip.open(fname)
content = infile.read()
return import_dobs_string(content, noempty, full_output, separator_insertion=separator_insertion) return import_dobs_string(content, noempty, full_output, separator_insertion=separator_insertion)
@ -745,16 +737,21 @@ def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=N
ad['layout'] = layout ad['layout'] = layout
data = '' data = ''
counters = [0 for o in obsl] counters = [0 for o in obsl]
offsets = [o.r_values[repname] - o.value if repname in o.r_values else 0 for o in obsl]
for ci in idx: for ci in idx:
data += '%d ' % ci data += '%d ' % ci
for oi in range(len(obsl)): for oi in range(len(obsl)):
o = obsl[oi] o = obsl[oi]
if repname in o.idl: if repname in o.idl:
if counters[oi] < 0: if counters[oi] < 0:
num = offsets[oi]
if num == 0:
data += '0 ' data += '0 '
else:
data += '%1.16e ' % (num)
continue continue
if o.idl[repname][counters[oi]] == ci: if o.idl[repname][counters[oi]] == ci:
num = o.deltas[repname][counters[oi]] num = o.deltas[repname][counters[oi]] + offsets[oi]
if num == 0: if num == 0:
data += '0 ' data += '0 '
else: else:
@ -763,7 +760,11 @@ def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=N
if counters[oi] >= len(o.idl[repname]): if counters[oi] >= len(o.idl[repname]):
counters[oi] = -1 counters[oi] = -1
else: else:
num = offsets[oi]
if num == 0:
data += '0 ' data += '0 '
else:
data += '%1.16e ' % (num)
else: else:
data += '0 ' data += '0 '
data += '\n' data += '\n'
@ -773,31 +774,31 @@ def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=N
allcov = {} allcov = {}
for o in obsl: for o in obsl:
for name in o.cov_names: for cname in o.cov_names:
if name in allcov: if cname in allcov:
if not np.array_equal(allcov[name], o.covobs[name].cov): if not np.array_equal(allcov[cname], o.covobs[cname].cov):
raise Exception('Inconsistent covariance matrices for %s!' % (name)) raise Exception('Inconsistent covariance matrices for %s!' % (cname))
else: else:
allcov[name] = o.covobs[name].cov allcov[cname] = o.covobs[cname].cov
pd['cdata'] = [] pd['cdata'] = []
for name in cov_names: for cname in cov_names:
cd = {} cd = {}
cd['id'] = name cd['id'] = cname
covd = {'id': 'cov'} covd = {'id': 'cov'}
if allcov[name].shape == (): if allcov[cname].shape == ():
ncov = 1 ncov = 1
covd['layout'] = '1 1 f' covd['layout'] = '1 1 f'
covd['#data'] = '%1.14e' % (allcov[name]) covd['#data'] = '%1.14e' % (allcov[cname])
else: else:
shape = allcov[name].shape shape = allcov[cname].shape
assert (shape[0] == shape[1]) assert (shape[0] == shape[1])
ncov = shape[0] ncov = shape[0]
covd['layout'] = '%d %d f' % (ncov, ncov) covd['layout'] = '%d %d f' % (ncov, ncov)
ds = '' ds = ''
for i in range(ncov): for i in range(ncov):
for j in range(ncov): for j in range(ncov):
val = allcov[name][i][j] val = allcov[cname][i][j]
if val == 0: if val == 0:
ds += '0 ' ds += '0 '
else: else:
@ -810,8 +811,8 @@ def create_dobs_string(obsl, name, spec='dobs v1.0', origin='', symbol=[], who=N
ds = '' ds = ''
for i in range(ncov): for i in range(ncov):
for o in obsl: for o in obsl:
if name in o.covobs: if cname in o.covobs:
val = o.covobs[name].grad[i] val = o.covobs[cname].grad[i]
if val != 0: if val != 0:
ds += '%1.14e ' % (val) ds += '%1.14e ' % (val)
else: else:

View file

@ -323,15 +323,16 @@ def test_dobsio():
o5 /= co2[0] o5 /= co2[0]
o5.tag = 2 * otag o5.tag = 2 * otag
tt1 = pe.Obs([np.random.rand(100)], ['t|r1'], idl=[range(2, 202, 2)]) tt1 = pe.Obs([np.random.rand(100), np.random.rand(100)], ['t|r1', 't|r2'], idl=[range(2, 202, 2), range(22, 222, 2)])
tt2 = pe.Obs([np.random.rand(100)], ['t|r2'], idl=[range(2, 202, 2)])
tt3 = pe.Obs([np.random.rand(102)], ['qe|r1']) tt3 = pe.Obs([np.random.rand(102)], ['qe|r1'])
tt = tt1 + tt2 + tt3 tt = tt1 + tt3
tt.tag = 'Test Obs: Ä' tt.tag = 'Test Obs: Ä'
ol = [o2, o3, o4, do, o5, tt] tt4 = pe.Obs([np.random.rand(100), np.random.rand(100)], ['t|r1', 't|r2'], idl=[range(1, 101, 1), range(2, 202, 2)])
ol = [o2, o3, o4, do, o5, tt, tt4, np.log(tt4 / o5**2), np.exp(o5 + np.log(co3 / tt3 + o4) / tt)]
print(ol) print(ol)
fname = 'test_rw' fname = 'test_rw'
@ -341,9 +342,6 @@ def test_dobsio():
os.remove(fname + '.xml.gz') os.remove(fname + '.xml.gz')
[o.gamma_method() for o in rl] [o.gamma_method() for o in rl]
for i in range(len(ol)):
assert (ol[i] - rl[i].is_zero())
for i in range(len(ol)): for i in range(len(ol)):
if isinstance(ol[i], pe.Obs): if isinstance(ol[i], pe.Obs):
o = ol[i] - rl[i] o = ol[i] - rl[i]
@ -353,6 +351,9 @@ def test_dobsio():
for j in range(len(or1)): for j in range(len(or1)):
o = or1[j] - or2[j] o = or1[j] - or2[j]
assert(o.is_zero()) assert(o.is_zero())
if isinstance(ol[i], pe.Obs):
for name in ol[i].r_values:
assert(np.isclose(ol[i].r_values[name], rl[i].r_values[name]))
def test_reconstruct_non_linear_r_obs(tmp_path): def test_reconstruct_non_linear_r_obs(tmp_path):