Re: BZSUPPORT: Blitz++ and Python

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 454-5900

Harry Zuzan writes:
 > Hi,
 > 
 > I have been using Blitz for while now and at the same time I
 > have been using Python for prototyping some of my numeric work.
 > 
 > What I want to do now is write C++ extensions to Python, incorporating
 > Blitz code for the heavier numeric work.  This involves adding to the
 > to Python's numeric library.  I can do it in C easily enough but I would
 > rather work in C++ and use Blitz when it suits me.
 > 
 > Q.  Has anybody on this list been extending Python with Blitz?
 > 
 > Harry Zuzan
 > ha...@st...

Yes, this works well using http://CXX.sourceforge.net/.  I've
appended a simple example (using MTL, but it works just the
same with blitz, an mtl::external_vec is just MTL's way
of instantiating with NeverDeleteData.  Some notes:

1) get the CXX_5.0 tagged version.

cvs -d:pserver:ano...@cv...:/cvsroot/CXX login
cvs -z3 -d:pserver:ano...@cv...:/cvsroot/CXX co -r CXX_5_0 CXX

2) if want two-dimensional numpy arrays, you'll need my patch to
   CXX_Array (not submitted yet, so mail me).

Here's rbins.hh, rbins.cc, testrbins.py and the output

_______rbins.hh______

#ifndef __rbins__h
#define __rbins__h

//----------------------------------*-C++-*----------------------------------//
// Copyright 1998 The Regents of the University of California. 
// All rights reserved. See LEGAL.LLNL for full text and disclaimer.
//---------------------------------------------------------------------------//

#include "Python.h"
#include "CXX_Objects.h"
#include "CXX_Extensions.h"

static std::string rbins_h_ID=("$Id: rbins.hh,v 1.3 2000/07/17 20:52:45 phil Exp phil $");

extern "C" void initrbins();

class rbins_module : public Py::ExtensionModule<rbins_module>
{
public:
  rbins_module();
  virtual ~rbins_module() {}
private:
  Py::Object ex_rbins (const Py::Tuple &args, const Py::Dict &kws);
};

#endif

______________rbins.cc____________

#include "Python.h"
#include "CXX_Array.h"
#include "CXX_Objects.h"
#include "Numeric/arrayobject.h"
#include "rbins.hh"
#include <cmath>
#include <mtl/mtl.h>
#include <mtl/linalg_vec.h>
#include <mtl/matrix.h>

using namespace Py;
using namespace std;

//
// pair the memory of each Py::Array with a
// corresponding mtl object, and manipulate the
// Py::Array data using mtl
//

typedef mtl::external_vec<float> MTLvectorf; 
typedef mtl::external_vec<int> MTLvectori; 

typedef mtl::matrix< int,
  mtl::rectangle<>, 
  mtl::dense<mtl::external>,
  mtl::row_major>::type MTLmatrixi; 

void initrbins()
{
  static rbins_module *rbins = new rbins_module;
}

static string rbins_cc_ID=("$Id: rbins.cc,v 1.7 2000/07/15 19:16:46 phil Exp phil $");

rbins_module::rbins_module(): ExtensionModule<rbins_module>( "rbins" )
{
  import_array();
  add_keyword_method("rbins", &rbins_module::ex_rbins, 
		     "rbins(numrows,numcols,numbins=numbins | binsize=binsize) -> dictionary\n");
  static string description;
  description += "\nrbins takes an image size (numrows,numcols)\n";
  description += "and either the numbins or the binsize and \n";
  description += "returns a dictionary with the bin centers and edges and \n";
  description += "an integer (numrows,numcols) image containing the bin number\n";
  description += "of every pixel\n";

  initialize( description.c_str() );

  Dict d( moduleDictionary() );
  d["rcsid"]=String(rbins_cc_ID);
}

