Passing Big Complex Arrays From Python To C++ - What's My Best Option?
Solution 1:
An easy solution that I used many times is to build your "C++ side" as a dll (=shared object on Linux/OS X), provide a simple, C-like entrypoint (straight integers, pointers & co., no STL stuff) and pass the data through ctypes.
This avoids boost/SIP/Swig/... build nightmares, can be kept zero-copy (with ctypes you can pass a straight pointer to your numpy data) and allow you to do whatever you want (especially on the build-side - no friggin' distutils, no boost, no nothing - build it with whatever can build a C-like dll) on the C++ side. It has also the nice side-effect of having your C++ algorithm callable from other languages (virtually any language has some way to interface with C libraries).
Here's a quick artificial example. The C++ side is just:
extern"C" {
doublesum_it(double *array, int size){
double ret = 0.;
for(int i=0; i<size; ++i) {
ret += array[i];
}
return ret;
}
}
This has to be compiled to a dll (on Windows) or a .so (on Linux), making sure to export the sum_it function (automatic with gcc, requires a .def file with VC++).
On the Python side, we can have a wrapper like
import ctypes
import os
import sys
import numpy as np
path = os.path.dirname(__file__)
cdll = ctypes.CDLL(os.path.join(path, "summer.dll"if sys.platform.startswith("win") else"summer.so"))
_sum_it = cdll.sum_it
_sum_it.restype = ctypes.c_double
defsum_it(l):
ifisinstance(l, np.ndarray) and l.dtype == np.float64 andlen(l.shape)==1:
# it's already a numpy array with the right features - go zero-copy
a = l.ctypes.data
else:
# it's a list or something else - try to create a copy
arr_t = ctypes.c_double * len(l)
a = arr_t(*l)
return _sum_it(a, len(l))
which makes sure that the data is marshaled correctly; then invoking the function is as trivial as
import summer
import numpy as np
# from a list (with copy)print summer.sum_it([1, 2, 3, 4.5])
# from a numpy array of the right type - zero-copyprint summer.sum_it(np.array([3., 4., 5.]))
See the ctypes documentation for more information on how to use it. See also the relevant documentation in numpy.
For complex numbers, the situation is slightly more complicated, as there's no builtin for it in ctypes; if we want to use std::complex<double> on the C++ side (which is pretty much guaranteed to work fine with the numpy complex layout, namely a sequence of two doubles), we can write the C++ side as:
extern"C" {
std::complex<double> sum_it_cplx(std::complex<double> *array, int size){
std::complex<double> ret(0., 0.);
for(int i=0; i<size; ++i) {
ret += array[i];
}
return ret;
}
}
Then, on the Python side, we have to replicate the c_complex layout to retrieve the return value (or to be able to build complex arrays without numpy):
classc_complex(ctypes.Structure):# Complex number, compatible with std::complex layout
_fields_ = [("real", ctypes.c_double), ("imag", ctypes.c_double)]
def__init__(self, pycomplex):
# Init from Python complexself.real = pycomplex.real
self.imag = pycomplex.imag
defto_complex(self):
# Convert to Python complexreturnself.real + (1.j) * self.imag
Inheriting from ctypes.Structure enables the ctypes marshalling magic, which is performed according to the _fields_ member; the constructor and extra methods are just for ease of use on the Python side.
Then, we have to tell ctypes the return type
_sum_it_cplx = cdll.sum_it_cplx
_sum_it_cplx.restype = c_complex
and finally write our wrapper, in a similar fashion to the previous one:
defsum_it_cplx(l):
ifisinstance(l, np.ndarray) and l.dtype == np.complexandlen(l.shape)==1:
# the numpy array layout for complexes (sequence of two double) is already# compatible with std::complex (see https://stackoverflow.com/a/5020268/214671)
a = l.ctypes.data
else:
# otherwise, try to build our c_complex
arr_t = c_complex * len(l)
a = arr_t(*(c_complex(r) for r in l))
ret = _sum_it_cplx(a, len(l))
return ret.to_complex()
Testing it as above
# from a complex list (with copy)print summer.sum_it_cplx([1. + 0.j, 0 + 1.j, 2 + 2.j])
# from a numpy array of the right type - zero-copyprint summer.sum_it_cplx(np.array([1. + 0.j, 0 + 1.j, 2 + 2.j]))
yields the expected results:
(3+3j)
(3+3j)
Solution 2:
I see the OP is over a year old now, but I recently addressed a similar problem using the native Python-C/C++ API and its Numpy-C/C++ extension, and since I personally don't enjoy using ctypes for various reasons (e.g., complex number workarounds, messy code), nor Boost, wanted to post my answer for future searchers.
Documentation for the Python-C API and Numpy-C API are both quite extensive (albeit a little overwhelming at first). But after one or two successes, writing native C/C++ extensions becomes very easy.
Here is an example C++ function that can be called from Python. It integrates a 3D numpy array of either real or complex (numpy.double or numpy.cdouble) type. The function will be imported through a DLL (.so) via the module cintegrate.so.
#include"Python.h"#include"numpy/arrayobject.h"#include<math.h>static PyObject * integrate3(PyObject * module, PyObject * args){
PyObject * argy=NULL; // Regular Python/C API
PyArrayObject * yarr=NULL; // Extended Numpy/C APIdouble dx,dy,dz;
// "O" format -> read argument as a PyObject type into argy (Python/C API)if (!PyArg_ParseTuple(args, "Oddd", &argy,&dx,&dy,&dz)
{
PyErr_SetString(PyExc_ValueError, "Error parsing arguments.");
returnNULL;
}
// Determine if it's a complex number array (Numpy/C API)int DTYPE = PyArray_ObjectType(argy, NPY_FLOAT);
int iscomplex = PyTypeNum_ISCOMPLEX(DTYPE);
// parse python object into numpy array (Numpy/C API)
yarr = (PyArrayObject *)PyArray_FROM_OTF(argy, DTYPE, NPY_ARRAY_IN_ARRAY);
if (yarr==NULL) {
Py_INCREF(Py_None);
return Py_None;
}
//just assume this for 3 dimensional array...you can generalize to N dimsif (PyArray_NDIM(yarr) != 3) {
Py_CLEAR(yarr);
PyErr_SetString(PyExc_ValueError, "Expected 3 dimensional integrand");
returnNULL;
}
npy_intp * dims = PyArray_DIMS(yarr);
npy_intp i,j,k,m;
double * p;
//initialize variable to hold result
Py_complex result = {.real = 0, .imag = 0};
if (iscomplex) {
for (i=0;i<dims[0];i++)
for (j=0;j<dims[1];j++)
for (k=0;k<dims[1];k++) {
p = (double*)PyArray_GETPTR3(yarr, i,j,k);
result.real += *p;
result.imag += *(p+1);
}
} else {
for (i=0;i<dims[0];i++)
for (j=0;j<dims[1];j++)
for (k=0;k<dims[1];k++) {
p = (double*)PyArray_GETPTR3(yarr, i,j,k);
result.real += *p;
}
}
//multiply by step size
result.real *= (dx*dy*dz);
result.imag *= (dx*dy*dz);
Py_CLEAR(yarr);
//copy result into returnable type with new referenceif (iscomplex) {
returnPy_BuildValue("D", &result);
} else {
returnPy_BuildValue("d", result.real);
}
};
Simply put that into a source file (we'll call it cintegrate.cxx along with the module definition stuff, inserted at the bottom:
static PyMethodDef cintegrate_Methods[] = {
{"integrate3", integrate3, METH_VARARGS,
"Pass 3D numpy array (double or complex) and dx,dy,dz step size. Returns Reimman integral"},
{NULL, NULL, 0, NULL} /* Sentinel */
};
staticstructPyModuleDefmodule = {
PyModuleDef_HEAD_INIT,
"cintegrate", /* name of module */NULL, /* module documentation, may be NULL */-1, /* size of per-interpreter state of the module,
or -1 if the module keeps state in global variables. */
cintegrate_Methods
};
Then compile that via setup.py much like Walter's boost example with just a couple obvious changes- replacing file.cc there with our file cintegrate.cxx, removing boost dependencies, and making sure the path to "numpy/arrayobject.h" is included.
In python then you can use it like:
import cintegrate
import numpy as np
arr = np.random.randn(4,8,16) + 1j*np.random.randn(4,8,16)
# arbitrary step size dx = 1., y=0.5, dz = 0.25
ans = cintegrate.integrate3(arr, 1.0, 0.5, .25)
This specific code hasn't been tested but is mostly copied from working code.
Solution 3:
Note added in edit.
As mentioned in the comments, python itself, being an interpreted language, has little potential for computational efficiency. So in order to make python scripts efficient, one must use modules which aren't all interpreted, but under the hood call compiled (and optimized) code written in, say, C/C++. This is exactly what numpy does for you, in particular for operations on whole arrays.
Therefore, the first step towards efficient python scripts is the usage of numpy. Only the second step is to try to use your own compiled (and optimized) code. Therefore, I have assumed in my example below that you were using numpy to store the array of complex numbers. Everything else would be ill-advised.
There are various ways in which you can access python's original data from within a C/C++ program. I personally have done this with boost.Python, but must warn you that the documentation and support are lousy at best: you're pretty much on your own (and stack overflow, of course).
For example your C++ file may look like this
// file.cc#include<boost/python.hpp>#include<boost/python/numpy.hpp>namespace p = boost::python;
namespace n = p::numpy;
n::ndarray func(const n::ndarray&input, double control_variable){
/*
your code here, see documentation for boost python
you pass almost any python variable, doesn't have to be numpy stuff
*/
}
BOOST_PYTHON_MODULE(module_name)
{
Py_Initialize();
n::initialize(); // only needed if you use numpy in the interface
p::def("function", func, "doc-string");
}
to compile this, you may use a python script such as
# setup.pyfrom distutils.core import setup
from distutils.extension import Extension
module_name = Extension(
'module_name',
extra_compile_args=['-std=c++11','-stdlib=libc++','-I/some/path/','-march=native'],
extra_link_args=['-stdlib=libc++'],
sources=['file.cc'],
libraries=['boost_python','boost_numpy'])
setup(
name='module_name',
version='0.1',
ext_modules=[module_name])
and run it as python setup.py build, which will create an appropriate .so file in a sub-directory of build, which you can import from python.
Post a Comment for "Passing Big Complex Arrays From Python To C++ - What's My Best Option?"