Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Feat] Meas class for gaussian error propagation. #245

Open
wants to merge 3 commits into
base: develop
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions pyerrors/input/json.py
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,8 @@ def _gen_data_d_from_list(ol):
def _gen_cdata_d_from_list(ol):
dl = []
for name in ol[0].cov_names:
if _is_uuid4_hex(name):
raise ValueError("Cannot safely serialize an Obs derived from a Meas object with a uuid as name. Consider recreating the Meas with an explict name.")
ed = {}
ed['id'] = name
ed['layout'] = str(ol[0].covobs[name].cov.shape).lstrip('(').rstrip(')').rstrip(',')
Expand Down Expand Up @@ -216,6 +218,11 @@ def _jsonifier(obj):
return json.dumps(d, indent=indent, ensure_ascii=False, default=_jsonifier, write_mode=json.WM_COMPACT)


def _is_uuid4_hex(s):
uuid4_hex_pattern = re.compile(r'^[0-9a-f]{32}$')
return bool(uuid4_hex_pattern.match(s))


def dump_to_json(ol, fname, description='', indent=1, gz=True):
"""Export a list of Obs or structures containing Obs to a .json(.gz) file.
Dict keys that are not JSON-serializable such as floats are converted to strings.
Expand Down
46 changes: 44 additions & 2 deletions pyerrors/obs.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@
from scipy.stats import skew, skewtest, kurtosis, kurtosistest
import numdifftools as nd
from itertools import groupby
from typing import Optional, Union
import uuid
from .covobs import Covobs

# Improve print output of numpy.ndarrays containing Obs objects.
Expand Down Expand Up @@ -1354,6 +1356,9 @@ def __init__(self, N):
final_result[i_val]._value = new_val
final_result[i_val].reweighted = reweighted

if not final_result[i_val].idl:
final_result[i_val].gm()

if multi == 0:
final_result = final_result.item()

Expand Down Expand Up @@ -1810,7 +1815,7 @@ def covobs_to_obs(co):
co : Covobs
Covobs to be embedded into the Obs
"""
o = Obs([], [], means=[])
o = Obs(samples=[], names=[], means=[])
o._value = co.value
o.names.append(co.name)
o._covobs[co.name] = co
Expand All @@ -1822,14 +1827,51 @@ def covobs_to_obs(co):
means = [means]

for i in range(len(means)):
ol.append(covobs_to_obs(Covobs(means[i], cov, name, pos=i, grad=grad)))
ol.append(covobs_to_obs(Covobs(float(means[i]), cov, name, pos=i, grad=grad)))
if ol[0].covobs[name].N != len(means):
raise Exception('You have to provide %d mean values!' % (ol[0].N))
if len(ol) == 1:
return ol[0]
return ol


class Meas(Obs):
"""Class for a scalar measurement.

Convenience wrapper for scalar measurements.
"""

def __init__(self, value: Union[float, int], dvalue: Union[float, int], name: Optional[str] = None):
""" Initialize Meas object.

Parameters
----------
value : float
Mean value of the measurement.
dvalue : list or array
Standard error of the measurement.
name : Optional[str]
Optional name identifier for the measurement. If none is specified, a random uuid
string is used instead.
"""
if not isinstance(value, (float, int)):
raise TypeError(f"value has to be a flaot or int, not {type(value)}")
if not isinstance(dvalue, (float, int)):
raise TypeError(f"dvalue has to be a float or int, not {type(dvalue)}")
super().__init__(samples=[], names=[], means=[])
if name is None:
name = uuid.uuid4().hex
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hey @fjosw,
the class looks good to me. Thanks for the idea and addition. I just wanted to ask if this is something that we want. I do not mind it, and if there is a specific use case for you, we can, of course, leave it in.

else:
if not isinstance(name, str):
raise TypeError(f"name has to be a str, not {type(name)}")

co = Covobs(float(value), float(dvalue) ** 2, name)
self._value = co.value
self.names.append(co.name)
self._covobs[co.name] = co
self._dvalue = np.sqrt(co.errsq())


def _determine_gap(o, e_content, e_name):
gaps = []
for r_name in e_content[e_name]:
Expand Down
11 changes: 11 additions & 0 deletions tests/json_io_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -409,3 +409,14 @@ def assert_equal_Obs(to, ro):
print(kw, "does not match.")
return False
return True


def test_meas_uuid(tmp_path):
meas = pe.Meas(0.3, 0.2)

for obs in [meas, meas + 1, meas * pe.pseudo_Obs(0.1, 0.001, "obs|r1")]:
with pytest.raises(ValueError):
jsonio.dump_to_json([obs], "test_file", indent=1, description='[This file should not be writable]')

name_meas = pe.Meas(0.3, 0.2, name="my name")
jsonio.dump_to_json([name_meas], "test_file", indent=1, description='[This file should be writable]')
34 changes: 34 additions & 0 deletions tests/obs_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -1458,3 +1458,37 @@ def test_missing_replica():
for op in [[O1O2, O1O2b], [O1O2O3, O1O2O3b]]:
assert np.isclose(op[1].value, op[0].value)
assert np.isclose(op[1].dvalue, op[0].dvalue, atol=0, rtol=5e-2)


def test_meas():
meas1 = pe.Meas(1.0, 0.1)
meas2 = pe.Meas(2, 1)

assert meas1 + meas2 == meas2 + meas1
assert meas1 - meas2 == -(meas2 - meas1)
assert meas1 * meas2 == meas2 * meas1

identical_sum = meas1 + meas1
assert identical_sum == 2 * meas1
assert np.isclose(identical_sum.dvalue, 2 * meas1.dvalue)

meas1_new = pe.Meas(1.0, 0.1)
not_identical_sum = meas1 + meas1_new
assert not_identical_sum.value == (2 * meas1).value
assert not_identical_sum != 2 * meas1
assert np.isclose(not_identical_sum.dvalue, np.sqrt(meas1.dvalue ** 2 + meas1_new.dvalue ** 2))

assert meas2 * meas2 == meas2 ** 2


def test_square_cov_obs():
cov = pe.cov_Obs(1, 0.1 ** 2, "testing")
cov2 = cov ** 2


def test_covobs_equal_meas():
value = 1.1
standard_error = 0.23
meas = pe.Meas(value, standard_error)
covo = pe.cov_Obs(value, standard_error ** 2, meas.names[0])
assert covo == meas