[Numpy-discussion] Options for wrapping C and C++ code for use with Numeric

Philip Austin paustin at eos.ubc.ca
Fri Nov 4 09:21:23 CST 2005


I'm planning to update my boost numeric wrappers
(http://www.eos.ubc.ca/research/clouds/num_util.html) to scipy_core as
soon as things settle down.  Here's an example of some fortran
program wrapped by boost that shows a couple of nice features:
i) mirroring of python types in C++, with indexing, etc., 
ii) transparent memory management (all increfs and decrefs are handled by
boost shared pointers) and iii) docstrings.  I'd be happy to contribute to a
page with some SWIG and boost examples for similar code fragments.


#define PY_ARRAY_UNIQUE_SYMBOL PyArrayHandle
#include <num_util.h>
#include <iostream> 
#include "boost/scoped_array.hpp"

namespace { const char* rcsid = ("@(#) $Id: thermwrap.C,v 1.2 2005/09/05 22:47:57 phil Exp $:"); };

using namespace std;
namespace py = boost::python;
namespace nbpl = num_util;

typedef int wave_unit;

extern "C" {
  void ptqv_(float& p,float& t, float& qv, float& psp, float& tsp, float& qsp, float& thsp, 
                   float& thesp);
  void ptpsp_(float& p,float& t, float& qv, float& psp, float& tsp, float& qsp, 
              float& thpt,float& thes, float& thsp, float& thesp);
  void mccla_(const char* atmos, float* z, float* p, float* t,float* rvden, 
                 float* o3den, float* den, int& strnlength);
  void thinv_(float& p,float& t,float& thsp);
  void theinv_(float& p,float& t, float& qsp, float& thsp, float& thesp);
}


py::dict ptqv(float p, float t, float qv)
{
  float psp,tsp,qsp,thsp,thesp;
  ptqv_(p,t,qv,psp,tsp,qsp,thsp,thesp);
  py::dict result;
  result["psp"]=psp;
  result["tsp"]=tsp;
  result["qsp"]=qsp;
  result["thsp"]=thsp;
  result["thesp"]=thesp;
  return result;
}

py::dict ptpsp(float p, float t, float psp)
{
  float qv,tsp,qsp,thpt,thes,thsp,thesp;
  ptpsp_(p,t,qv,psp,tsp,qsp,thpt,thes,thsp,thesp);
  py::dict result;
  result["qv"]=qv;
  result["tsp"]=tsp;
  result["qsp"]=qsp;
  result["thpt"]=thpt;
  result["thes"]=thes;
  result["thsp"]=thsp;
  result["thesp"]=thesp;
  return result;
}

py::dict thinv(float p, float thsp)
{
  float t;
  thinv_(p,t,thsp);
  py::dict result;
  result["t"]=t;
  return result;
}

py::dict theinv(float p, float t, float the)
{
  //t input is first guess
  float qv,thsp;
  theinv_(p,t,qv,thsp,the);
  py::dict result;
  result["t"]=t;
  result["qsp"]=qv;
  result["thsp"]=thsp;
  return result;
}

py::numeric::array mcclatchey(string atmos)
{
  int strlength=atmos.length();
  int length=33;
  boost::scoped_array<float> z(new float[length]);
  boost::scoped_array<float> p(new float[length]);
  boost::scoped_array<float> t(new float[length]);
  boost::scoped_array<float> rvden(new float[length]);
  boost::scoped_array<float> o3den(new float[length]);
  boost::scoped_array<float> den(new float[length]);
  mccla_(atmos.c_str(), z.get(), p.get(), t.get(), rvden.get(), o3den.get(), den.get(),strlength );
  std::vector<int> dims;
  dims.push_back(6);
  dims.push_back(33);
  py::numeric::array standsound(nbpl::makeNum(dims, PyArray_DOUBLE));
  double* sndPtr=(double*) nbpl::data(standsound);
  int index;
  for(int i=0;i<33;++i){
    index=i;
    sndPtr[index]=z[i];
    index=33 + i;
    sndPtr[index]=p[i];
    index=2*33 + i;
    sndPtr[index]=t[i];
    index=3*33 + i;
    sndPtr[index]=rvden[i];
    index=4*33 + i;
    sndPtr[index]=o3den[i];
    index=5*33 + i;
    sndPtr[index]=den[i];
  }
  return standsound;
}

BOOST_PYTHON_MODULE(thermwrap)
{
  using namespace boost::python;
  import_array();
  numeric::array::set_module_and_type("Numeric", "ArrayType");
  scope().attr("RCSID") = rcsid;
  scope().attr("__doc__") = "wrappers for BettsThermo.f routines: ptqv, ptpsp";
  string docstring("input: pressure (hPa), temp (K), qv (kg/kg)\n");
  docstring += "output psp (hPa), tsp (K), qsp (kg/kg), thsp (K), thesp (K)";
  def("ptqv", ptqv,docstring.c_str());
  docstring="input: pressure (hPa), temp (K), psp (hPa)\n";
  docstring += "output: qv (kg/kg), tsp (K), qsp (kg/kg), thpt (theta(p,t)), thes (K), \n";
  docstring +="thsp (theta(psl,tsl), thesp equiv. potential temp.(theta-e at sl)";
  def("ptpsp",ptpsp,docstring.c_str());
  def("thinv",thinv,"thinv");
  def("theinv",theinv,"theinv");
  def("mcclatchey",mcclatchey,"mcclatchey");
}









More information about the Numpy-discussion mailing list