[C++-SIG] blitz++ and NumPy, revision 1

[Phil Austin]

>>>>> "Paul" == Paul F Dubois <dubois1 at llnl.gov> writes:

    Paul> a. In C++ certain things are not inherited. These include
    Paul> the constructors and copy constructors, and I think the
    Paul> assignment operator. In the CXX documentation it gives you a
    Paul> list in the comments of the ones you have to reimplement in
    Paul> a child. 

     PA> Thanks, I'll play around with this after I get this pile of finals
         graded.


Actually, I could live with something like the following, which just
pairs the Blitz and Py:Arrays together in a structure.
Corrections/suggestions are welcome.

Regards, Phil



___________________

#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
{
public:
  Py::Array Pn;
  blitz::Array<T,Ndim> Bz;
  TiedArray(blitz::TinyVector<int,Ndim> dimens);
};


template<> class TiedArray<float,2>
{
public:
  Py::Array Pn;
  blitz::Array<float,2> Bz;
  TiedArray(blitz::TinyVector<int,2> dimens):Pn(dimens,PyArray_FLOAT),Bz((float*) Pn.to_C(),dimens) {};  
  TiedArray(Py::seqref<Py::Object> pyob):Pn(pyob),Bz((float*) Pn.to_C(),
                            blitz::shape(Pn.dimension(1),Pn.dimension(2)))
    {
      assert(Pn.rank()==2);
      assert(Pn.species()==PyArray_FLOAT);
    };
};


static PyObject *
ex_strucout(PyObject* self, PyObject* args)
{
  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();
}