This is intended as an example to demonstrate the use of overloading in object oriented programming. This was written as a Jupyter notebook (aka IPython) in Python 3. To run in Python 2, simply rename the variables that have unicode names, and replace truediv with div.
While there are several nice Python libraries that support uncertainty (for example, the powerful uncertainties package and the related units and uncertainties package pint), they usually use standard error combination rules. For a beginning physics class, often 'maximum error' combination is used. Here, instead of using a standard deviation based error and using combination rules based on uncorrelated statistical distributions, we assume a simple maximum error and simply add errors.
To implement this, let's build a Python class and use overloading to implement algebraic operations.
import unittest
import math
The main rules are listed below.
If $C=A+B$ or $C=A-B$, then $$ \delta C = \delta A + \delta B. \tag{1} $$
If $C=AB$ or $C=A/B$, then $$ \delta C = C_0 \left( \frac{\delta A}{A_0} + \frac{\delta B}{B_0} \right). \tag{2} $$
Following from that, if $C=A^n$ then $$ \delta C = A_0^n \left( n \frac{\delta A}{A_0} \right). \tag{3} $$
class LabUnc(object):
@staticmethod
def combine(a, b):
return a + b
rounding_rule = 1
'This is the number to round at for display, lab rule is 1, particle physics uses 3.54'
def __init__(self, number, uncertainty=0):
self.n = number
self.s = abs(uncertainty)
@property
def ndigits(self):
v = math.ceil(-math.log10(self.s) + math.log10(self.rounding_rule))
return v if v > 0 else 0
@property
def max(self):
return self.n + self.s
@property
def min(self):
return self.n - self.s
def __repr__(self):
return '{0}({1.n}, {1.s})'.format(type(self).__name__, self)
def __str__(self):
return (format(self.n,'0.'+str(self.ndigits) + 'f')
+ ' ± ' + format(self.s,'0.'+str(self.ndigits) + 'f'))
def _repr_html_(self):
return str(self)
def __eq__(self,other):
return abs(self.n - other.n)<.0000001 and abs(self.s - other.s)<.0000001
def __add__(self, other): # (1)
return self.__class__(self.n+other.n, self.combine(self.s, other.s))
def __sub__(self, other): # (1)
return self.__class__(self.n-other.n, self.combine(self.s, other.s))
def __mul__(self, other): # (2)
C = self.n * other.n
δC = C * self.combine(self.s/self.n, other.s/other.n)
return type(self)(C, δC)
def __truediv__(self, other): # (2)
C = self.n / other.n
δC = C * self.combine(self.s/self.n, other.s/other.n)
return self.__class__(C, δC)
def __pow__(self, power): # (3)
C = self.n**power
δC = C * (power * self.s/self.n)
return self.__class__(C, δC)
Now that we have a (hopefully!) working class, let's put together a quick unittest to see if it does what we expect.
class TestLabUncCombine(unittest.TestCase):
def testAdd(self):
self.assertEqual(LabUnc(10,1) + LabUnc(10,2), LabUnc(20,3))
self.assertEqual(LabUnc(20,.1) + LabUnc(10,.1), LabUnc(30,.2))
def testSub(self):
self.assertEqual(LabUnc(10,1) - LabUnc(10,2), LabUnc(0,3))
self.assertEqual(LabUnc(20,.1) - LabUnc(10,.1), LabUnc(10,.2))
def testMul(self):
self.assertEqual(LabUnc(10,1) * LabUnc(10,2), LabUnc(100,30))
self.assertEqual(LabUnc(20,.1) * LabUnc(10,.1), LabUnc(200,3))
def testTrueDiv(self):
self.assertEqual(LabUnc(10,1) / LabUnc(10,2), LabUnc(1,.3))
self.assertEqual(LabUnc(20,.1) / LabUnc(10,.1), LabUnc(2,.03))
def testPow(self):
self.assertEqual(LabUnc(10,1)**2, LabUnc(100,20))
self.assertEqual(LabUnc(20,.1)**3, LabUnc(8000,120))
Since we are putting our TestCase in a IPython notebook, we'll need to run it manually. This would be simpler if it was in a separate testing file.
suite = unittest.TestLoader().loadTestsFromTestCase(TestLabUncCombine)
unittest.TextTestRunner(verbosity=2).run(suite);
We can now say that our class works! Now, let's do some example calculations.
A = LabUnc(3, .005)
B = LabUnc(1.11, .05)
C = LabUnc(.004, .001)
D = LabUnc(2.02, .08)
print('A × B =', A * B)
print('A × C =', A * C)
print('A / B =', A / B)
print('A / D =', A / D)
print('A⁵ =', A**5)
a = LabUnc(12,2)
b = LabUnc(3,.5)
c = LabUnc(2,.2)
d = (a/b) - c
print('d₀ =', d.n)
print('d_max =', a.max / b.min - c.min)
print('d_min =', a.min / b.max - c.max)
print('δd =', d.s)
print('d₀ + δd =', d.max)
print('d₀ - δd =', d.min)
Bonus: inheritance¶
As a quick example of inheritance, let's set up the traditional uncertainty method using our class and inheritance. You may have noticed that the previous class was built with a combine method rather than simply adding uncertainties. Let's use the standard deviation rules to replace this combine with the one from traditional error analysis.
If $C=A+B$ or $C=A-B$, then $$ \delta C = \sqrt{\delta A^2 + \delta B^2}. \tag{1} $$
If $C=AB$ or $C=A/B$, then $$ \delta C = C_0 \sqrt{ \left( \frac{\delta A}{A_0}\right)^2 + \left(\frac{\delta B}{B_0}\right)^2 }. \tag{2} $$
class StdUnc(LabUnc):
@staticmethod
def combine(a, b):
return math.sqrt(a**2 + b**2)
Now, we can test this; since there are libraries that handle uncertainty this way, we can compare to those. (I'll be using slightly more advanced testing methods here.) Notice that I used n and s for the value and uncertainty; this was to allow duck typing for the comparison method for the uncertainties library classes.
from uncertainties import ufloat
import itertools
from operator import add, sub, mul, truediv
class TestStdUncCombine(unittest.TestCase):
def setUp(self):
cases = ((10,1),(10,2),(20,.1),(10,.1),(.1234,.02))
self.pairs = itertools.permutations(cases, 2)
def run_operation_on_each(self, operation):
for a,b in self.pairs:
self.assertEqual(operation(StdUnc(*a),StdUnc(*b)),
operation(ufloat(*a),ufloat(*b)))
def testAdd(self):
self.run_operation_on_each(add)
def testSub(self):
self.run_operation_on_each(sub)
def testMul(self):
self.run_operation_on_each(mul)
def testTrueDiv(self):
self.run_operation_on_each(truediv)
suite = unittest.TestLoader().loadTestsFromTestCase(TestStdUncCombine)
unittest.TextTestRunner(verbosity=2).run(suite);
No comments:
Post a Comment