Wednesday, 24 April 2013

Beautify your math code. Part I


I've discovered few typical mistakes that people often do in matlab/gnu octave/python code. Maybe mistake is strong word, not mistakes, but suboptimal usage of environment capabilities (yes, I know mistake sounds better :)). Since I am familiar with python more that mathab, I'll provide few snippets written on python with numpy and scipy libraries, but idea remains the same in matlab also.

Grid functions.


We need setup geometry and some initial data:

import numpy as np
from scitools.numpyutils import meshgrid

h = 0.1
L = H = 1.
x = np.arange(0, L, h)
y = np.arange(0, H, h)
u_analytical = lambda x, y: 2*x + np.exp(y)
 
Ugly way to fill function values on a grid:
u1 = np.zeros( (len(x), len(y) ) )
for i, x_val in enumerate( x ):
    for j, y_val in enumerate( y ):
        u1[i][j] = u_analytical(x_val, y_val)
 
Good way to fill function values on a grid:
(xx, yy) = meshgrid(x, y, sparse=False, indexing='ij')
u1_ = u_analytical(xx, yy) 

Sparse matrices

First you should use sparse matrices, they rock! Second you should carefully choose sparse matrix format. As example why it should be done I'll quote documentation:

Each sparse format has certain advantages and disadvantages. For instance, adding new non-zero entries to a lil_matrix is fast, however changing the sparsity pattern of a csr_matrix requires a significant amount of work. On the other hand, operations such as matrix-vector multiplication and matrix-matrix arithmetic are much faster with csr_matrix than lil_matrix. A good strategy is to construct matrices using one format and then convert them to another that is better suited for efficient computation.
It is not very hard to choose correct sparse matrix format, especially in numerical computations, where in most cases you know all valuable information about your matrix before computation.

The saddest thing is that I saw bad written code in examples given by teachers in their numerical courses.