5

2017/06/13 EDIT: I tried using boost as was suggested, but after spending more than 3 days trying to get it to compile and link, and failing, I decided that the stupid painful way was probably the fastest and less painfull.... so now my code just saves a mess of gigantic text files (splitting arrays and the complex/ imaginary parts of the numbers across files) that C++ then reads. Elegant... no.... effective... yes.


I have some scientific code, currently written in Python, that is being slowed down by a numerical 3d integration step within a loop. To overcome this I am re-writing this particular step in C++. (Cython etc is not an option).

Long story short: I want to transfer several very large arrays of complex numbers from the python code to the C++ integrator as conveniently and painlessly as possible. I could do this manually and painfully using text or binary files - but before I embark on this, I was wondering if I have any better options?

I'm using visual studio for C++ and anaconda for python (not my choice!)

Is there any file format or method that would make it quick and convenient to save an array of complex numbers from python and then recreate it in C++?

Many thanks, Ben

Ben
  • 353
  • 1
  • 6
  • 14
  • Why not try to make the algorithm more efficient? If it is inefficient in Python, it will still be in inefficient in C++ (in theory). –  Jun 07 '17 at 14:44
  • 1
    Is doing your integration in Python with Numpy an option? It's quite fast at processing arrays because it uses compiled code rather than the Python interpreter. – PM 2Ring Jun 07 '17 at 14:52
  • I've tried with Numpy and come to the conclusion that, for me, writing this one step in C++ will ultimately save me time. What I really want is low level memory access and I am finding python/numpy a little opaque in this regard. Maybe it is possible, but I have limited time to get this done, and I have previous experience writing C++ code - so I'd rather just do it this way than spend a lot of time learning to do it in Numpy. Once the data is imported I know what to do in C++. – Ben Jun 07 '17 at 15:02
  • 1
    @InternetAussie There is no efficiency in python. numpy gains efficiency from being implemented in C, not python. The OP simply wants to exploit the same. – Walter Jun 07 '17 at 15:03
  • Also, this step is going to be repeated many hundreds of thousands of times. Any speed increase, even if it is only 20% or something is potentially a big plus. I'm pretty set on using C++. Thanks for the suggestions though! – Ben Jun 07 '17 at 15:03
  • Ben: So your question is really about "Passing big arrays of complex numbers from Python to C++", right? – martineau Jun 07 '17 at 15:06
  • Fair enough. So what's the datatype of your "arrays of complex numbers from the python code"? Are they plain Python lists? If so, Numpy can convert them efficiently to arrays that are contiguous in memory, using either C or Fortran array conventions. And of course it can also do the reverse conversion. So Numpy may be useful for you, even if you don't want to use it for performing the integration. – PM 2Ring Jun 08 '17 at 09:20
  • 1
    I think you're **asking the wrong question**. You *assume* that the solution to your problem is to pass python data to C++, but perhaps (and likely) the solution is to use numpy. I strongly suggest, you ask another question, where you sketch the actual problem, i.e. give an MVCE (in python). – Walter Jun 09 '17 at 09:21
  • My arrays are either numpy arrays or mpmath arrays - either way, I can convert them to anything else before they are needed. – Ben Jun 09 '17 at 13:01
  • Walter - you may be right, but as it stands, I'm not familiar with all the ins and outs of numpy and I know how to do this in C++, once I have the data. I would like to use C++ so I can know exactly what is going on at the lowest levels, rather than having to study numpy and work out what it's doing behind the scenes. The C++ part of this is no issue for me. I just want to import the data conveniently. – Ben Jun 09 '17 at 13:04

3 Answers3

9

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" {
double sum_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

def sum_it(l):
    if isinstance(l, np.ndarray) and l.dtype == np.float64 and len(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-copy
print 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):

class c_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 complex
        self.real = pycomplex.real
        self.imag = pycomplex.imag

    def to_complex(self):
        # Convert to Python complex
        return self.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:

def sum_it_cplx(l):
    if isinstance(l, np.ndarray) and l.dtype == np.complex and len(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-copy
print summer.sum_it_cplx(np.array([1. + 0.j, 0 + 1.j, 2 + 2.j]))

yields the expected results:

(3+3j)
(3+3j)
Matteo Italia
  • 123,740
  • 17
  • 206
  • 299
  • Could you give a simple example (or point to one)? What is `ctypes`? – Walter Jun 13 '17 at 10:49
  • @Walter: ctypes is the standard module to interact from Python code with C libraries; see [its documentation](https://docs.python.org/2/library/ctypes.html). Also, I added a simple example (with compatibility both with regular lists - with the unavoidable copy - and zero-copy compatibility with numpy arrays). – Matteo Italia Jun 13 '17 at 16:40
  • Glad it worked! Did you obtain the speedup you were after in the end? – Matteo Italia Jun 16 '17 at 05:17
  • This solution works great for real arrays but it does not answer the question of how to do it starting from complex numpy arrays in python. Any ideas how to this for *complex* arrays? – co9olguy Nov 09 '17 at 20:32
  • 1
    @co9olguy: I just updated the example on my git, I'll add it to the question as well. – Matteo Italia Nov 10 '17 at 08:28
  • 1
    @peedurrr: it's actually on Bitbucket - see the first link in the post: https://bitbucket.org/mitalia/ctypes_dll_example – Matteo Italia Jul 02 '20 at 13:15
  • 1
    @MatteoItalia : getting segmentation fault while processing numpy array – Devil Nov 28 '22 at 06:18
6

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 API
    double 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.");
        return NULL;
    }

    // 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 dims
    if (PyArray_NDIM(yarr) != 3) {
        Py_CLEAR(yarr);
        PyErr_SetString(PyExc_ValueError, "Expected 3 dimensional integrand");
        return NULL;
    }

    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 reference
    if (iscomplex) {
        return Py_BuildValue("D", &result);
    } else {
        return Py_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 */
};


static struct PyModuleDef module = {
   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.

v.chaplin
  • 617
  • 6
  • 8
  • If you get segmentation fault, you might need to add `import_array();` before `int DTYPE = PyArray_ObjectType(argy, NPY_FLOAT); ` (see https://stackoverflow.com/a/66245936) – David Střelák Dec 02 '22 at 13:52
1

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.py

from 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.

Walter
  • 44,150
  • 20
  • 113
  • 196
  • @martineau why not? – Walter Jun 07 '17 at 15:14
  • My comment was before you edited it—obviously,wise guy. – martineau Jun 07 '17 at 15:32
  • 1
    Note that this answer assumes that the input are ndarrays, which OP explicitly said are not (well, OP said they aren't using numpy). – Andras Deak -- Слава Україні Jun 08 '17 at 12:03
  • @AndrasDeak I assume nothing. I merely demonstrate how one may pass data from python to C++, exemplifying with an `numpy.ndarray`, but that is only an example. However, the OP would be daft to not use a `numpy.ndarray`. Note also that the OP does not mention numpy and the comment made does not imply that numpy arrays are not used. – Walter Jun 09 '17 at 09:09
  • I am using numpy! I've just found it to be a little slow - or perhaps I should say hard (for me with my current knowledge) to optimise. This is why I would rather use C++. – Ben Jun 09 '17 at 13:10
  • I am trying to implement your solution now. – Ben Jun 09 '17 at 13:11
  • I had to give up on this. Compiling boost was a nightmare, and even when I managed that, my code wouldnt compile without errors. I gave up after three days and just resorted to a million horrible text files. Thanks for the tip about boost though - I may come back to it for future projects where I'm not stuck using visual studio – Ben Jun 13 '17 at 07:59