[MATRIX-SIG] [PYTHON MATRIX-SIG] SWIG, NumPy and C++ templated arrays

[Phil Austin]

Someone on the SWIG mailing list asked about SWIG and numeric Python.
I've put together an example of our initial attempt using MV++
arrays (see ftp://ftp.geog.ubc.ca/pub/swig/example.tar.gz (11 Kbytes)).
I'm including below the README and interface files from the example.
Any comments/questions appreciated.



	 SWIG-python example for templated C++ arrays and the
		     numerical python extension.

	     Michael Cherkassoff (mcherk@geog.ubc.ca) and
		   Phil Austin (phil@geog.ubc.ca).


This directory contains example code which glues Python NumPy arrays
to C++ templated arrays (all necessary headers from Roldan Poza's MV++
library are included).  This approach currently works for us, and we
will soon be using it wholesale as we move to either MV++'s replacement
(http://math.nist.gov/tnt/lib/) or one of the new expression-template
based array classes (like Blitz++
http://monet.uwaterloo.ca/blitz/oon.html, or A++
http://www.c3.lanl.gov/~dquinlan/A++P++.html ).

Below I describe how the example works, but for the savants, we'd
appreciate help with three unresolved issues with the interface file
(plus, of course, any other comments on either the python or c++).

1) Because SWIG doesn't handle templates, we fool it with a second
set of typedefs (see cutmod.i and cutmodfake.h)--is there a more
elegant approach?

2) We were unable to 'do the right thing' in the typemap to 
go from $source to a PyArrayObject*. That is, in cutmod.i (see
comments therein):

aa=(PyArrayObject *)PyArray_ContiguousFromObject($source, PyArray_Float,1,1);

gives a zero pointer, and so we resorted to:

aa = (PyArrayObject *) $source;

3) We're still confused about reference counting.  Should we be
calling PyDECREF(aa) in cutmod.i?



The code:

 Cutmod.cpp defines a C++ function "cut", which cuts data into a
 series of histogram bins.  Cutmod.i gives the SWIG interface, and
 test.py contains a simple test.  The cut function (modeled on cut in
 Splus) takes either two python lists or two arrays (data and bins),
 and returns a list of lists containing the indexes of the data
 elements that fit into the specified bin intervals.

i.e. in Python, typing 

_______________
from Numeric import *
import cutmod


data=array( [10, 50, 24, 45, 70, 59, 21, 33, 44, 65,
	     13, 11, 99, 15, 14, 78, 93, 45, 18, 32,
	     37, 16, 18, 84, 75, 48, 33, 17, 61, 90])


bins =array( [10, 20, 30, 40, 50, 60, 70, 80, 90, 100])

output=cutmod.cut(data,bins)

_________________

produces:

output
[[0, 11, 10, 14, 13, 21, 27, 18, 22], [6, 2], [19, 7, 26, 20], 
[8, 3, 17, 25], [1, 5], [28, 9], [4, 24, 15], [23], [29, 16, 12]]


to install, just untar the file and make with

make cutmode.so

invoke python and

import test

the output should be test.output

tested only on Linux 2.0.27 using gcc 2.7.2.1


Appendix:

For those unfamiliar with MV++:



Here's the way I would implement cut in Python:


def cut(data,bins): #input: data (numarray), bins (numarray) to cut data. 
                    #output: binlist (list of lists), indexes in each bin
        binlength=bins.shape[0]
	datalength=data.shape[0]
#populate an empty list of lists
	binlist=[0]*(binlength-1)
	for i in range(len(binlist)):
	    binlist[i]=[]
	index=range(bins.shape[0])
	sortedIndex=argsort(data)
	i=0
#move through the bins, dropping the sorted indices in the appropriate bin
	for leftbinmarker in range(binlength-1):
	    if i > (datalength-1): return binlist
	    while data[sortedIndex[i]] < bins[leftbinmarker+1]:
		binlist[leftbinmarker].append(sortedIndex[i])
		i=i+1
                if i > (datalength-1): return binlist

cutmod.cc is a straight translation of this code using
MV++ vectors instead of lists, with a little extra error checking.


==> cutmod.i <==
// Code for SWIG to generate wrapper for C++ function that
// is to be called from Python.

// On the Python's end the function takes two arguments, each of
// them either 1-D NumPy array or a Python list of floats. It returns a
// list of lists of integers.

// On the C++ end the function takes two MV++ vectors of floats.
// It returns a MV++ vector of vectors of integers.

/* Here are the contents of cutmod.h:

#include "mvmtp.h"
typedef MV_Vector<int> VectorInt;
typedef MV_Vector<float> VectorFloat;
typedef MV_Vector<VectorInt> VectorList;

******************************************************/

%module cutmod

%{
#include "cutmod.h"       // this is going to be included without
#include "arrayobject.h"  // changes into wrapper code 
%}

%include "cutmodfake.h"   // this is for SWIG to look at.
			  // cutmod.h cannot be put here since
			  // SWIG doesn't like templates. And without
			  // something here it generates
			  // PyArg_ParseTuple(args,"ss",&_obj0,&_obj1)
			  // instead of desired
			  // PyArg_ParseTuple(args,"OO",&_obj0,&_obj1)

%typemap(python,in) VectorFloat {
	int i, size;
	enum TYPE {List, Array};
	TYPE type;
	PyArrayObject * aa;
	float * bb;

	if (PyList_Check($source))
	   {
	   size = PyList_Size($source);
	   type = List;
	   }
	else if (PyArray_Check($source))
	   {
	   size = PyArray_Size($source);
	   type = Array;
	   aa = (PyArrayObject *) $source;    // This is a hack, but it
	   bb = (float *) aa->data;           // works. The proper way:
			// with PyArray_ContiguousFromObject() does not. 8-(

/* The proper code would have been:

aa=(PyArrayObject *)PyArray_ContiguousFromObject($source, PyArray_Float,1,1);

but this gives zero pointer for some reason obscure for me  */
	   }
	else
	   {
	   PyErr_SetString(PyExc_TypeError, "neither List nor Array");
	   return NULL;
	   }

	VectorFloat a(size);     // Cannot declare $target since it is
				 // already being done by SWIG.
	for (i=0; i<size; i++) 
	   {
	   switch(type) 
	      {
	      case List:
	         a(i) = PyFloat_AsDouble(PyList_GetItem($source,i));
		 break;
	      case Array:
		 a(i) = bb[i];
		 break;
	      }
	   }
	$target = a;

// I'm not quite sure but there may be needed something like
// PyDECREF(aa); here. So far things haven't crushed either way -
// with or without it

}

%typemap(python,out) VectorList {
	int size = $source.size();
	$target = PyList_New(size);
	for (int i=0; i<size; i++)
		{
		int sizei = $source(i).size();
		PyObject * listi = PyList_New(sizei);
		for(int j=0; j<sizei; j++)
			PyList_SetItem(listi,j,PyInt_FromLong($source(i)(j)));
		PyList_SetItem($target,i,listi);
		}
// Same thing about DECREFing listi inside or outside the cycle.
}

extern VectorList cut(VectorFloat,VectorFloat);


Phil Austin		INTERNET: phil@geog.ubc.ca
(604) 822-2175		FAX:	  (604) 822-6150

http://www.geog.ubc.ca/~phil
Associate Professor
Atmospheric Sciences Programme
Geography #217
University of British Columbia
1984 W Mall
Vancouver, BC  V6T 1Z2
CANADA


_______________
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
_______________

_______________
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
_______________