within a framework of a simulation of some biophysical model I have a C++ class that implements my model with a member function that needs return a PyArrayObject*. The class is defined in the header file ng_networks.h which reads something like:
#include <stdio.h>
#include <math.h>
#include <Python.h>
#include <numpy/arrayobject.h> /* Numpy as seen from C */
struct ngn{...}; /* Structure of model parameter */
class ngn_rate{ /* Model Class */
int NEQ;
ngn pars;
double *y,*dy;
public:
ngn_rate(); // Constructor
~ngn_rate(); // Destructor
PyObject* bifsystem();
}
ngn_rate::ngn_rate(){
// Set ngn pars default values
...
// Allocate memory for model variables and RHS of equations
y = (double*)calloc(NEQ,sizeof(double));
dy = (double*)calloc(NEQ,sizeof(double));
}
ngn_rate::~ngn_rate{
free(y);
free(dy);
}
PyObject* ngn_rate::bifsystem(){
long int NUMEL = NEQ; // PyArray creation function requires (long int*)
// Does some operations on y,dy
...
// Build PyObject* from dy
// Create Python Array Object...
PyObject* out_array = PyArray_SimpleNew(1,&NUMEL,NPY_DOUBLE);
// ... and make C pointer to point to it
double* dy_aux = (double*)((PyArrayObject*)out_array)->data;
// Copy dy into PyObject out_array
for(int i=0;i<NUMEL;i++) dy_aux[i] = dy[i];
return out_array;
}
As you may guess this class is eventually called from a Python module. With this regard I am using scipy.weave to interface my C code with Python. So the calling Python module looks something like:
def ngn_rate_py():
support_code = """
#include <Python.h>
#include "ng_networks.h"
"""
source_files = [...] # A list of C/CPP source file names
libs = ['m']
dirs = [...] # A list of #include dirs
# My C code to interface with Python
code = """
//Initialize Model
ngn_rate network();
//Return dy_dt
return_val = network.bifsystem();
"""
vars = []
dy = weave.inline(code,
vars,
support_code = support_code,
sources = source_files,
libraries = libs,
library_dirs = dirs,
include_dirs = dirs,
runtime_library_dirs = dirs,
type_converters = converters.blitz,
compiler = 'gcc',
extra_compile_args = ['-std=c++11'],
force = 1)
return dy
When I run the above module unfortunately it produces a segmentation fault. After some debugging and trial-and-error I figured out that the problem is caused for some reason by the initialization of PyObject* out_array in ng_networks.h. Indeed when I create the PyObject* in C code in weave it works perfectly: that is I modify the ngn_rate::bifsystem() class member function so that it returns a double* and then I build the PyObject* from the latter within the weave interface:
class ngn_rate{ /* Model Class */
...
public:
...
double* bifsystem();
}
double* ngn_rate::bifsystem(){
long int NUMEL = NEQ; // PyArray creation function requires (long int*)
// Does some operations on y,dy
...
return dy;
}
And then in the weave interface:
def ngn_rate_py():
support_code = """
#include <Python.h>
#include "ng_networks.h"
"""
code = """
//Initialize Model
ngn_rate network();
//Create temporary dy
double *dy = network.bifsystem();
//Create PyObject*
PyObject* out_array = PyArray_SimpleNew(1, &NUMEL, NPY_DOUBLE);
double* dy_aux = (double*) ((PyArrayObject*) out_array)->data;
//Assign dy to PyObject
for(int i=0;i<NUMEL;i++) dy_aux[i]=dy[i];
return_val = out_array;
I cannot figure out why the above works whereas I get a segmentation fault if I make PyObject* be returned by my class. Noteworthy is that in some other code I had just a static/non-member C-function that returned a PyObject* that when called from weave worked just fine. So my guess is that there are some issues with PyObject* used within C class object. But I cannot figure out what. And although I can work with the above code creating the PyObject* within the weave interface, for portability I'd rather have it directly provided by my class ngn_rate.
Thanks in advance for your feedback.
M
