通过Cython将C++向量传递给Numpy,无需复制并自动处理内存管理

Question

处理大型矩阵(NxM,1K <= N <= 20K&10K <= M <= 200K),我经常需要通过Cython将Numpy矩阵传递给C++来完成工作,这可以按预期工作而无需复制.

但是,有时我需要在C++中启动和预处理矩阵并将其传递给Numpy(Python 3.6).让我们假设矩阵是线性化的(因此大小为N*M,它是一维矩阵 - col/row major在这里无关紧要).根据这里的信息:在没有数据副本的情况下在Python中公开C计算数组并修改它以实现C++兼容性,我能够传递C++数组.

问题是如果我想使用标准向量而不是启动数组,我会得到分段错误.例如,考虑以下文件:

fast.h

#include <iostream>
#include <vector>

using std::cout; using std::endl; using std::vector;
int* doit(int length);

fast.cpp

#include "fast.h"
int* doit(int length) {
    // Something really heavy
    cout << "C++: doing it fast " << endl; 

    vector<int> WhyNot;

    // Heavy stuff - like reading a big file and preprocessing it
    for(int i=0; i<length; ++i)
        WhyNot.push_back(i); // heavy stuff

    cout << "C++: did it really fast" << endl;
    return &WhyNot[0]; // or WhyNot.data()
}

faster.pyx

cimport numpy as np
import numpy as np
from libc.stdlib cimport free
from cpython cimport PyObject, Py_INCREF

np.import_array()

cdef extern from "fast.h":
    int* doit(int length)

cdef class ArrayWrapper:
    cdef void* data_ptr
    cdef int size

    cdef set_data(self, int size, void* data_ptr):
        self.data_ptr = data_ptr
        self.size = size

    def __array__(self):
        print ("Cython: __array__ called")
        cdef np.npy_intp shape[1]
        shape[0] = <np.npy_intp> self.size
        ndarray = np.PyArray_SimpleNewFromData(1, shape,
                                               np.NPY_INT, self.data_ptr)
        print ("Cython: __array__ done")
        return ndarray

    def __dealloc__(self):
        print("Cython: __dealloc__ called")
        free(<void*>self.data_ptr)
        print("Cython: __dealloc__ done")


def faster(length):
    print("Cython: calling C++ function to do it")
    cdef int *array = doit(length)
    print("Cython: back from C++")
    cdef np.ndarray ndarray
    array_wrapper = ArrayWrapper()
    array_wrapper.set_data(length, <void*> array)
    print("Ctyhon: array wrapper set")
    ndarray = np.array(array_wrapper, copy=False)
    ndarray.base = <PyObject*> array_wrapper
    Py_INCREF(array_wrapper)
    print("Cython: all done - returning")
    return ndarray 

setup.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy 

ext_modules = [Extension(
    "faster", 
    ["faster.pyx", "fast.cpp"], 
    language='c++',
    extra_compile_args=["-std=c++11"],
    extra_link_args=["-std=c++11"]
)]

setup(
    cmdclass = {'build_ext': build_ext}, 
    ext_modules = ext_modules,
    include_dirs=[numpy.get_include()]
)

如果你用它构建

python setup.py build_ext --inplace

并运行Python 3.6解释器,如果输入以下内容,经过几次尝试后会出现seg错误.

>>> from faster import faster
>>> a = faster(1000000)
Cython: calling C++ function to do it
C++: doing it fast
C++: did it really fast
Cython: back from C++
Ctyhon: array wrapper set
Cython: __array__ called
Cython: __array__ done
Cython: all done - returning
>>> a = faster(1000000)
Cython: calling C++ function to do it
C++: doing it fast
C++: did it really fast
Cython: back from C++
Ctyhon: array wrapper set
Cython: __array__ called
Cython: __array__ done
Cython: all done - returning
Cython: __dealloc__ called
Segmentation fault (core dumped)

有几点需要注意:

• 如果你使用数组而不是vector(在fast.cpp中),这就像魅力一样!
• 如果你打电话faster(1000000)把结果放到除此之外的其他东西上variable a就行了.

如果您输入较小的数字,就像faster(10)您将获得更详细的信息,如:

Cython: calling C++ function to do it
C++: doing it fast
C++: did it really fast
Cython: back from C++
Ctyhon: array wrapper set
Cython: __array__ called
Cython: __array__ done
Cython: all done - returning
Cython: __dealloc__ called <--- Perhaps this happened too early or late?
*** Error in 'python': double free or corruption (fasttop): 0x0000000001365570 ***
======= Backtrace: =========
More info here ....

令人费解的是为什么阵列不会发生这种情况？无论!

我经常使用矢量,并希望能够在这些场景中使用它们.

Answer 1

我觉得@ FlorianWeimer的答案提供了一个体面的解决方案(分配vector和传递到你的C++函数),但它应该能够返回从一个载体doit,通过使用移动构造函数避免拷贝.

from libcpp.vector cimport vector

cdef extern from "<utility>" namespace "std" nogil:
  T move[T](T) # don't worry that this doesn't quite match the c++ signature

cdef extern from "fast.h":
    vector[int] doit(int length)

# define ArrayWrapper as holding in a vector
cdef class ArrayWrapper:
    cdef vector[int] vec
    cdef Py_ssize_t shape[1]
    cdef Py_ssize_t strides[1]

    # constructor and destructor are fairly unimportant now since
    # vec will be destroyed automatically.

    cdef set_data(self, vector[int]& data):
       self.vec = move(data)

    # now implement the buffer protocol for the class
    # which makes it generally useful to anything that expects an array
    def __getbuffer__(self, Py_buffer *buffer, int flags):
        # relevant documentation http://cython.readthedocs.io/en/latest/src/userguide/buffer.html#a-matrix-class
        cdef Py_ssize_t itemsize = sizeof(self.vec[0])

        self.shape[0] = self.vec.size()
        self.strides[0] = sizeof(int)
        buffer.buf = <char *>&(self.vec[0])
        buffer.format = 'i'
        buffer.internal = NULL
        buffer.itemsize = itemsize
        buffer.len = self.v.size() * itemsize   # product(shape) * itemsize
        buffer.ndim = 1
        buffer.obj = self
        buffer.readonly = 0
        buffer.shape = self.shape
        buffer.strides = self.strides
        buffer.suboffsets = NULL

然后,您应该能够将其用作:

cdef vector[int] array = doit(length)
cdef ArrayWrapper w
w.set_data(array) # "array" itself is invalid from here on
numpy_array = np.asarray(w)

Answer 2

从中返回时doit,WhyNot对象超出范围,并释放数组元素.这意味着它&WhyNot[0]不再是有效的指针.您需要将WhyNot对象存储在其他位置,可能在调用者提供的位置.

一种方法是doit分成三个函数,doit_allocate它们分配向量并返回指向它的指针,doit如前所述(但是带有一个参数,该参数接收指向, and解除分配矢量的预分配向量doit_free` 的指针.

像这样的东西:

vector<int> *
doit_allocate()
{
    return new vector<int>;
}

int *
doit(vector<int> *WhyNot, int length)
{
    // Something really heavy
    cout << "C++: doing it fast " << endl; 

    // Heavy stuff - like reading a big file and preprocessing it
    for(int i=0; i<length; ++i)
        WhyNot->push_back(i); // heavy stuff

    cout << "C++: did it really fast" << endl;
    return WhyNot->front();
}

void
doit_free(vector<int> *WhyNot)
{
    delete WhyNot;
}