[Matrix-SIG] numpy interface

[Phil Austin]

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
(http://monet.uwaterloo.ca/blitz)
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
{
public:
  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);
      Bz.reference(tempArray);
    };  
  TiedArray(Py::seqref<Py::Object> pyob):Pn(pyob)
    {
      assert(Pn.rank()==2);
      assert(Pn.species()==PyArray_FLOAT);
      blitz::Array<float,2> tempArray((float*) Pn.to_C(),blitz::shape(Pn.dimension(1),
	      Pn.dimension(2)));
      Bz.reference(tempArray);
    };
};


static PyObject *
ex_strucout(PyObject* self, PyObject* args)
//
// calculate the 2nd order structure function given
// a 2-d NumPy array of floats
//
{
  Py_Initialize();
  import_array();
  PyImport_ImportModule("Numeric");

  Py::Tuple t(args);
  TiedArray<float,2> tau(t[0]);
  assert(tau.Pn.dimension(1)==tau.Pn.dimension(2));
  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);

  sf.Bz=0.;
  count.Bz=0.;

  for(int i=0;i<N;i++){
    for(int j=0;j<N;j++){
      if(!isnan(tau.Bz(i,j)))
	{
	  for(int k=0;k<N;k++){
	    for(int l=0;l<N;l++){
	      if(!isnan(tau.Bz(k,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.){
	sf.Bz(i,j)=sf.Bz(i,j)/count.Bz(i,j);
      }
    }
  }
  Py::Dict output;
  output["count"]=count.Pn;
  output["struct"]=sf.Pn;
  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

test=arange(4,typecode=Float32)
test.shape=(2,2)
print strucfun.strucout(test)