-
Notifications
You must be signed in to change notification settings - Fork 336
/
tl_uval.py
238 lines (201 loc) · 7.33 KB
/
tl_uval.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
# Copyright (c) 2018 Technical University of Munich
# Author: Martin Becker
#
# This program is free software; you can redistribute it and/or modify it
# under the terms and conditions of the GNU General Public License,
# version 2, as published by the Free Software Foundation.
#
# This program is distributed in the hope it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
# more details.
#
import math
import logging
import operator
from tl_io import warn
log = logging.getLogger(__name__)
TEMPVAL = 'anon'
div_op = operator.div if 'div' in operator.__dict__ else None # type: ignore
def combine_uval(ulist):
"""
Combine multiple measurements of the same event into one measurement.
Uses weighted average.
"""
combined = None
if ulist is not None:
combined = ulist[0]
for oth in ulist[1:]:
combined.update(oth)
return combined
class UVal:
"""
Measurement value annotated with uncertainty. Supports binary operators for error propagation.
"""
def __init__(self, name, value, stddev=0., samples=1, mux=100., comment="", computed=False):
self.name = name
self.comment = comment
self.value = value
self.stddev = stddev
self.samples = samples
self.computed = computed
self.is_ratio = False
self.multiplex = mux
def __repr__(self):
return "{} [{} +- {}]*{}".format(self.name, self.value, self.stddev, self.samples)
def format_value(self, unit):
if self.value is None:
return ""
if self.is_ratio:
return "{:>16.1f} ".format(self.value * 100.)
elif unit == "Count" or unit == "Clocks":
return "{:16,.0f} ".format(self.value)
elif self.value > 1000:
return "{:16,.1f} ".format(self.value)
else:
return "{:16.2f}".format(self.value)
def format_value_raw(self):
if self.value is None:
return ""
if self.is_ratio:
return "{:>13.1f}".format(self.value * 100.)
else:
return "{:13.1f}".format(self.value)
def format_uncertainty(self):
"""string representation of measurement uncertainty"""
vs = ""
if self.stddev is not None:
if self.is_ratio:
if self.value != 0.:
v = self.stddev * 100.
else:
v = 0.
vs += "{:.1f}".format(v)
else:
vs += "{:6,.1f}".format(self.stddev)
return vs
def format_mux(self):
vs = ""
if self.multiplex and self.multiplex == self.multiplex:
vs = "[{:4.1f}%]".format(self.multiplex)
return vs
@staticmethod
def _merge_mux(lhs, rhs):
return min(lhs.multiplex, rhs.multiplex)
def update(self, other):
"""merge data from other event into this"""
assert isinstance(other, UVal), "wrong type"
# --
# calc weighted average
n = self.samples + other.samples
res = (1. / n) * (self.samples * self + other.samples * other)
# apply 'res' to this
self.samples = n
self.value = res.value
self.stddev = res.stddev
self.multiplex = UVal._merge_mux(self, other)
######################
# operators
######################
def ensure_uval(binop): # type: ignore
"""decorator to ensure binary operators are both UVals"""
def wrapper(self, v): # type: ignore
if isinstance(v, UVal):
return binop(self, v) # type: ignore
elif isinstance(v, (float, int)):
return binop(self, UVal(TEMPVAL, value=v, stddev=0)) # type: ignore
else:
return NotImplemented
return wrapper
@ensure_uval # type: ignore
def __sub__(self, other):
return UVal._calc(operator.sub, self, other)
@ensure_uval # type: ignore
def __add__(self, other):
return UVal._calc(operator.add, self, other)
@ensure_uval # type: ignore
def __mul__(self, other):
return UVal._calc(operator.mul, self, other)
@ensure_uval # type: ignore
def __div__(self, other):
return UVal._calc(operator.div, self, other) # type: ignore
@ensure_uval # type: ignore
def __truediv__(self, other):
return UVal._calc(operator.truediv, self, other)
@ensure_uval # type: ignore
def __lt__(self, other):
return self.value < other.value
@ensure_uval # type: ignore
def __le__(self, other):
return self.value <= other.value
@ensure_uval # type: ignore
def __eq__(self, other):
return self.value == other.value
@ensure_uval # type: ignore
def __ne__(self, other):
return not self.__eq__(other)
@ensure_uval # type: ignore
def __ge__(self, other):
return self.value >= other.value
@ensure_uval # type: ignore
def __gt__(self, other):
return self.value > other.value
@ensure_uval # type: ignore
def __rsub__(self, other):
"""other - self"""
return UVal._calc(operator.sub, other, self)
@ensure_uval # type: ignore
def __radd__(self, other):
"""other + self"""
return UVal._calc(operator.add, other, self)
@ensure_uval # type: ignore
def __rmul__(self, other):
"""other * self"""
return UVal._calc(operator.mul, other, self)
@ensure_uval # type: ignore
def __rdiv__(self, other):
"""other / self"""
return UVal._calc(operator.div, other, self) # type: ignore
@ensure_uval # type: ignore
def __rtruediv__(self, other):
"""other / self"""
return UVal._calc(operator.truediv, other, self)
def __nonzero__(self): # python 2
return self.value != 0.0
def __bool__(self): # python 3
return self.value != 0.0
#########################
# uncertainty propagator
#########################
@staticmethod
def _calc(op, lhs, rhs, cov=0.):
"""Compute the result of 'lhs [op] rhs' and propagate standard deviations"""
A = lhs.value
B = rhs.value
a = lhs.stddev
b = rhs.stddev
# new value
f = op(float(A), B)
if isinstance(f, float) and f.is_integer():
f = int(f)
# uncertainty
if op in (operator.mul, operator.truediv, div_op):
sgn = 1 if op == operator.mul else -1
if A != 0 and B != 0:
u = abs(f) * math.sqrt(pow(float(a)/A, 2) + pow(float(b)/B, 2) + sgn*2.*cov/(A*B))
elif op == operator.mul:
u = 0.
elif op == div_op or op == operator.truediv:
u = 0.
if A != 0:
warn("Error prop failed because of DIV/0: {} {} {}".format(lhs, op, rhs))
elif op in (operator.add, operator.sub):
sgn = 1 if op == operator.add else -1
u = math.sqrt(pow(a, 2) + pow(b, 2) + sgn*2.*cov)
else:
u = None
log.error("Unsupported operation for uncertainty propagator in {} {} {}".format(lhs, op, rhs))
# --
ret = UVal(TEMPVAL, value=f, stddev=u, mux=UVal._merge_mux(lhs, rhs), computed=True)
log.debug("{} {} {} => {}", lhs, op, rhs, ret)
return ret