Object rbins_module::ex_rbins(const Tuple &args, const Dict &kws)
{

//verify that there are exactly two argments (numrows and numcols)
//will throw IndexError if not

  args.verify_length(2);

//verify that each can be safely converted to int
//(will throw if not)

  int numrows=Int(args[0]);
  int numcols=Int(args[1]);

  int numbins;
  float binsize;
  int minrad=0;
  int maxrad=(numrows-1)*(numrows-1) + (numcols-1)*(numcols-1);

//now go through keywords dictionary, looking for either numbins or
//binsize

  bool foundkey = false;
  List keylist = kws.keys();
  for(List::iterator i = keylist.begin(); i != keylist.end(); ++i){
	{
	  if(foundkey){
		throw KeyError("need either numbins or binsize specified (not both)");
	  }
	  string keyword=String(*i);
	  if (keyword == "numbins"){
		numbins=Int(kws["numbins"]);
		binsize = (maxrad-minrad)/((float)numbins);
		foundkey = true;
	  }
	  else if(keyword == "binsize"){
		binsize = Float(kws["binsize"]);
		numbins = ceil((maxrad - minrad)/binsize);
		foundkey=true;
	  }
	  else{
		throw KeyError("only  numbins or binsize keywords accepted");
	  }
	}
  }
  if(!foundkey){
	cout << "selecting a default bin number of 10" << endl;
	numbins=10;
	binsize = (maxrad-minrad)/((float)numbins);
  }

  //
  // let Python manage the memory, but use paired mtl arrays to manipulate it
  // (note the mtl typedef specifies external memory)
  //

  Array binedges(numbins+1,PyArray_FLOAT);
  MTLvectorf mtledges((float*) binedges.to_C(),numbins+1);
  Array bincenters(numbins,PyArray_FLOAT);
  MTLvectorf mtlcenters((float*) bincenters.to_C(),numbins);

  mtledges[0] = minrad;
  for(int i=1;i <= numbins; ++i){
	mtledges[i] = minrad + i*binsize;
	mtlcenters[i-1]=mtledges[i-1] + (mtledges[i] - mtledges[i-1])/2.;
  }

  //
  // create a 2-d Py::Array to hold the indexes
  //

  vector<int> dimens;
  dimens.push_back(numrows); dimens.push_back(numcols);
  Array binarray(dimens,PyArray_INT);
  MTLmatrixi mtlbinarray((int*) binarray.to_C(),numrows,numcols);

  //
  // calculate the bin index for every image index 
  //

  int rowsq,index,distance;
  for(int row = 0; row < numrows; ++row){
	rowsq=row*row;
	for(int col=0; col < numcols; ++col){
	  distance = rowsq + col*col;
	  index = floor((distance - minrad)/binsize);
	  mtlbinarray(row,col)=index;
	}
  }

  Dict out;
  out["centers"]=bincenters;
  out["edges"]=binedges;
  out["bins"]=binarray;
  return out;
}

____________testrbins.py___________

import rbins,sys

print rbins.__doc__

print rbins.rbins.__doc__

print "\n",rbins.rcsid

numrows,numcols=(5,5)

print "\nindexes:"

for row in range(numrows):
	print " "
	for col in range(numcols):
		sys.stdout.write("%3d " % (row**2 + col**2))

out=rbins.rbins(numrows,numcols,numbins=4)

print "\n\nbin edges: ",out['edges'],"\n"
print "bin centers: ",out['centers'],"\n"

bins=out['bins']
print "bins:"
for row in range(bins.shape[0]):
	print bins[row,:]

_____________output___________

rbins takes an image size (numrows,numcols)
and either the numbins or the binsize and 
returns a dictionary with the bin centers and edges and 
an integer (numrows,numcols) image containing the bin number
of every pixel

rbins(numrows,numcols,numbins=numbins | binsize=binsize) -> dictionary

$Id: rbins.cc,v 1.7 2000/07/15 19:16:46 phil Exp phil $

indexes:

  0   1   4   9  16  
  1   2   5  10  17  
  4   5   8  13  20  
  9  10  13  18  25  
 16  17  20  25  32 

bin edges:  array([0.0, 8.0, 16.0, 24.0, 32.0], 'f') 

bin centers:  array([4.0, 12.0, 20.0, 28.0], 'f') 

bins:
array([0, 0, 0, 1, 2], 'i')
array([0, 0, 0, 1, 2], 'i')
array([0, 0, 1, 1, 2], 'i')
array([1, 1, 1, 2, 3], 'i')
array([2, 2, 2, 3, 4], 'i')

--------------------- blitz-support list --------------------------------
* To subscribe/unsubscribe: use the handy web form at
http://oonumerics.org/blitz/lists.html