This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.

Summary

Chapter topics

Classes

A class contains attributes, which are variables (data attributes) and functions (method attributes, also called just methods). A first rough overview of a class can be to just list the attributes, e.g., in a UML diagram.

Below is a sample class with three data attributes (m, M, and G) and three methods (a constructor, force, and visualize). The class represents the gravity force between two masses. This force is computed by the force method, while the visualize method plots the force as a function of the distance between the masses.

class Gravity(object):
    """Gravity force between two physical objects."""

    def __init__(self, m, M):
        self.m = m           # mass of object 1
        self.M = M           # mass of object 2
        self.G = 6.67428E-11 # gravity constant, m**3/kg/s**2

    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):
        from scitools.std import plot, linspace
        r = linspace(r_start, r_stop, n)
        g = self.force(r)
        title='Gravity force: m=%g, M=%g' % (self.m, self.M)
        plot(r, g, title=title)

Note that to access attributes inside the force method, and to call the force method inside the visualize method, we must prefix with self. Also recall that all methods must take self, "this" instance, as first argument, but the argument is left out in calls. The assignment of a data attributes to a local variable (e.g., G = self.G) inside methods is not necessary, but here it makes the mathematical formula easier to read and compare with standard mathematical notation.

This class (found in file Gravity.py) can be used to find the gravity force between the Moon and the Earth:

mass_moon = 7.35E+22;  mass_earth = 5.97E+24
gravity = Gravity(mass_moon, mass_earth)
r = 3.85E+8  # Earth-Moon distance in meters
Fg = gravity.force(r)
print 'force:', Fg

Special methods

A collection of special methods, with two leading and trailing underscores in the method names, offers special syntax in Python programs.

The table below provides an overview of the most important special methods.

Terminology

The important computer science topics in this document are

• classes
• attributes
• methods
• constructor (__init__)
• special methods (__add__, __str__, __ne__, etc.)

Example: interval arithmetic

Input data to mathematical formulas are often subject to uncertainty, usually because physical measurements of many quantities involve measurement errors, or because it is difficult to measure a parameter and one is forced to make a qualified guess of the value instead. In such cases it could be more natural to specify an input parameter by an interval \( [a,b] \), which is guaranteed to contain the true value of the parameter. The size of the interval expresses the uncertainty in this parameter. Suppose all input parameters are specified as intervals, what will be the interval, i.e., the uncertainty, of the output data from the formula? This section develops a tool for computing this output uncertainty in the cases where the overall computation consists of the standard arithmetic operations.

To be specific, consider measuring the acceleration of gravity by dropping a ball and recording the time it takes to reach the ground. Let the ground correspond to \( y=0 \) and let the ball be dropped from \( y=y_0 \). The position of the ball, \( y(t) \), is then $$ \begin{equation*} y(t) = y_0 - \frac{1}{2}gt^2 \thinspace . \end{equation*} $$ If \( T \) is the time it takes to reach the ground, we have that \( y(T)=0 \), which gives the equation \( \frac{1}{2}gT^2=y_0 \), with solution $$ \begin{equation*}g = 2y_0T^{-2} \thinspace . \end{equation*} $$ In such experiments we always introduce some measurement error in the start position \( y_0 \) and in the time taking (\( T \)). Suppose \( y_0 \) is known to lie in \( [0.99,1.01] \) m and \( T \) in \( [0.43, 0.47] \) s, reflecting a 2% measurement error in position and a 10% error from using a stop watch. What is the error in \( g \)? With the tool to be developed below, we can find that there is a 22% error in \( g \).

Problem

Assume that two numbers \( p \) and \( q \) are guaranteed to lie inside intervals, $$ \begin{equation*} p=[a,b], \quad q = [c, d] \thinspace . \end{equation*} $$ The sum \( p+q \) is then guaranteed to lie inside an interval \( [s,t] \) where \( s=a+c \) and \( t=b+d \). Below we list the rules of interval arithmetic, i.e., the rules for addition, subtraction, multiplication, and division of two intervals:

• \( 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)] \) provided that \( [c,d] \) does not contain zero

• The acceleration of gravity, \( g=2y_0T^{-2} \), given a 2% uncertainty in \( y_0 \): \( y_0=[0.99,1.01] \), and a 10% uncertainty in \( T \): \( T=[T_m\cdot 0.95, T_m\cdot 1.05] \), with \( T_m=0.45 \).
• The volume of a sphere, \( V=\frac{4}{3}\pi R^3 \), given a 20% uncertainty in \( R \): \( R=[R_m\cdot 0.9, R_m\cdot 1.1] \), with \( R_m =6 \).

