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.
Theory
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} $$
Code
Now, we’ll make this a class (using Python 3 features, like unicode variables):
import math
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) < 0.0000001 and abs(self.s - other.s) < 0.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)
Note on attrs (click to expand)
If I were doing this today, I would probably use the amazing attrs library. That makes it easy and simple to make a class that holds data.
Tests
Now that we have a (hopefully!) working class, let’s put together a quick unittest to see if it does what we expect. I normal would use pytest, but I’ll stick to the standard library for this post, both for simplicity and because this started out as a Jupyter notebook:
import unittest
class TestLabUncCombine(unittest.TestCase):
def testAdd(self):
self.assertEqual(LabUnc(10, 1) + LabUnc(10, 2), LabUnc(20, 3))
self.assertEqual(LabUnc(20, 0.1) + LabUnc(10, 0.1), LabUnc(30, 0.2))
def testSub(self):
self.assertEqual(LabUnc(10, 1) - LabUnc(10, 2), LabUnc(0, 3))
self.assertEqual(LabUnc(20, 0.1) - LabUnc(10, 0.1), LabUnc(10, 0.2))
def testMul(self):
self.assertEqual(LabUnc(10, 1) * LabUnc(10, 2), LabUnc(100, 30))
self.assertEqual(LabUnc(20, 0.1) * LabUnc(10, 0.1), LabUnc(200, 3))
def testTrueDiv(self):
self.assertEqual(LabUnc(10, 1) / LabUnc(10, 2), LabUnc(1, 0.3))
self.assertEqual(LabUnc(20, 0.1) / LabUnc(10, 0.1), LabUnc(2, 0.03))
def testPow(self):
self.assertEqual(LabUnc(10, 1) ** 2, LabUnc(100, 20))
self.assertEqual(LabUnc(20, 0.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)
testAdd (__main__.TestLabUncCombine) ... ok
testMul (__main__.TestLabUncCombine) ... ok
testPow (__main__.TestLabUncCombine) ... ok
testSub (__main__.TestLabUncCombine) ... ok
testTrueDiv (__main__.TestLabUncCombine) ... ok
----------------------------------------------------------------------
Ran 5 tests in 0.010s
OK
Examples
We can now say that our class works! Now, let’s do some example calculations.
A = LabUnc(3, 0.005)
B = LabUnc(1.11, 0.05)
C = LabUnc(0.004, 0.001)
D = LabUnc(2.02, 0.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 × B = 3.3 ± 0.2
A × C = 0.012 ± 0.003
A / B = 2.7 ± 0.1
A / D = 1.49 ± 0.06
A⁵ = 243 ± 2
a = LabUnc(12, 2)
b = LabUnc(3, 0.5)
c = LabUnc(2, 0.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)
d₀ = 2.0
d_max = 3.8
d_min = 0.657142857142857
δd = 1.5333333333333332
print("d₀ + δd =", d.max)
print("d₀ - δd =", d.min)
d₀ + δd = 3.533333333333333
d₀ - δd = 0.4666666666666668
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.
Theory
If \(C=A+B\) or \(C=A-B\), then
$$ \delta C = \sqrt{\delta A^2 + \delta B^2}. \tag{4} $$
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{5} $$
Code
class StdUnc(LabUnc):
@staticmethod
def combine(a, b):
return math.sqrt(a**2 + b**2)
Tests
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, 0.1), (10, 0.1), (0.1234, 0.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)
testAdd (__main__.TestStdUncCombine) ... ok
testMul (__main__.TestStdUncCombine) ... ok
testSub (__main__.TestStdUncCombine) ... ok
testTrueDiv (__main__.TestStdUncCombine) ... ok
----------------------------------------------------------------------
Ran 4 tests in 0.049s
OK