Special methods¶

class MyClass:
    def __init__(self, a, b):
        self.a = a
        self.b = b

p1 = MyClass(2, 5)
p2 = MyClass(-1, 10)

p3 = p1 + p2
p4 = p1 - p2
p5 = p1*p2
p6 = p1**7 + 4*p3

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-1ce7ba7a902d> in <module>()
      7 p2 = MyClass(-1, 10)
      8 
----> 9 p3 = p1 + p2
     10 p4 = p1 - p2
     11 p5 = p1*p2

TypeError: unsupported operand type(s) for +: 'MyClass' and 'MyClass'

Special methods allow nice syntax and are recognized by double leading and trailing underscores¶

def __init__(self, ...)
def __call__(self, ...)
def __add__(self, other)

# Python syntax
y = Y(4)
print y(2)
z = Y(6)
print y + z

# What's actually going on
Y.__init__(y, 4)
print Y.__call__(y, 2)
Y.__init__(z, 6)
print Y.__add__(y, z)

We shall learn about many more such special methods

Example on a call special method¶

Replace the value method by a call special method:

class Y:
    def __init__(self, v0):
        self.v0 = v0
        self.g = 9.81

    def __call__(self, t):
        return self.v0*t - 0.5*self.g*t**2
    
    def value(self, t):
        return self.v0*t - 0.5*self.g*t**2

def func(t, v0=3, g=9.81):
    return v0*t - 0.5*g*t**2

Now we can write

t = 0.1
y = Y(3)
v = y.value(t)
print('v from y.value(t):', v)
v = y(t) # same as v = y.__call__(0.1) or Y.__call__(y, 0.1)
print('v from y.__call__(t):', v)

v from y.value(t): 0.25095
v from y.__call__(t): 0.25095

Note:

The instance y behaves and looks as a function!
The value(t) method does the same, but __call__ allows nicer syntax for computing function values

Representing a function by a class revisited¶

Given a function with $n+1$ parameters and one independent variable,

$$ f(x; p_0,\ldots,p_n) $$

it is wise to represent f by a class where $p_0,\ldots,p_n$ are attributes and __call__(x) computes $f(x)$

class MyFunc:
    def __init__(self, p0, p1, p2, ..., pn):
        self.p0 = p0
        self.p1 = p1
        ...
        self.pn = pn

    def __call__(self, x):
        return ...

Can we automatically differentiate a function?¶

Given some mathematical function in Python, say

def f(x):
    return x**3

can we make a class Derivative and write

dfdx = Derivative(f)

so that dfdx behaves as a function that computes the derivative of f(x)?

print(dfdx(2))   # computes 3*x**2 for x=2

12.000060000261213

Automagic differentiation; solution¶

Method.

We use numerical differentiation "behind the curtain":

$$ f'(x) \approx {f(x+h)-f(x)\over h} $$

for a small (yet moderate) $h$, say $h=10^{-5}$

Implementation.

class Derivative:
    def __init__(self, f, h=1e-5):
        self.f = f
        self.h = float(h)

    def __call__(self, x):
        f, h = self.f, self.h      # make short forms
        return (f(x+h) - f(x))/h

Automagic differentiation; demo¶

from math import pi, sin, cos
df = Derivative(sin)
x = pi
df(x)

-0.9999999999898844

cos(x)  # exact

-1.0

Another function¶

$g(t) = t^3$

$g'(t) = 3t^2$

$g'(1) = 3$

def g(t):
    return t**3

dg = Derivative(g)
t = 1
dg(t)  # compare with 3 (exact)

3.000030000110953

Automagic differentiation; useful in Newton's method¶

Newton's method solves nonlinear equations $f(x)=0$, but the method requires $f'(x)$

def Newton(f, xstart, dfdx, epsilon=1E-6):
    ...
    return x, no_of_iterations, f(x)

Suppose $f'(x)$ requires boring/lengthy derivation, then class Derivative is handy:

def f(x):
    return 100000*(x - 0.9)**2 * (x - 1.1)**3

df = Derivative(f)
xstart = 1.01
Newton(f, xstart, df, epsilon=1E-5)