Solution

The new type is naturally realized as a class IntervalMath whose data consist of the lower and upper bound of the interval. Special methods are used to implement arithmetic operations and printing of the object. Having understood class Vec2D from the section Example: Class for vectors in the plane, it should be straightforward to understand the class below:

class IntervalMath(object):
    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 __div__(self, other):
        a, b, c, d = self.lo, self.up, other.lo, other.up
        # [c,d] cannot contain zero:
        if c*d <= 0:
            raise ValueError\ 
                  ('Interval %s cannot be denominator because '\ 
                  'it contains zero' % other)
        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)

The code of this class is found in the file IntervalMath.py. A quick demo of the class can go as

I = IntervalMath
a = I(-3,-2)
b = I(4,5)
expr = 'a+b', 'a-b', 'a*b', 'a/b'
for e in expr:
    print '%s =' % e, eval(e)

The output becomes

a+b = [1, 3]
a-b = [-8, -6]
a*b = [-15, -8]
a/b = [-0.75, -0.4]

This gives the impression that with very short code we can provide a new type that enables computations with interval arithmetic and thereby with uncertain quantities. However, the class above has severe limitations as shown next.

Consider computing the uncertainty of \( aq \) if \( a \) is expressed as an interval \( [4,5] \) and \( q \) is a number (float):

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

This does not work so well:

  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'

The problem is that a*q is a multiplication between an IntervalMath object a and a float object q. The __mul__ method in class IntervalMath is invoked, but the code there tries to extract the lo attribute of q, which does not exist since q is a float.

We can extend the __mul__ method and the other methods for arithmetic operations to allow for a number as operand - we just convert the number to an interval with the same lower and upper bounds:

    def __mul__(self, other):
        if isinstance(other, (int, float)):
            other = IntervalMath(other, 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))

Looking at the formula \( g=2y_0T^{-2} \), we run into a related problem: now we want to multiply \( 2 \) (int) with \( y_0 \), and if \( y_0 \) is an interval, this multiplication is not defined among int objects. To handle this case, we need to implement an __rmul__(self, other) method for doing other*self, as explained in the section Special methods for "right" operands:

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

Similar methods for addition, subtraction, and division must also be included in the class.

Returning to \( g=2y_0T^{-2} \), we also have a problem with \( T^{-2} \) when \( T \) is an interval. The expression T**(-2) invokes the power operator (at least if we do not rewrite the expression as 1/(T*T)), which requires a __pow__ method in class IntervalMath. We limit the possibility to have integer powers, since this is easy to compute by repeated multiplications:

    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')

Another natural extension of the class is the possibility to convert an interval to a number by choosing the midpoint of the interval:

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

float(a) calls a.__float__(), which we implement as

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

A __repr__ method returning the right syntax for recreating the present instance is also natural to include in any class:

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

We are now in a position to test out the extended class IntervalMath.

>>> g = 9.81
>>> y_0 = I(0.99, 1.01)      # 2% uncertainty
>>> Tm = 0.45                # mean T
>>> T = I(Tm*0.95, Tm*1.05)  # 10% uncertainty
>>> print T
[0.4275, 0.4725]
>>> g = 2*y_0*T**(-2)
>>> g
IntervalMath(8.86873, 11.053)
>>> # Compute with mean values
>>> T = float(T)
>>> y = 1
>>> g = 2*y_0*T**(-2)
>>> print '%.2f' % g
9.88

Another formula, the volume \( V=\frac{4}{3}\pi R^3 \) of a sphere, shows great sensitivity to uncertainties in \( R \):

>>> Rm = 6
>>> R = I(Rm*0.9, Rm*1.1)   # 20 % error
>>> V = (4./3)*pi*R**3
>>> V
IntervalMath(659.584, 1204.26)
>>> print V
[659.584, 1204.26]
>>> print float(V)
931.922044761
>>> # Compute with mean values
>>> R = float(R)
>>> V = (4./3)*pi*R**3
>>> print V
904.778684234

Here, a 20% uncertainty in \( R \) gives almost 60% uncertainty in \( V \), and the mean of the \( V \) interval is significantly different from computing the volume with the mean of \( R \).

The complete code of class IntervalMath is found in IntervalMath.py. Compared to the implementations shown above, the real implementation in the file employs some ingenious constructions and help methods to save typing and repeating code in the special methods for arithmetic operations. You can read more about interval arithmetics on Wikipedia.

• ← Prev
• Next →