[Matrix-SIG] numpy interface

Phil Austin phil@geog.ubc.ca
Tue, 20 Apr 1999 10:48:54 -0700

Christian Tismer writes:
 > I don't believe that you need SWIG to call NumPy from C, since
 > NumPy's C API can be used like Python's C API.

Right -- and you might also want to have a look at Paul Dubois'
CXX_Array.h (http://xfiles.llnl.gov/CXX_Objects/cxx.htm), which
provides a very clean C++ interface to PyArrays (albeit with 1-d
subscripting only), dicts, tupples, etc.  

From the C++ side, your code would look something
like this:

Py::Array count(pyob);   //construct a 1-d array Pn with pyob's data
Py::Dict output;  //construct a C++ version of a Python dictionary
output["count"]=count; //put the count array in the output dictionary
return Py::new_reference_to(output);  //send the dictionary with the 
                                      // count array back to Python

To get around the 1-d subscripting limitation I use a class I
call TiedArray, which creates a Py::Array and a Blitz++ array
that share the same data.  I can then use Blitz notation on the
C++ side, while easily shipping the data back to Python.  

Here's an example I posted about 4 months ago which calculates
something called the 2 dimensional isotropic 
2nd order structure function for a
satellite image.  I've modified the Py::Array constructor
to support Blitz TinyVectors -- I'm happy to send you this
two-line modification if you're a Blitz user.

#include "Python.h"
#include "CXX_Objects.h"
#include "CXX_Extensions.h"
#include <blitz/array.h>
#include "CXX_Array.h"
#include <math.h>

template<typename T,int Ndim> class TiedArray;

template<> class TiedArray<float,2> 
//take an existing PyObject or a N-length TinyVector and
//construct an N-dimensional array shared by Python and Blitz.
//Note that the Python array controls the memory deallocation;
//the Blitz array just holds a reference to the Pn's data
  Py::Array Pn;
  blitz::Array<float,2> Bz;
  TiedArray(blitz::TinyVector<int,2> dimens):Pn(dimens,PyArray_FLOAT) //added to CXX_Array by PA: 98/12/30
      blitz::Array<float,2> tempArray((float*) Pn.to_C(),dimens);
  TiedArray(Py::seqref<Py::Object> pyob):Pn(pyob)
      blitz::Array<float,2> tempArray((float*) Pn.to_C(),blitz::shape(Pn.dimension(1),

static PyObject *
ex_strucout(PyObject* self, PyObject* args)
// calculate the 2nd order structure function given
// a 2-d NumPy array of floats

  Py::Tuple t(args);
  TiedArray<float,2> tau(t[0]);
  int N=tau.Pn.dimension(1);
  int Nbig=2*N-1;
  blitz::TinyVector<int,2> tinyDimens(Nbig,Nbig);
  TiedArray<float,2> sf(tinyDimens);
  TiedArray<float,2> count(tinyDimens);


  for(int i=0;i<N;i++){
    for(int j=0;j<N;j++){
	  for(int k=0;k<N;k++){
	    for(int l=0;l<N;l++){
		  sf.Bz(k-i+N,l-j+N)=sf.Bz(k-i+N,l-j+N) + pow((tau.Bz(k,l) - tau.Bz(i,j)),2);
		  count.Bz(k-i+N,l-j+N)=count.Bz(k-i+N,l-j+N) + 1.;
  for (int i=0; i<Nbig; i++) {
    for (int j=0; j<Nbig; j++) {
      if(count.Bz(i,j) > 0.){
  Py::Dict output;
  return Py::new_reference_to(output);

extern "C" void initstrucfun();

static Py::ExtensionModule* strucfun;

void initstrucfun()
    // experimental initialization stuff
    strucfun = new Py::ExtensionModule("strucfun");
    strucfun->add("strucout", ex_strucout, "calculate the structure function");
    Py::Dict d = strucfun->initialize();

___________compile above into strucfun.so and use this from python

from Numeric import *
import strucfun

print strucfun.strucout(test)