Moving from Perl to Python

[Phil Austin]

jonathon <jblake at stamp-coin.com> writes:

> On 25 Sep 1999, Harry George wrote:
> 
> 	Other interesting/useful language translations might be
> 	
> 	_Numerical Fortran_/_Numerical C_/_Numerical Pascal_  	
> 			<< Or has this been done?  >>

Thanks to Paul Dubois and collaborators, my guess is that the
difficulty threshold for wrapping Fortran code using either CXX or
pyfort (http://xfiles.llnl.gov/python.htm) has become so low that
reimplementation probably isn't worth it (expecially since netlib
usually contains a free Fortran/C/C++ free alternative that could be
substituted).  For those who haven't seen the new technology, here's a
small example: the complete Python/CXX wrapper for the NR routines fasper.f
(Section 13.8 of the 2nd edition) and gasdev.f:

#include "Python.h"
#include "CXX_Objects.h"
#include "CXX_Extensions.h"
#include "CXX_Array.h"
#include <math.h>
#include <iostream>
#include <assert.h>

using namespace std;

extern "C" {
  void fasper_(float* x,float* y,int& n,float& ofac,float& hifac,float* wk1,
	       float* wk2,int& nwk,int& nout,int& jmax,float& prob);
  float gasdev_(int& idum);
}

static PyObject *
ex_fasper(PyObject* self, PyObject* args)
{
  Py_Initialize();
  import_array();
  PyImport_ImportModule("Numeric");
  Py::Tuple inargs(args);  
  Py::Array x(inargs[0]);
  assert(x.species()==PyArray_FLOAT);
  Py::Array y(inargs[1]);
  assert(y.species()==PyArray_FLOAT);
  int jint = Py::Int(inargs[2]);
  float ofacfloat = Py::Float(inargs[3]);
  float hifacfloat = Py::Float(inargs[4]);
  int MPint = Py::Int(inargs[5]);
  Py::Array wk1(MPint,PyArray_FLOAT);
  Py::Array wk2(MPint,PyArray_FLOAT);
  int nout;
  int jmax;
  float prob;
  fasper_( (float*) x.to_C(), (float*) y.to_C(), jint, ofacfloat,  
              hifacfloat ,(float*) wk1.to_C(),
	      (float*) wk2.to_C(), MPint ,nout,jmax,prob);
  Py::Dict output;
  output["wk1"]=wk1;
  output["wk2"]=wk2;
  output["nout"]=Py::Int(nout);
  output["jmax"]=Py::Int(jmax);
  output["prob"]=Py::Float(prob);
  return Py::new_reference_to(output);
}


static PyObject *
ex_gasdev(PyObject* self, PyObject* args)
{
  Py_Initialize();
  import_array();
  PyImport_ImportModule("Numeric");
  Py::Tuple inargs(args);  
  int idumint = Py::Int(inargs[0]);
  float ranNum=gasdev_(idumint);
  Py::Tuple output(2);
  output[0]=Py::Int(idumint);
  output[1]=Py::Float(ranNum);
  return Py::new_reference_to(output);
}


static Py::ExtensionModule* lombper;

extern "C" void initlombper();

void initlombper()
{
    // experimental initialization stuff
    lombper = new Py::ExtensionModule("lombper");
    lombper->add("fasper", ex_fasper, "calculate the Lomb periodogram");
    lombper->add("gasdev", ex_gasdev, "return a gaussian deviate");
    Py::Dict d = lombper->initialize();
}