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

costrafo411
Last change on this file since ba91f71 was b699768, checked in by Piotr Rozyczko <piotr.rozyczko@…>, 9 years ago

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