Source Code for Module reflectometry.reduction.uncertainty

  1  """ 
  2   
  3  Uncertainty propagation class, and log() and exp() functions. 
  4   
  5  Based on scalars or numpy vectors, this class allows you to store and 
  6  manipulate values+uncertainties, with propagation of gaussian error for 
  7  addition, subtraction, multiplication, division, power, exp() and log(). 
  8   
  9  Storage properties are determined by the numbers used to set the value 
 10  and uncertainty.  Be sure to use floating point uncertainty vectors 
 11  for inplace operations since numpy does not do automatic type conversion. 
 12  Normal operations can use mixed integer and floating point.  In place 
 13  operations (a *= b, etc.) create at most one extra copy for each operation. 
 14  c = a*b by contrast uses four intermediate vectors, so shouldn't be used 
 15  for huge arrays. 
 16  """ 
 17   
 18  from __future__ import division 
 19   
 20  import numpy 
 21  import err1d 
 22  from formatnum import format_uncertainty 
 23   
 24  __all__ = ['Uncertainty'] 
 25   
 26  # TODO: rename to Measurement and add support for units? 
 27  # TODO: C implementation of *,/,**? 
 28 -class Uncertainty(object): 
 29      # Make standard deviation available 
 30 -    def _getdx(self): return numpy.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(numpy.abs(self.x),self.variance) 
148   
149 -    def __str__(self): 
150          #return str(self.x)+" +/- "+str(numpy.sqrt(self.variance)) 
151          if numpy.isscalar(self.x): 
152              return format_uncertainty(self.x,numpy.sqrt(self.variance)) 
153          else: 
154              return [format_uncertainty(v,dv) 
155                      for v,dv in zip(self.x,numpy.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   
206 -def log(val): return self.log() 
207 -def exp(val): return self.exp() 
208   
209 -def 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(numpy.array([1,2,3,4,5]),numpy.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(numpy.array([8,7]),numpy.array([4,5])) 
295      assert z[2].x == 8 and z[2].variance == 4 
296      A = Uncertainty(numpy.array([a.x]*2),numpy.array([a.variance]*2)) 
297      B = Uncertainty(numpy.array([b.x]*2),numpy.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   
317  if __name__ == "__main__": test() 
318