Numerical Python Notes

Some random notes on using numarray and python for scientific computing. The numarray package allows for efficient use of arrays in python.

Passing arrays from python to C and preserving multidimensional indexing

The API for numarray allows for writing extension modules in C to process numarray array. When dealing with multidimensional arrays, it is convinent to be able to retain multidimensional indexing into the array. Since the data in a well-defined numarray array is contiguous, it is possible to do this indexing by defining a double ** pointer and allocating each row to point into the proper offset into the numarray data segment.

The sample code illustrates this. To build, use the setup.py as python setup.py build.

This defines the dump method to print out the contents of a single double-precision numarray array:

>>> import arraylib
>>> from numarray import *
>>> a = arange(64, type=Float64)
>>> a.shape = (16,4)
>>> arraylib.dump(a)
0.000000 1.000000 2.000000 3.000000
4.000000 5.000000 6.000000 7.000000
8.000000 9.000000 10.000000 11.000000
12.000000 13.000000 14.000000 15.000000
16.000000 17.000000 18.000000 19.000000
20.000000 21.000000 22.000000 23.000000
24.000000 25.000000 26.000000 27.000000
28.000000 29.000000 30.000000 31.000000
32.000000 33.000000 34.000000 35.000000
36.000000 37.000000 38.000000 39.000000
40.000000 41.000000 42.000000 43.000000
44.000000 45.000000 46.000000 47.000000
48.000000 49.000000 50.000000 51.000000
52.000000 53.000000 54.000000 55.000000
56.000000 57.000000 58.000000 59.000000
60.000000 61.000000 62.000000 63.000000
1
>>> print a
[[  0.   1.   2.   3.]
 [  4.   5.   6.   7.]
 [  8.   9.  10.  11.]
 [ 12.  13.  14.  15.]
 [ 16.  17.  18.  19.]
 [ 20.  21.  22.  23.]
 [ 24.  25.  26.  27.]
 [ 28.  29.  30.  31.]
 [ 32.  33.  34.  35.]
 [ 36.  37.  38.  39.]
 [ 40.  41.  42.  43.]
 [ 44.  45.  46.  47.]
 [ 48.  49.  50.  51.]
 [ 52.  53.  54.  55.]
 [ 56.  57.  58.  59.]
 [ 60.  61.  62.  63.]]

After the argument list is parsed in the extension module, we define a double ** pointer and allocate space to hold the row pointers.

  cArray = (double **) malloc(sizeof(double *)*nx);
  
  for (i = 0; i < nx; i++) {
    cArray[i] = (double *) NA_OFFSETDATA(pArray) + i*ny;
  }
  

where here, pArray is the PyArrayObject pointing to the passed numarray array (see the sample code for full details). Now we can index the array as cArray[i][j] in C.