[Feat] Meas class for gaussian error propagation.

This commit is contained in:
Fabian Joswig 2024-11-17 12:18:49 +01:00
parent 0ce765a99d
commit f2a98bb5bd
2 changed files with 70 additions and 2 deletions

View file

@ -9,6 +9,8 @@ import matplotlib.pyplot as plt
from scipy.stats import skew, skewtest, kurtosis, kurtosistest from scipy.stats import skew, skewtest, kurtosis, kurtosistest
import numdifftools as nd import numdifftools as nd
from itertools import groupby from itertools import groupby
from typing import Optional, Union
import uuid
from .covobs import Covobs from .covobs import Covobs
# Improve print output of numpy.ndarrays containing Obs objects. # Improve print output of numpy.ndarrays containing Obs objects.
@ -1354,6 +1356,9 @@ def derived_observable(func, data, array_mode=False, **kwargs):
final_result[i_val]._value = new_val final_result[i_val]._value = new_val
final_result[i_val].reweighted = reweighted final_result[i_val].reweighted = reweighted
if not final_result[i_val].idl:
final_result[i_val].gm()
if multi == 0: if multi == 0:
final_result = final_result.item() final_result = final_result.item()
@ -1810,7 +1815,7 @@ def cov_Obs(means, cov, name, grad=None):
co : Covobs co : Covobs
Covobs to be embedded into the Obs Covobs to be embedded into the Obs
""" """
o = Obs([], [], means=[]) o = Obs(samples=[], names=[], means=[])
o._value = co.value o._value = co.value
o.names.append(co.name) o.names.append(co.name)
o._covobs[co.name] = co o._covobs[co.name] = co
@ -1822,7 +1827,7 @@ def cov_Obs(means, cov, name, grad=None):
means = [means] means = [means]
for i in range(len(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): if ol[0].covobs[name].N != len(means):
raise Exception('You have to provide %d mean values!' % (ol[0].N)) raise Exception('You have to provide %d mean values!' % (ol[0].N))
if len(ol) == 1: if len(ol) == 1:
@ -1830,6 +1835,43 @@ def cov_Obs(means, cov, name, grad=None):
return ol 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
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): def _determine_gap(o, e_content, e_name):
gaps = [] gaps = []
for r_name in e_content[e_name]: for r_name in e_content[e_name]:

View file

@ -1458,3 +1458,29 @@ def test_missing_replica():
for op in [[O1O2, O1O2b], [O1O2O3, O1O2O3b]]: for op in [[O1O2, O1O2b], [O1O2O3, O1O2O3b]]:
assert np.isclose(op[1].value, op[0].value) assert np.isclose(op[1].value, op[0].value)
assert np.isclose(op[1].dvalue, op[0].dvalue, atol=0, rtol=5e-2) 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