Automagic differentiation; test function¶

How can we test class Derivative?
Method 1: compute $(f(x+h)-f(x))/h$ by hand for some $f$ and $h$
Method 2: utilize that linear functions are differentiated exactly by our numerical formula, regardless of $h$

Test function based on method 2:

def test_Derivative():
    # The formula is exact for linear functions, regardless of h
    #f = lambda x: a*x + b
    def f(x):
        return a*x + b
    a = 3.5; b = 8
    dfdx = Derivative(f, h=0.5)
    diff = abs(dfdx(4.5) - a)
    assert diff < 1E-14, 'bug in class Derivative, diff=%s' % diff
    
retur = test_Derivative()
if retur == None:
    print('Success!')

Success!

Automagic differentiation; explanation of the test function¶

Use of lambda functions:

f = lambda x: a*x + b

is equivalent to

def f(x):
    return a*x + b

Lambda functions are convenient for producing quick, short code

Use of closure:

f = lambda x: a*x + b
a = 3.5; b = 8
dfdx = Derivative(f, h=0.5)
dfdx(4.5)

Looks straightforward...but

How can Derivative.__call__ know a and b when it calls our f(x) function?
Local functions inside functions remember (have access to) all local variables in the function they are defined (!)
f can access a and b in test_Derivative even when called from __call__ in class `Derivative
f is known as a closure in computer science

Automagic integration; problem setting¶

Given a function $f(x)$, we want to compute

$$ F(x; a) = \int_a^x f(t)dt $$

Automagic integration; technique¶

$$ F(x; a) = \int_a^x f(t)dt $$

Technique: Midpoint rule or Trapezoidal rule, here the latter:

$$ \int_a^x f(t)dt = h\left({1\over2}f(a) + \sum_{i=1}^{n-1} f(a+ih) + {1\over2}f(x)\right) $$

Desired application code:

from math import exp, sin
def f(x):
    return exp(-x**2)*sin(10*x)

a = 0; n = 200
F = Integral(f, a, n)
x = 1.2
print(F(x))

0.08596973533274166

Automagic integration; implementation¶

def trapezoidal(f, a, x, n):
    h = (x-a)/float(n)
    I = 0.5*f(a)
    for i in range(1, n):
        I += f(a + i*h)
    I += 0.5*f(x)
    I *= h
    return I

Class Integral holds f, a and n as attributes and has a call special method for computing the integral:

class Integral:
    def __init__(self, f, a, n=100):
        self.f, self.a, self.n = f, a, n

    def __call__(self, x):
        return trapezoidal(self.f, self.a, x, self.n)

Automagic integration; test function¶

How can we test class Integral?
Method 1: compute by hand for some $f$ and small $n$
Method 2: utilize that linear functions are integrated exactly by our numerical formula, regardless of $n$

Test function based on method 2:

def test_Integral():
    f = lambda x: 2*x + 5
    F = lambda x: x**2 + 5*x - (a**2 + 5*a)
    a = 2
    I = Integral(f, a, n=4)
    x = 6
    diff = abs(I(x) - (F(x) - F(a)))
    assert diff < 1E-15, 'bug in class Integral, diff=%s' % diff
test_Integral()

Special method for printing¶

In Python, we can usually print an object a by print a, works for built-in types (strings, lists, floats, ...)
Python does not know how to print objects of a user-defined class, but if the class defines a method __str__, Python will use this method to convert an object to a string

Example:

class Y:
    def __init__(self, v0, g=9.81):
        self.v0 = v0
        self.g = g
    
    def __call__(self, t):
        return self.v0*t - 0.5*self.g*t**2

    def __str__(self):
        return 'v0*t - 0.5*g*t**2; v0=%g' % self.v0

Demo:

y_instance = Y(1.5)
y_instance(0.2)
print(str(y_instance))

v0*t - 0.5*g*t**2; v0=1.5

y

<__main__.Y at 0x10a7a95c0>

Class for polynomials; functionality¶

A polynomial can be specified by a list of its coefficients. For example, $1 - x^2 + 2x^3$ is

$$ 1\cdot x^0 + 0\cdot x^1 - 1\cdot x^2 + 2\cdot x^3 $$

and the coefficients can be stored as [1, 0, -1, 2]

Desired application code:

p1 = Polynomial([1, -1])
print(p1)

p2 = Polynomial([0, 1, 0, 0, -6, -1])
p3 = p1 + p2
print(p3.coeff)

print(p3)

p2.differentiate()
print(p2)

How can we make class Polynomial?

Class Polynomial; basic code¶

class Polynomial1:
    def __init__(self, coefficients):
        self.coeff = coefficients

    def __call__(self, x):
        s = 0
        # i = 0, 1, 2, 3, ...
        for i in range(len(self.coeff)):
            s += self.coeff[i]*x**i
        return s

Class Polynomial; addition¶

class Polynomial2(Polynomial1):
    # __init__...
    # ...
    def __add__(self, other):
        # return self + other
        #
        # self.coeff  = [x, x, x, x, x]
        # other.coeff = [x, x, x]
        # or
        # self.coeff = [x, x]
        # other.coeff = [x, x, x, x]
        #
        # start with the longest list and add in the other:
        if len(self.coeff) > len(other.coeff):
            coeffsum = self.coeff[:]  # copy!
            for i in range(len(other.coeff)):
                coeffsum[i] += other.coeff[i]
        else:
            coeffsum = other.coeff[:] # copy!
            for i in range(len(self.coeff)):
                coeffsum[i] += self.coeff[i]
        return Polynomial(coeffsum)

Class Polynomial; multiplication¶

Mathematics:

Multiplication of two general polynomials:

$$ \left(\sum_{i=0}^Mc_ix^i\right)\left(\sum_{j=0}^N d_jx^j\right) = \sum_{i=0}^M \sum_{j=0}^N c_id_j x^{i+j} $$

The coeff. corresponding to power $i+j$ is $c_i\cdot d_j$. The list r of coefficients of the result: r[i+j] = c[i]*d[j] (i and j running from 0 to $M$ and $N$, resp.)

Implementation:

class Polynomial3(Polynomial2):
    # __init__...
    # ...
    # __add__...
    def __mul__(self, other):
        M = len(self.coeff) - 1
        N = len(other.coeff) - 1
        coeff = [0]*(M+N+1)  # or zeros(M+N+1)
        for i in range(0, M+1):
            for j in range(0, N+1):
                coeff[i+j] += self.coeff[i]*other.coeff[j]
        return Polynomial(coeff)

Class Polynomial; differentation¶

Mathematics:

Rule for differentiating a general polynomial:

$$ {d\over dx}\sum_{i=0}^n c_ix^i = \sum_{i=1}^n ic_ix^{i-1} $$

If c is the list of coefficients, the derivative has a list of coefficients, dc, where dc[i-1] = i*c[i] for i running from 1 to the largest index in c. Note that dc has one element less than c.

Implementation:

class Polynomial4(Polynomial3):
    # ...
    def differentiate(self):    # change self
        for i in range(1, len(self.coeff)):
            self.coeff[i-1] = i*self.coeff[i]
        del self.coeff[-1]

    def derivative(self):       # return new polynomial
        dpdx = Polynomial(self.coeff[:])  # copy
        dpdx.differentiate()
        return dpdx

Class Polynomial; pretty print¶

class Polynomial(Polynomial4):
    # ...
    def __str__(self):
        s = ''
        for i in range(0, len(self.coeff)):
            if self.coeff[i] != 0:
                s += ' + %g*x^%d' % (self.coeff[i], i)
        # fix layout (lots of special cases):
        s = s.replace('+ -', '- ')
        s = s.replace(' 1*', ' ')
        s = s.replace('x^0', '1')
        s = s.replace('x^1 ', 'x ')
        s = s.replace('x^1', 'x')
        if s[0:3] == ' + ':  # remove initial +
            s = s[3:]
        if s[0:3] == ' - ':  # fix spaces for initial -
            s = '-' + s[3:]
        return s

Class for polynomials; usage¶

Consider

$$ p_1(x)= 1-x,\quad p_2(x)=x - 6x^4 - x^5 $$

and their sum

$$ p_3(x) = p_1(x) + p_2(x) = 1 -6x^4 - x^5 $$

p1 = Polynomial([1, -1])
print(p1)

1 - x

p2 = Polynomial([0, 1, 0, 0, -6, -1])
p3 = p1 + p2
print(p3.coeff)
print(p3)

[1, 0, 0, 0, -6, -1]
1 - 6*x^4 - x^5

print(p2)

x - 6*x^4 - x^5

p2.differentiate()
print(p2)

1 - 24*x^3 - 5*x^4

The programmer is in charge of defining special methods!¶

How should, e.g., __add__(self, other) be defined? This is completely up to the programmer, depending on what is meaningful by object1 + object2.

An anthropologist was asking a primitive tribesman about arithmetic. When the anthropologist asked, What does two and two make? the tribesman replied, Five. Asked to explain, the tribesman said, If I have a rope with two knots, and another rope with two knots, and I join the ropes together, then I have five knots.

Special methods for arithmetic operations¶

c = a + b    #  c = a.__add__(b)

c = a - b    #  c = a.__sub__(b)

c = a*b      #  c = a.__mul__(b)

c = a/b      #  c = a.__div__(b) # Python 2.x
c = a/b      #  c = a.__truediv__(b) # Python 3.x

c = a**e     #  c = a.__pow__(e)

Special methods for comparisons¶

a == b       #  a.__eq__(b)

a != b       #  a.__ne__(b)

a < b        #  a.__lt__(b)

a <= b       #  a.__le__(b)

a > b        #  a.__gt__(b)

a >= b       #  a.__ge__(b)

Class for vectors in the plane¶

Mathematical operations for vectors in the plane:

$$ \begin{align*} (a,b) + (c,d) &= (a+c, b+d)\\ (a,b) - (c,d) &= (a-c, b-d)\\ (a,b)\cdot(c,d) &= ac + bd\\ (a,b) &= (c, d)\hbox{ if }a=c\hbox{ and }b=d \end{align*} $$

Desired application code:

u = Vec2D(0,1)
v = Vec2D(1,0)
print(u + v)

a = u + v
w = Vec2D(1,1)
a == w

print(u - v)

print(u*v)

Class for vectors; implementation¶

 class Vec2D:
    def __init__(self, x, y):
        self.x = x;  self.y = y

    def __add__(self, other):
        return Vec2D(self.x+other.x, self.y+other.y)

    def __sub__(self, other):
        return Vec2D(self.x-other.x, self.y-other.y)

    def __mul__(self, other):
        return self.x*other.x + self.y*other.y

    def __abs__(self):
        return math.sqrt(self.x**2 + self.y**2)

    def __eq__(self, other):
        return self.x == other.x and self.y == other.y

    def __str__(self):
        return '(%g, %g)' % (self.x, self.y)

    def __ne__(self, other):
        return not self.__eq__(other)  # reuse __eq__

The repr special method: `eval(repr(p))` creates `p`¶

class MyClass:
    def __init__(self, a, b):
        self.a, self.b = a, b

    def __str__(self):
        """Return string with pretty print."""
        return 'a=%s, b=%s' % (self.a, self.b)

    def __repr__(self):
        """Return string such that eval(s) recreates self."""
        return 'MyClass(%s, %s)' % (self.a, self.b)

m = MyClass(1, 5)
print(m)      # calls m.__str__()

a=1, b=5

str(m)       # calls m.__str__()

s = repr(m)  # calls m.__repr__()
s

'MyClass(1, 5)'

m2 = eval(s) # same as m2 = MyClass(1, 5)
m2           # calls m.__repr__()

MyClass(1, 5)

m2 == m

False

Example with repr print method¶

class Person:
    """Stores a Person object which has a name."""

    def __init__(self, name):
        """Store name."""
        self.name = name

    def __str__(self):
        """Pretty print."""
        return str(self.name)

    def __repr__(self):
        """Print code for regenerating this instance. AND DEBUGGING!"""
        return 'Person(%r)' % self.name # %r brukes med __repr__

class Car:
    """Car that can take hold Person(s)"""
    
    def __init__(self, people):
        """Initiate car with people, either a person or list og persons."""
        if isinstance(people, list):
            self.people = people # people is already of list of Persons
        elif isinstance(people, Person):
            self.people = [people] # Make a list of the one Person given
            
    def __str__(self):
        """Pretty print"""
        # Join all the persons in self.people with ', '
        names = ', '.join(str(person) for person in self.people) 
        return 'The people in my car are %s.' % names
    
    def __repr__(self):
        """Print code for regenerating this instance. AND DEBUGGING!"""
        return 'Car(%r)' % self.people

andreas = Person('Andreas')
peter = Person('Peter')
car = Car([andreas, peter])
print(car)

The people in my car are Andreas, Peter.

print(repr(car)) # a lot easier to debug!

Car([Person('Andreas'), Person('Peter')])

new_car = eval(repr(car))
print(new_car)
print(new_car == car)

The people in my car are Andreas, Peter.
False

str(car) == str(new_car)

True

repr(car) == repr(new_car)

True

Class for complex numbers; functionality¶

Python already has a class complex for complex numbers, but implementing such a class is a good pedagogical example on class programming (especially with special methods).

Usage:

u = Complex(2,-1)
v = Complex(1)     # zero imaginary part
w = u + v
print(w)

(3, -1)

w != u

True

u*v

Complex(2, -1)

u < v

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-85-a9de6d3b6e83> in <module>()
----> 1 u < v

TypeError: '<' not supported between instances of 'Complex' and 'Complex'

print(w + 4)

(7, -1)

print(4 - w)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-87-322eebf29fa0> in <module>()
----> 1 print(4 - w)

TypeError: unsupported operand type(s) for -: 'int' and 'Complex'

Class for complex numbers; implementation (part 1)¶

class Complex:
    def __init__(self, real, imag=0.0):
        self.real = real
        self.imag = imag

    def __add__(self, other):
        return Complex(self.real + other.real,
                       self.imag + other.imag)

    def __sub__(self, other):
        return Complex(self.real - other.real,
                       self.imag - other.imag)

    def __mul__(self, other):
        return Complex(self.real*other.real - self.imag*other.imag,
                       self.imag*other.real + self.real*other.imag)

    def __truediv__(self, other): # __div__ for python 2.x
        ar, ai, br, bi = self.real, self.imag, \
                         other.real, other.imag # short forms
        r = float(br**2 + bi**2)
        return Complex((ar*br+ai*bi)/r, (ai*br-ar*bi)/r)

    def __abs__(self):
        return sqrt(self.real**2 + self.imag**2)

    def __neg__(self):   # defines -c (c is Complex)
        return Complex(-self.real, -self.imag)

    def __eq__(self, other):
        return self.real == other.real and \
               self.imag == other.imag

    def __ne__(self, other):
        return not self.__eq__(other)

    def __str__(self):
        return '(%g, %g)' % (self.real, self.imag)

    def __repr__(self):
        return 'Complex' + str(self)

    def __pow__(self, power):
        raise NotImplementedError(
          'self**power is not yet impl. for Complex')

Refining the special methods for arithmetics¶

Can we add a real number to a complex number?

u = Complex(1, 2)
w = u + 4.5 # actually works, since float object have .real and .imag attributes.
print(w)

(5.5, 2)

Problem: we have assumed that other is Complex. Remedy:

class Complex_robust(Complex):
    
    def __add__(self, other):
        if isinstance(other, (float,int)):
            other = Complex(other)
        return Complex(self.real + other.real,
                       self.imag + other.imag)

    # or
    '''
    def __add__(self, other):
        if isinstance(other, (float,int)):
            return Complex(self.real + other, self.imag)
        else:
            return Complex(self.real + other.real,
                           self.imag + other.imag)
    '''

Special methods for "right" operands; addition¶

What if we try this:

u = Complex_robust(1, 2)
w = 4.5 + u
print(w)

Problem: Python's float objects cannot add a Complex.

Remedy: if a class has an __radd__(self, other) special method, Python applies this for other + self

class Complex_more_robust(Complex_robust):
    
    def __radd__(self, other):
        """Rturn other + self."""
        # other + self = self + other:
        return self.__add__(other)

u = Complex_more_robust(1, 2)
w = 4.5 + u
'''
4.5.__add__(u) -> raise ValueError __add__ i float er ikke definert
for objekter av type Complex
u.__radd__(4.5)
'''
print(w)

(5.5, 2)

Special methods for "right" operands; subtraction¶

Right operands for subtraction is a bit more complicated since $a-b \neq b-a$:

class Complex_most_robust(Complex_more_robust):
    
    def __sub__(self, other):
        if isinstance(other, (float,int)):
            other = Complex(other)
        return Complex(self.real - other.real,
                       self.imag - other.imag)

    def __rsub__(self, other):
        if isinstance(other, (float,int)):
            other = Complex(other)
        return other.__sub__(self)

What's in a class?¶

 class A:
    """A class for demo purposes."""
    def __init__(self, value):
        self.v = value

Any instance holds its attributes in the self.__dict__ dictionary (Python automatically creates this dict)

a = A([1,2])
print(a.__dict__)  # all attributes

{'v': [1, 2]}

dir(a)            # what's in object a?

['__class__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 'v']

a.__doc__         # programmer's documentation of A

'A class for demo purposes.'

Ooops - we can add new attributes as we want!¶

a.myvar = 10            # add new attribute (!)
a.__dict__

{'myvar': 10, 'v': [1, 2]}

dir(a)

b = A(-1)
b.__dict__              # b has no myvar attribute

{'v': -1}

dir(b)

Summary of defining a class¶

Example on a defining a class with attributes and methods:

import matplotlib.pyplot as plt
import numpy as np

class Gravity:
    """Gravity force between two objects."""
    def __init__(self, m, M):
        self.m = m
        self.M = M
        self.G = 6.67428E-11 # gravity constant

    def force(self, r):
        G, m, M = self.G, self.m, self.M
        return G*m*M/r**2

    def visualize(self, r_start, r_stop, n=100):
        r = np.linspace(r_start, r_stop, n)
        g = self.force(r)
        title = 'm=%g, M=%g' % (self.m, self.M)
        plt.plot(r, g)
        plt.title(title)
        plt.show()

Summary of using a class¶

Example on using the class:

mass_moon = 7.35E+22
mass_earth = 5.97E+24

# make instance of class Gravity:
gravity = Gravity(mass_moon, mass_earth)

r = 3.85E+8  # earth-moon distance in meters
Fg = gravity.force(r)   # call class method
print('Fg =', Fg)
gravity.visualize(1, 24)

Summary of special methods¶

c = a + b implies c = a.__add__(b)
There are special methods for a+b, a-b, a*b, a/b, a**b, -a, if a:, len(a), str(a) (pretty print), repr(a) (recreate a with eval), etc.
With special methods we can create new mathematical objects like vectors, polynomials and complex numbers and write "mathematical code" (arithmetics)
The call special method is particularly handy: v = c(5) means v = c.__call__(5)
Functions with parameters should be represented by a class with the parameters as attributes and with a call special method for evaluating the function

Summarizing example: interval arithmetics for uncertainty quantification in formulas¶

Uncertainty quantification:

Consider measuring gravity $g$ by dropping a ball from $y=y_0$ to $y=0$ in time $T$:

$$ g = 2y_0T^{-2} $$

What if $y_0$ and $T$ are uncertain? Say $y_0\in [0.99,1.01]$ m and $T\in [0.43, 0.47]$ s. What is the uncertainty in $g$?

The uncertainty can be computed by interval arithmetics¶

Interval arithmetics.

Rules for computing with intervals, $p=[a,b]$ and $q=[c,d]$:

$p+q = [a + c, b + d]$
$p-q = [a - d, b - c]$
$pq = [\min(ac, ad, bc, bd), \max(ac, ad, bc, bd)]$
$p/q = [\min(a/c, a/d, b/c, b/d), \max(a/c, a/d, b/c, b/d)]$ ($[c,d]$ cannot contain zero)

Obvious idea: make a class for interval arithmetics!

Class for interval arithmetics¶

class IntervalMath:
    def __init__(self, lower, upper):
        self.lo = float(lower)
        self.up = float(upper)

    def __add__(self, other):
        a, b, c, d = self.lo, self.up, other.lo, other.up
        return IntervalMath(a + c, b + d)

    def __sub__(self, other):
        a, b, c, d = self.lo, self.up, other.lo, other.up
        return IntervalMath(a - d, b - c)

    def __mul__(self, other):
        a, b, c, d = self.lo, self.up, other.lo, other.up
        return IntervalMath(min(a*c, a*d, b*c, b*d),
                            max(a*c, a*d, b*c, b*d))

    def __truediv__(self, other): # From python 3.x
        a, b, c, d = self.lo, self.up, other.lo, other.up
        if c*d <= 0: return None
        return IntervalMath(min(a/c, a/d, b/c, b/d),
                            max(a/c, a/d, b/c, b/d))
    def __str__(self):
        return '[%g, %g]' % (self.lo, self.up)

Demo of the new class for interval arithmetics¶

Code:

I = IntervalMath   # abbreviate
a = I(-3,-2)
b = I(4,5)
a*b

expr = 'a+b', 'a-b', 'a*b', 'a/b'   # test expressions
for e in expr:
    print(e, '=', eval(e))

Shortcomings of the class¶

This code

a = I(4,5)
q = 2
b = a*q

leads to

      File "IntervalMath.py", line 15, in __mul__
        a, b, c, d = self.lo, self.up, other.lo, other.up
    AttributeError: 'float' object has no attribute 'lo'

Problem: IntervalMath times int is not defined.

Remedy: (cf. class Complex)

class IntervalArithmetics:
    ...
    def __mul__(self, other):
        if isinstance(other, (int, float)):      # NEW
            other = IntervalMath(other, other)   # NEW
        a, b, c, d = self.lo, self.up, other.lo, other.up
        return IntervalMath(min(a*c, a*d, b*c, b*d),
                            max(a*c, a*d, b*c, b*d))

(with similar adjustments of other special methods)

More shortcomings of the class¶

Try to compute g = 2*y0*T**(-2): multiplication of int (2) and IntervalMath (y0), and power operation T**(-2) are not defined

class IntervalArithmetics:
    ...
    def __rmul__(self, other):
        if isinstance(other, (int, float)):
            other = IntervalMath(other, other)
        return other*self

    def __pow__(self, exponent):
        if isinstance(exponent, int):
            p = 1
            if exponent > 0:
                for i in range(exponent):
                    p = p*self
            elif exponent < 0:
                for i in range(-exponent):
                    p = p*self
                p = 1/p
            else:   # exponent == 0
                p = IntervalMath(1, 1)
            return p
        else:
            raise TypeError('exponent must int')

Adding more functionality to the class: rounding¶

"Rounding" to the midpoint value:

a = IntervalMath(5,7)
float(a)

is achieved by

class IntervalArithmetics:
    ...
    def __float__(self):
        return 0.5*(self.lo + self.up)

Adding more functionality to the class: repr and str methods¶

class IntervalArithmetics:
    ...
    def __str__(self):
        return '[%g, %g]' % (self.lo, self.up)

    def __repr__(self):
        return '%s(%g, %g)' % \
          (self.__class__.__name__, self.lo, self.up)

Demonstrating the class: $g=2y_0T^{-2}$¶

g = 9.81
y_0 = I(0.99, 1.01)
Tm = 0.45                 # mean T
T = I(Tm*0.95, Tm*1.05)   # 10% uncertainty
print T

g = 2*y_0*T**(-2)
g

# computing with mean values:
T = float(T)
y = 1
g = 2*y_0*T**(-2)
print '%.2f' % g

Demonstrating the class: volume of a sphere¶

R = I(6*0.9, 6*1.1)   # 20 % error
V = (4./3)*pi*R**3
V

print V

print float(V)

# compute with mean values:
R = float(R)
V = (4./3)*pi*R**3
print V

20% uncertainty in $R$ gives almost 60% uncertainty in $V$