source: sasview/src/sas/sascalc/data_util/uncertainty.py @ 41d6187

Last change on this file since 41d6187 was 574adc7, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

convert sascalc to python 2/3 syntax
  • Property mode set to 100644
File size: 11.1 KB
Line 
1r"""
2Uncertainty propagation class for arithmetic, log and exp.
3
4Based on scalars or numpy vectors, this class allows you to store and
5manipulate values+uncertainties, with propagation of gaussian error for
6addition, subtraction, multiplication, division, power, exp and log.
7
8Storage properties are determined by the numbers used to set the value
9and uncertainty.  Be sure to use floating point uncertainty vectors
10for inplace operations since numpy does not do automatic type conversion.
11Normal operations can use mixed integer and floating point.  In place
12operations such as *a \*= b* create at most one extra copy for each operation.
13By contrast, *c = a\*b* uses four intermediate vectors, so shouldn't be used
14for huge arrays.
15"""
16
17from __future__ import division
18
19import numpy as np
20
21from .import err1d
22from .formatnum import format_uncertainty
23
24__all__ = ['Uncertainty']
25
26# TODO: rename to Measurement and add support for units?
27# TODO: C implementation of *,/,**?
28class Uncertainty(object):
29    # Make standard deviation available
30    def _getdx(self): return np.sqrt(self.variance)
31    def _setdx(self,dx):
32        # Direct operation
33        #    variance = dx**2
34        # Indirect operation to avoid temporaries
35        self.variance[:] = dx
36        self.variance **= 2
37    dx = property(_getdx,_setdx,doc="standard deviation")
38
39    # Constructor
40    def __init__(self, x, variance=None):
41        self.x, self.variance = x, variance
42
43    # Numpy array slicing operations
44    def __len__(self):
45        return len(self.x)
46    def __getitem__(self,key):
47        return Uncertainty(self.x[key],self.variance[key])
48    def __setitem__(self,key,value):
49        self.x[key] = value.x
50        self.variance[key] = value.variance
51    def __delitem__(self, key):
52        del self.x[key]
53        del self.variance[key]
54    #def __iter__(self): pass # Not sure we need iter
55
56    # Normal operations: may be of mixed type
57    def __add__(self, other):
58        if isinstance(other,Uncertainty):
59            return Uncertainty(*err1d.add(self.x,self.variance,other.x,other.variance))
60        else:
61            return Uncertainty(self.x+other, self.variance+0) # Force copy
62    def __sub__(self, other):
63        if isinstance(other,Uncertainty):
64            return Uncertainty(*err1d.sub(self.x,self.variance,other.x,other.variance))
65        else:
66            return Uncertainty(self.x-other, self.variance+0) # Force copy
67    def __mul__(self, other):
68        if isinstance(other,Uncertainty):
69            return Uncertainty(*err1d.mul(self.x,self.variance,other.x,other.variance))
70        else:
71            return Uncertainty(self.x*other, self.variance*other**2)
72    def __truediv__(self, other):
73        if isinstance(other,Uncertainty):
74            return Uncertainty(*err1d.div(self.x,self.variance,other.x,other.variance))
75        else:
76            return Uncertainty(self.x/other, self.variance/other**2)
77    def __pow__(self, other):
78        if isinstance(other,Uncertainty):
79            # Haven't calcuated variance in (a+/-da) ** (b+/-db)
80            return NotImplemented
81        else:
82            return Uncertainty(*err1d.pow(self.x,self.variance,other))
83
84    # Reverse operations
85    def __radd__(self, other):
86        return Uncertainty(self.x+other, self.variance+0) # Force copy
87    def __rsub__(self, other):
88        return Uncertainty(other-self.x, self.variance+0)
89    def __rmul__(self, other):
90        return Uncertainty(self.x*other, self.variance*other**2)
91    def __rtruediv__(self, other):
92        x,variance = err1d.pow(self.x,self.variance,-1)
93        return Uncertainty(x*other,variance*other**2)
94    def __rpow__(self, other): return NotImplemented
95
96    # In-place operations: may be of mixed type
97    def __iadd__(self, other):
98        if isinstance(other,Uncertainty):
99            self.x,self.variance \
100                = err1d.add_inplace(self.x,self.variance,other.x,other.variance)
101        else:
102            self.x+=other
103        return self
104    def __isub__(self, other):
105        if isinstance(other,Uncertainty):
106            self.x,self.variance \
107                = err1d.sub_inplace(self.x,self.variance,other.x,other.variance)
108        else:
109            self.x-=other
110        return self
111    def __imul__(self, other):
112        if isinstance(other,Uncertainty):
113            self.x, self.variance \
114                = err1d.mul_inplace(self.x,self.variance,other.x,other.variance)
115        else:
116            self.x *= other
117            self.variance *= other**2
118        return self
119    def __itruediv__(self, other):
120        if isinstance(other,Uncertainty):
121            self.x,self.variance \
122                = err1d.div_inplace(self.x,self.variance,other.x,other.variance)
123        else:
124            self.x /= other
125            self.variance /= other**2
126        return self
127    def __ipow__(self, other):
128        if isinstance(other,Uncertainty):
129            # Haven't calcuated variance in (a+/-da) ** (b+/-db)
130            return NotImplemented
131        else:
132            self.x,self.variance = err1d.pow_inplace(self.x, self.variance, other)
133        return self
134
135    # Use true division instead of integer division
136    def __div__(self, other): return self.__truediv__(other)
137    def __rdiv__(self, other): return self.__rtruediv__(other)
138    def __idiv__(self, other): return self.__itruediv__(other)
139
140
141    # Unary ops
142    def __neg__(self):
143        return Uncertainty(-self.x,self.variance)
144    def __pos__(self):
145        return self
146    def __abs__(self):
147        return Uncertainty(np.abs(self.x),self.variance)
148
149    def __str__(self):
150        #return str(self.x)+" +/- "+str(np.sqrt(self.variance))
151        if np.isscalar(self.x):
152            return format_uncertainty(self.x,np.sqrt(self.variance))
153        else:
154            return [format_uncertainty(v,dv)
155                    for v,dv in zip(self.x,np.sqrt(self.variance))]
156    def __repr__(self):
157        return "Uncertainty(%s,%s)"%(str(self.x),str(self.variance))
158
159    # Not implemented
160    def __floordiv__(self, other): return NotImplemented
161    def __mod__(self, other): return NotImplemented
162    def __divmod__(self, other): return NotImplemented
163    def __mod__(self, other): return NotImplemented
164    def __lshift__(self, other): return NotImplemented
165    def __rshift__(self, other): return NotImplemented
166    def __and__(self, other): return NotImplemented
167    def __xor__(self, other): return NotImplemented
168    def __or__(self, other): return NotImplemented
169
170    def __rfloordiv__(self, other): return NotImplemented
171    def __rmod__(self, other): return NotImplemented
172    def __rdivmod__(self, other): return NotImplemented
173    def __rmod__(self, other): return NotImplemented
174    def __rlshift__(self, other): return NotImplemented
175    def __rrshift__(self, other): return NotImplemented
176    def __rand__(self, other): return NotImplemented
177    def __rxor__(self, other): return NotImplemented
178    def __ror__(self, other): return NotImplemented
179
180    def __ifloordiv__(self, other): return NotImplemented
181    def __imod__(self, other): return NotImplemented
182    def __idivmod__(self, other): return NotImplemented
183    def __imod__(self, other): return NotImplemented
184    def __ilshift__(self, other): return NotImplemented
185    def __irshift__(self, other): return NotImplemented
186    def __iand__(self, other): return NotImplemented
187    def __ixor__(self, other): return NotImplemented
188    def __ior__(self, other): return NotImplemented
189
190    def __invert__(self): return NotImplmented  # For ~x
191    def __complex__(self): return NotImplmented
192    def __int__(self): return NotImplmented
193    def __long__(self): return NotImplmented
194    def __float__(self): return NotImplmented
195    def __oct__(self): return NotImplmented
196    def __hex__(self): return NotImplmented
197    def __index__(self): return NotImplmented
198    def __coerce__(self): return NotImplmented
199
200    def log(self):
201        return Uncertainty(*err1d.log(self.x,self.variance))
202
203    def exp(self):
204        return Uncertainty(*err1d.exp(self.x,self.variance))
205
206def log(val): return self.log()
207def exp(val): return self.exp()
208
209def test():
210    a = Uncertainty(5,3)
211    b = Uncertainty(4,2)
212
213    # Scalar operations
214    z = a+4
215    assert z.x == 5+4 and z.variance == 3
216    z = a-4
217    assert z.x == 5-4 and z.variance == 3
218    z = a*4
219    assert z.x == 5*4 and z.variance == 3*4**2
220    z = a/4
221    assert z.x == 5./4 and z.variance == 3./4**2
222
223    # Reverse scalar operations
224    z = 4+a
225    assert z.x == 4+5 and z.variance == 3
226    z = 4-a
227    assert z.x == 4-5 and z.variance == 3
228    z = 4*a
229    assert z.x == 4*5 and z.variance == 3*4**2
230    z = 4/a
231    assert z.x == 4./5 and abs(z.variance - 3./5**4 * 4**2) < 1e-15
232
233    # Power operations
234    z = a**2
235    assert z.x == 5**2 and z.variance == 4*3*5**2
236    z = a**1
237    assert z.x == 5**1 and z.variance == 3
238    z = a**0
239    assert z.x == 5**0 and z.variance == 0
240    z = a**-1
241    assert z.x == 5**-1 and abs(z.variance - 3./5**4) < 1e-15
242
243    # Binary operations
244    z = a+b
245    assert z.x == 5+4 and z.variance == 3+2
246    z = a-b
247    assert z.x == 5-4 and z.variance == 3+2
248    z = a*b
249    assert z.x == 5*4 and z.variance == (5**2*2 + 4**2*3)
250    z = a/b
251    assert z.x == 5./4 and abs(z.variance - (3./5**2 + 2./4**2)*(5./4)**2) < 1e-15
252
253    # ===== Inplace operations =====
254    # Scalar operations
255    y = a+0; y += 4
256    z = a+4
257    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
258    y = a+0; y -= 4
259    z = a-4
260    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
261    y = a+0; y *= 4
262    z = a*4
263    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
264    y = a+0; y /= 4
265    z = a/4
266    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
267
268    # Power operations
269    y = a+0; y **= 4
270    z = a**4
271    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
272
273    # Binary operations
274    y = a+0; y += b
275    z = a+b
276    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
277    y = a+0; y -= b
278    z = a-b
279    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
280    y = a+0; y *= b
281    z = a*b
282    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
283    y = a+0; y /= b
284    z = a/b
285    assert y.x == z.x and abs(y.variance-z.variance) < 1e-15
286
287
288    # =============== vector operations ================
289    # Slicing
290    z = Uncertainty(np.array([1,2,3,4,5]),np.array([2,1,2,3,2]))
291    assert z[2].x == 3 and z[2].variance == 2
292    assert (z[2:4].x == [3,4]).all()
293    assert (z[2:4].variance == [2,3]).all()
294    z[2:4] = Uncertainty(np.array([8,7]),np.array([4,5]))
295    assert z[2].x == 8 and z[2].variance == 4
296    A = Uncertainty(np.array([a.x]*2),np.array([a.variance]*2))
297    B = Uncertainty(np.array([b.x]*2),np.array([b.variance]*2))
298
299    # TODO complete tests of copy and inplace operations for vectors and slices.
300
301    # Binary operations
302    z = A+B
303    assert (z.x == 5+4).all() and (z.variance == 3+2).all()
304    z = A-B
305    assert (z.x == 5-4).all() and (z.variance == 3+2).all()
306    z = A*B
307    assert (z.x == 5*4).all() and (z.variance == (5**2*2 + 4**2*3)).all()
308    z = A/B
309    assert (z.x == 5./4).all()
310    assert (abs(z.variance - (3./5**2 + 2./4**2)*(5./4)**2) < 1e-15).all()
311
312    # printing; note that sqrt(3) ~ 1.7
313    assert str(Uncertainty(5,3)) == "5.0(17)"
314    assert str(Uncertainty(15,3)) == "15.0(17)"
315    assert str(Uncertainty(151.23356,0.324185**2)) == "151.23(32)"
316
317if __name__ == "__main__": test()
Note: See TracBrowser for help on using the repository browser.