将大型复杂数组从Python传递到C++ - 我最好的选择是什么？

Question

2017/06/13编辑:我尝试按照建议使用提升,但在花了超过3天试图让它进行编译和链接,并且失败后,我认为愚蠢的痛苦方式可能是最快且不那么痛苦的.. ..所以现在我的代码只是保存了一堆巨大的文本文件(拆分数组和文件中数字的复杂/虚部),然后C++读取.优雅......不......有效......是的.

---

我有一些科学代码,目前用Python编写,在循环中通过数字3d集成步骤减慢速度.为了克服这个问题,我在C++中重写了这一特定步骤.(Cython等不是一个选项).

简而言之:我希望尽可能方便,轻松地将几个非常大的复数数组从python代码传输到C++集成器.我可以使用文本或二进制文件手动和痛苦地做到这一点 - 但在我开始之前,我想知道我是否有更好的选择？

我正在使用Visual Studio for C++和anaconda for python(不是我的选择!)

是否有任何文件格式或方法可以快速方便地从python中保存一组复数,然后在C++中重新创建它？

非常感谢,本

Answer 1

我多次使用的一个简单的解决方案是将你的"C++方面"构建为dll(= Linux/OS X上的共享对象),提供一个简单的,类似C的入口点(直整数,指针和co.,没有STL的东西)并传递数据ctypes.

这避免了boost/SIP/Swig/...构建的噩梦,可以保持零拷贝(使用ctypes可以直接指向你的numpy数据)并允许你做任何你想做的事情(特别是在构建方面 - 没有friggin'distutils,没有提升,没有任何东西 - 用C++方面的任何可以构建类似C的dll来构建它.它还具有使用其他语言调用C++算法的良好副作用(实际上任何语言都有某种方式可以与C库进行交互).

这是一个快速的人为例子.C++方面只是:

extern "C" {
double sum_it(double *array, int size) {
    double ret = 0.;
    for(int i=0; i<size; ++i) {
        ret += array[i];
    }
    return ret;
}
}

这必须编译为dll(在Windows上)或.so(在Linux上),确保导出sum_it函数(使用gcc自动,需要.defVC++文件).

在Python方面,我们可以有一个包装器

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))

这可以确保数据正确封送; 然后调用该函数就像微不足道的那样

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.]))

有关如何使用它的更多信息,请参阅ctypes文档.另请参阅numpy中的相关文档.

对于复杂的数字,情况稍微复杂一些,因为在ctypes中没有内置的东西; 如果我们想std::complex<double>在C++端使用(几乎可以保证在numpy复杂布局中工作正常,即两个双精度序列),我们可以将C++端编写为:

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;
}
}

然后,在Python方面,我们必须复制c_complex布局以检索返回值(或者能够构建没有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

继承ctypes.Structure允许ctypes编组魔法,这是根据_fields_成员执行的; 构造函数和额外的方法只是为了易于在Python端使用.

然后,我们必须告诉ctypes返回类型

_sum_it_cplx = cdll.sum_it_cplx
_sum_it_cplx.restype = c_complex

最后以与前一个类似的方式编写我们的包装器:

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()

如上所述进行测试

# 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]))

产生预期结果:

(3+3j)
(3+3j)

Answer 2

我看到 OP 现在已经有一年多了，但我最近使用本机 Python-C/C++ API 及其 Numpy-C/C++ 扩展解决了类似的问题，并且由于我个人不喜欢使用 ctypes，原因有多种（例如，复杂的数字解决方法、混乱的代码），也不是 Boost，想为未来的搜索者发布我的答案。

Python-C API 和 Numpy-C API 的文档都非常广泛（尽管一开始有点让人不知所措）。但在一两次成功之后，编写本机 C/C++ 扩展就变得非常容易。

下面是一个可以从 Python 调用的 C++ 函数示例。numpy.double它集成了实数或复数 (或) 类型的 3D numpy 数组numpy.cdouble。.so该函数将通过模块的DLL() 导入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);
    }

};

只需将其放入源文件中（我们将其cintegrate.cxx与模块定义内容一起调用，插入到底部：

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
};

然后通过类似于 Walter 的 boost 示例的方式进行编译，setup.py只需进行一些明显的更改 - 将file.cc其替换为我们的 file ，删除 boost 依赖项，并确保包含 的cintegrate.cxx路径。"numpy/arrayobject.h"

在 python 中，你可以像这样使用它：

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)

此特定代码尚未经过测试，但大部分是从工作代码复制的。