• C++_Extensions_that_use_NumPy_arrays

Writing C++ extensions that use NumPy arrays

Some basics

Passing numpy arrays to C++ using ctypes requires 3 steps:

1. In Python: make a ctypes object of the numpy array. For example:

import numpy as N
from numpyctypes import c_ndarray
c_myarray = c_ndarray(myarray, dtype = N.double, ndim = 3)

c_ndarray() checks if the given numpy array fulfills the requirements specified by the keyword parameters. The required numpyctypes.py module is attached to this wiki page. Use help(c_ndarray) to find out more.

2. In C++: the numpy arrays are converted to a C struct "numpyArray":

   1 #include <iostream>
   2 #include "ndarray.h"
   3 
   4 extern "C" {
   5  
   6 int myfunc(numpyArray<double> array)
   7 {
   8 // Do stuff
   9 }
  10 }

It would be nice if the numpy array could be converted right away to a C++ array class, but unfortunately this is not possible. ctypes only understands C variables, and a C struct is the closest thing to a C++ class. The C++ header file ndarray.h contains the definition of "numpyArray" and is attached to this wiki page.

3. In C++: convert the C struct "numpyArray" to a C++ "Ndarray" object.

   1 #include <iostream>
   2 #include "ndarray.h"
   3 
   4 extern "C" {
   5  
   6 int myfunc(numpyArray<double> array)
   7 {
   8     Ndarray<double,3> a(array);
   9    // Do stuff
  10 }
  11 }

Note that you need to know the number of dimensions (axes) of the array to define an "Ndarray" object. "Ndarray" is also defined in the C++ header file "ndarray.h". What's the gain? When using a multidimensional numpy array in a C/C++ extension, you can't simply use a[i][j[k] to access an element. Instead, you need to compute the location of the element using strides info, e.g. a[i*strides[0]+j*strides[1]+k*strides[2]]. The reason for this is the way numpy internally stores an array. Such a way to access an element is rather cumbersome and inconvenient. C++, however, allows to overload the [] operator to incorporate the strides information automatically. This is exactly what Ndarray does: a[i][j][k] is automatically converted to a[i*strides[0]+j*strides[1]+k*strides[2]] in a transparent way. The technique used is called template metaprogramming and it has the advantage that the conversion is done at compile time, not at runtime!

A worked out example

myextension.py

import numpy as N
from numpy.ctypeslib import load_library
from numpyctypes import c_ndarray

mylib = load_library('libmyextension', '.')       # '.' is the directory of the C++ lib  

def myfunc(array1, array2):
    arg1 = c_ndarray(array1, dtype=N.double, ndim = 3, shape = (4,3,2))
    arg2 = c_ndarray(array2, dtype=N.double, ndim = 3, shape = (4,3,2))
    return mylib.myfunc(arg1, arg2)

myextension.cpp

   1 #include <iostream>
   2 #include "ndarray.h"
   3 
   4 extern "C" {
   5  
   6 int myfunc(numpyArray<double> array1, numpyArray<double> array2)
   7 {
   8     Ndarray<double,3> a(array1);
   9     Ndarray<double,3> b(array2);
  10 
  11     double sum=0.0;
  12 
  13     for (int i = 0; i < a.getShape(0); i++)
  14     {
  15         for (int j = 0; j < a.getShape(1); j++)
  16         {
  17             for (int k = 0; k < a.getShape(2); k++)
  18             {
  19                 a[i][j][k] = 2.0 * b[i][j][k];
  20                 sum += a[i][j][k];
  21            }
  22         }
  23     }
  24     return sum;    
  25 }
  26 } // end extern "C"
  27 

myfunc() is extremely trivial, but it serves as an example. You need to compile myextension.cpp and make a shared library from it. The easiest way is to use Scons with the constructor file:

env = Environment()
env.Replace(CFLAGS=['-O2','-Wall','-ansi','-pedantic'])
env.SharedLibrary(target='myextension', source=['myextension.cpp'])

On a mac, this will create the myextension.dylib. Now you can use myfunc in a python program:

import numpy as N
from myextension import myfunc

a = N.zeros((4,3,2))
b = N.arange(4*3*2).reshape(4,3,2) * 1.0     # double array !
s = myfunc(a,b)
print s, a, b