$$ \newcommand{\tp}{\thinspace .} $$

 

 

 

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

Example: Class for vectors in the plane

This section explains how to implement two-dimensional vectors in Python such that these vectors act as objects we can add, subtract, form inner products with, and do other mathematical operations on. To understand the forthcoming material, it is necessary to have digested the section Special methods, in particular the sections Adding objects and Arithmetic operations and other special methods.

Some mathematical operations on vectors

Vectors in the plane are described by a pair of real numbers, \( (a, b) \). There are mathematical rules for adding and subtracting vectors, multiplying two vectors (the inner or dot or scalar product), the length of a vector, and multiplication by a scalar: $$ \begin{align} (a,b) + (c,d) &= (a+c, b+d), \tag{4}\\ (a,b) - (c,d) &= (a-c, b-d), \tag{5}\\ (a,b)\cdot(c,d) &= ac + bd, \tag{6}\\ ||(a,b)|| &= \sqrt{(a,b)\cdot(a,b)} \thinspace . \tag{7} \end{align} $$ Moreover, two vectors \( (a,b) \) and \( (c,d) \) are equal if \( a=c \) and \( b=d \).

Implementation

We may create a class for plane vectors where the above mathematical operations are implemented by special methods. The class must contain two data attributes, one for each component of the vector, called x and y below. We include special methods for addition, subtraction, the scalar product (multiplication), the absolute value (length), comparison of two vectors (== and !=), as well as a method for printing out a vector.

class Vec2D(object):
    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 __add__, __sub__, __mul__, __abs__, and __eq__ methods should be quite straightforward to understand from the previous mathematical definitions of these operations. The last method deserves a comment: here we simply reuse the equality operator __eq__, but precede it with a not. We could also have implemented this method as

    def __ne__(self, other):
        return self.x != other.x or self.y != other.y

Nevertheless, this implementation requires us to write more, and it has the danger of introducing an error in the logics of the boolean expressions. A more reliable approach, when we know that the __eq__ method works, is to reuse this method and observe that a != b means not (a == b).

A word of warning is in place regarding our implementation of the equality operator (== via __eq__). We test for equality of each component, which is correct from a mathematical point of view. However, each vector component is a floating-point number that may be subject to rounding errors both in the representation on the computer and from previous (inexact) floating-point calculations. Two mathematically equal components may be different in their inexact representations on the computer. The remedy for this problem is to avoid testing for equality, but instead check that the difference between the components is sufficiently small. The function numpy.allclose can be used for this purpose:

if a == b:

by

if numpy.allclose(a, b):

A more reliable equality operator can now be implemented:

class Vec2D(object):
    ...
    def __eq__(self, other):
        return numpy.allclose(self.x, other.x) and \ 
               numpy.allclose(self.y, other.y)

As a rule of thumb, you should never apply the == test to two float objects.

The special method __len__ could be introduced as a synonym for __abs__, i.e., for a Vec2D instance named v, len(v) is the same as abs(v), because the absolute value of a vector is mathematically the same as the length of the vector. However, if we implement

    def __len__(self):
        # Reuse implementation of __abs__
        return abs(self)  # equiv. to self.__abs__()

we will run into trouble when we compute len(v) and the answer is (as usual) a float. Python will then complain and tell us that len(v) must return an int. Therefore, __len__ cannot be used as a synonym for the length of the vector in our application. On the other hand, we could let len(v) mean the number of components in the vector:

    def __len__(self):
        return 2

This is not a very useful function, though, as we already know that all our Vec2D vectors have just two components. For generalizations of the class to vectors with \( n \) components, the __len__ method is of course useful.

Usage

Let us play with some Vec2D objects:

>>> u = Vec2D(0,1)
>>> v = Vec2D(1,0)
>>> w = Vec2D(1,1)
>>> a = u + v
>>> print a
(1, 1)
>>> a == w
True
>>> a = u - v
>>> print a
(-1, 1)
>>> a = u*v
>>> print a
0
>>> print abs(u)
1.0
>>> u == v
False
>>> u != v
True

When you read through this interactive session, you should check that the calculation is mathematically correct, that the resulting object type of a calculation is correct, and how each calculation is performed in the program. The latter topic is investigated by following the program flow through the class methods. As an example, let us consider the expression u != v. This is a boolean expression that is True since u and v are different vectors. The resulting object type should be bool, with values True or False. This is confirmed by the output in the interactive session above. The Python calculation of u != v leads to a call to

u.__ne__(v)

which leads to a call to

u.__eq__(v)

The result of this last call is False, because the special method will evaluate the boolean expression

0 == 1 and 1 == 0

which is obviously False. When going back to the __ne__ method, we end up with a return of not False, which evaluates to True.

Comment

For real computations with vectors in the plane, you would probably just use a Numerical Python array of length 2. However, one thing such objects cannot do is evaluating u*v as a scalar product. The multiplication operator for Numerical Python arrays is not defined as a scalar product (it is rather defined as \( (a,b)\cdot(c,d) = (ac, bd) \)). Another difference between our Vec2D class and Numerical Python arrays is the abs function, which computes the length of the vector in class Vec2D, while it does something completely different with Numerical Python arrays.