[Matrix-SIG] Re: blitz++ vs. NumPy?

[Phil Austin]

The following message is a courtesy copy of an article
that has been posted as well.



posting to comp.lang.python, henk13@my-dejanews.com writes:

> 
> I just wondered what the advantages and disadvantages would be of doing
> Python numerical computiation in a Blitz++ (http://seurat.uwaterloo.ca/blitz)
> framework rather than NumPy's Lapack/Blas framework (this is _not_ meant to
> offend the people at LLNL). The inherint OO-design of Blitz nicely extends
> the mental concept in which Python programs are developed. For instance, the
> Array class could be straightforwardly mapped to a Python Array extension
> type.
> 
> Is this some utopia or does it make sense?
> 
> Comments are welcome!
> 

Yes, the fit you describe is very close to completion thanks to Paul
Dubois,Todd Veldhuizen, Geoff Furnish and others.  Here's an example
of a Python extension I wrote in about 20 minutes last night, which
calculates the 2nd order structure function for a 256 x 256 pixel
satellite image of cloud liquid water path.  It uses Paul's CXX_Array.h
(see the source, http://xfiles.llnl.gov/CXX_Objects/cxx.htm and
particularly the last pages of
http://starship.skyport.net/~da/freebie). I've added a constructor
that creates an n-dimensional PyArray from dimensions specified by a
Blitz 1-d array.

Blitz arrays are reference-counted, and Todd's provided a constructor
that takes a raw C pointer (in this case to the PyArray data).
I can mirror the PyArrays with Blitz arrays of the
same dimensions, denoted by the "tied" prefix below and get the
template metaprograms, 2-dimensional subscripts etc., without paying
for a copy.  I return the structure function and the pixel counts to
Python in a dictionary.  The program strucfun.cxx below is compiled into
strucfun.so and then called from Python like this:

from Numeric import *
import strucfun

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


Some comments/questions:

1)  From my perspective the combination of Blitz and CXX_Objects pushes
NumPy across a threshold: extension writing has become something that
Fortran77 and Matlab programming scientists (like me) can do.  

2) Two things I wasn't able to figure out in 20 minutes are:

  a) How to add my Blitz-based Py::Array constructor to Paul's code
     using inheritence.  It looks like:

    explicit Array (blitz::Array<int,1> dimens, PyArray_TYPES t = PyArray_DOUBLE)
        : Sequence(FromAPI(PyArray_FromDims(dimens.size(), dimens.data(), t))) {
        validate();
    }

  b) Is there a general way to convert between a Blitz TinyVector of
     length N (where N isn't known at compile time)
     and a blitz::Array<int,1>?  I'd like to avoid specifying the
     dimensions to the Py::Array and the blitz::Array separately.

3) This worked the first time I tried it; everyone responsible for that
little miracle has my heartfelt thanks.

__________strucfun.cxx_____________

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

USING(namespace Py)
USING(namespace std)

blitz::Array<int,1> convertit( blitz::TinyVector<int,2> tinyDimens)
{
  blitz::Array<int,1> junk(2);
  for(int j=0;j < tinyDimens.length();++j){
    junk(j)=tinyDimens(j);
  }
  return junk;
}


static PyObject *
ex_strucout(PyObject* self, PyObject* args)
{
  Py_Initialize();
  import_array();
  PyImport_ImportModule("Numeric");

  Tuple t(args);
  Array tau(t[0]);
//
// change these to exceptions later
//
  assert(tau.rank()==2);
  assert(tau.dimension(1)==tau.dimension(2));
  assert(tau.species()==6);
  int N=tau.dimension(1);
  int Nbig=2*N-1;
  blitz::TinyVector<int,2> tinyDimens(Nbig,Nbig);
  Array sf(convertit(tinyDimens),PyArray_FLOAT);
  Array count(convertit(tinyDimens),PyArray_FLOAT);
  // make blitz arrays that point to the same data and have the
  // same shape
  blitz::Array<float,2> tiedTau((float*) sf.to_C(),tinyDimens);
  blitz::Array<float,2> tiedSf((float*) sf.to_C(),tinyDimens);
  blitz::Array<float,2> tiedCount((float*) count.to_C(),tinyDimens);
  tiedSf=0.;
  tiedCount=0.;

  for(int i=0;i<N;i++){
    for(int j=0;j<N;j++){
      if(!isnan(tiedTau(i,j)))
	{
	  for(int k=0;k<N;k++){
	    for(int l=0;l<N;l++){
	      if(!isnan(tiedTau(k,l)))
		{
		  tiedSf(k-i+N,l-j+N)=tiedSf(k-i+N,l-j+N) + pow((tiedTau(k,l) - tiedTau(i,j)),2);
		  tiedCount(k-i+N,l-j+N)=tiedCount(k-i+N,l-j+N) + 1.;
		}
	    }
	  }
	}
    }
  }
  for (int i=0; i<Nbig; i++) {
    for (int j=0; j<Nbig; j++) {
      if(tiedCount(i,j) > 0.){
	tiedSf(i,j)=tiedSf(i,j)/tiedCount(i,j);
      }
    }
  }
  Dict output;
  output["count"]=count;
  output["struct"]=sf;
  return new_reference_to(output);
}


extern "C" void initstrucfun();

static ExtensionModule* strucfun;

void initstrucfun()
{
    // experimental initialization stuff
    strucfun = new ExtensionModule("strucfun");
    strucfun->add("strucout", ex_strucout, "calculate the structure function");
}






-- 
------------------------------------------------------------
Phil Austin  (phil@geog.ubc.ca)
Department of Geography,                 Tel: (604) 822-2663
University of British Columbia, B.C.     Fax: (604) 822-6150