特征数组的返回数组：从C++到Python的矩阵，不需要拷贝

Question

我有一些C++代码，用来生成和处理^ {< CD1>}矩阵的数组。 最后，我想在python中使用这些矩阵，并认为这可能是pybind11的工作

基本上，我想要回到python中的是两个嵌套列表/numpy数组 mat_a(I, 4, 4)和mat_b(J, K, 4, 4)。 因为我必须在C++中做很多线性代数的事情，我想使用EGIN，我使用的数据结构是 std::array<std::array<Eigen::Matrix4f, 2>, 3>>> mat_b // for J=3, K=2。 现在的问题是如何将其高效地应用于python

此外，我想对多个输入执行这些计算，结果是mat_a(N, I, 4, 4)和mat_b(N, J, K, 4, 4)。每个{{CD8}}的计算是独立的，但我认为在C++中，在{{CD8}}上写这个循环可能更快。另一方面，如果我们在C++中只有固定大小的数组，任务就会变得更容易，这个循环也可以移动到Python。p>

下面是我的问题的一些伪代码（I=5，J=3，K=2）：

// example.cpp
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <pybind11/stl.h>
#include <pybind11/functional.h>
#include <pybind11/stl_bind.h>

#include <array>
#include <vector>
#include <Eigen/Dense>


Eigen::Matrix4f get_dummy(){
    Eigen::Matrix4f mat_a;
    mat_a << 1, 2, 3, 4,
             5, 6, 7, 8,
             9, 8, 7, 6,
             5, 4, 3, 2;
    return mat_a;
}

std::pair< std::vector<std::array<Eigen::Matrix4f, 5> >,
           std::vector<std::array<std::array<Eigen::Matrix4f, 2>, 3> > >  get_matrices(std::vector<float> & x){

    std::vector<std::array<Eigen::Matrix4f, 5> > mat_a(x.size());
    std::vector< std::array< std::array< Eigen::Matrix4f, 2>, 3> > mat_b(x.size());

    //    for (u_int i=0; i< x.size(); i++)
    //        do_stuff(x[i], mat_a[i], mat_b[i]);
    mat_a[0][0] = get_dummy();

    return std::make_pair(mat_a, mat_b);
    }


PYBIND11_MODULE(example, m) {
    m.def("get_dummy", &get_dummy, pybind11::return_value_policy::reference_internal);
    m.def("get_matrices", &get_matrices, pybind11::return_value_policy::reference_internal);
}

我通过以下方式编译代码：

c++ -O3 -Wall -shared -std=c++14 -fPIC `python3 -m pybind11 --includes` example.cpp -o example`python3-config --extension-suffix`

然后在python中使用它：

import numpy as np
import example

x = np.zeros(1000)

mat_a, mat_b = get_matrices(x)

print(np.shape(mat_a))
print(np.shape(mat_b))
print(mat_a[0][0])

如果我只想返回一个Eigen::Matrix，它工作得很快，而且我可以说它不需要复制。但是当我尝试用std::array/std::vector嵌套Eigen:Matrices时，pybind返回一个numpy数组的嵌套列表，而不是一个多维数组。 这和预期的一样，我对它的工作原理印象深刻，但对我来说似乎相当缓慢，尤其是随着阵列尺寸的增加

问题是我如何改进这一点以获得多维numpy数组而不进行不必要的复制

有些路我试过了，但没有成功（对我来说，这并不意味着它们一般都不成功；我就是想不出来）：

• 使用Eigen::Tensor代替Eigen:Matrix的数组
• 在Python中创建矩阵并通过引用< /LI>将其传递给C++
• 为阵列构建自定义包装器<；数组<；矩阵4f，K>；，J>

Answer 1

如果您不太喜欢Eigen，另一种可能是xtensor（ found here）。我使用过他们的python绑定，在此之前给出了一个直接与python（found here）通信的示例。这将具有能够处理更大、多维阵列的优势。线性代数不会那么圆滑（在这里很难打败Eigen），但会类似于在numpy（np.dot(A,B)）中所做的

如果您想坚持使用Eigen，请注意使用STL有一些技术细节。由于std::array不再能够包含固定数量的矩阵，当您移动到std::vector时，您将遇到对齐问题（无可否认，我并不完全理解）。很快将为您获得一个xtensor的工作实现

Answer 2

您最好的选择可能是在python端创建数据，以便对其进行refcounted和垃圾收集

test.py

import example
import numpy as np

array = np.zeros((3, 2, 4, 4), 'f4')

example.do_math(array, 3, 2)
print(array[0, 0])

示例.cpp

#define PY_SSIZE_T_CLEAN
#include <Python.h>

#include <Eigen/Dense>

Eigen::Matrix4f get_dummy() {
    Eigen::Matrix4f mat_a;
    mat_a << 1, 2, 3, 4,
             5, 6, 7, 8,
             9, 8, 7, 6,
             5, 4, 3, 2;
    return mat_a;
}

PyObject * example_meth_do_math(PyObject * self, PyObject * args, PyObject * kwargs) {
    static char * keywords[] = {"array", "rows", "cols", NULL};

    PyObject * array;
    int rows, cols;

    if (!PyArg_ParseTupleAndKeywords(args, kwargs, "Oii", keywords, &array, &rows, &cols)) {
        return NULL;
    }

    Py_buffer view = {};
    if (PyObject_GetBuffer(array, &view, PyBUF_SIMPLE)) {
        return NULL;
    }

    Eigen::Matrix4f * ptr = (Eigen::Matrix4f *)view.buf;

    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            ptr[i * cols + j] = get_dummy();
        }
    }

    PyBuffer_Release(&view);
    Py_RETURN_NONE;
}

PyMethodDef module_methods[] = {
    {"do_math", (PyCFunction)example_meth_do_math, METH_VARARGS | METH_KEYWORDS, NULL},
    {},
};

PyModuleDef module_def = {PyModuleDef_HEAD_INIT, "example", NULL, -1, module_methods};

extern "C" PyObject * PyInit_example() {
    PyObject * module = PyModule_Create(&module_def);
    return module;
}

setup.py

from setuptools import Extension, setup

ext = Extension(
    name='example',
    sources=['./example.cpp'],
    extra_compile_args=['-fpermissive'],
    include_dirs=['.'], # add the path of Eigen
    library_dirs=[],
    libraries=[],
)

setup(
    name='example',
    version='0.1.0',
    ext_modules=[ext],
)

从这里开始，添加第二个参数并使用这两个数组进行计算应该很简单

您可以使用python setup.py develop来构建它

如果要分发它，可以使用python setup.py bdist_wheel创建一个控制盘文件

我使用numpy来创建数据，这确保了数据的底层内存是C连续的

该示例保持简单，并使用Matrix4f指针迭代3x2矩阵数组。您可以随意将ptr转换为Eigen::Array<Eigen::Matrix4f>, 3, 2>。无法将其强制转换为std::vector，因为std::vector的内部数据包含指针

请注意std::vector<std::array<...>>内存中没有单个连续数组。改用Eigen::Array

编辑：

下面是一个使用Eigen{}{}的函数：

PyObject * example_meth_do_math(PyObject * self, PyObject * args, PyObject * kwargs) {
    static char * keywords[] = {"array", NULL};

    PyObject * array;

    if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O", keywords, &array)) {
        return NULL;
    }

    Py_buffer view = {};
    if (PyObject_GetBuffer(array, &view, PyBUF_SIMPLE)) {
        return NULL;
    }

    Eigen::Map<Eigen::Array<Eigen::Matrix4f, 2, 3>> array_map((Eigen::Matrix4f *)view.buf, 2, 3);

    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < 3; ++j) {
            array_map(i, j) = get_dummy();
        }
    }

    PyBuffer_Release(&view);
    Py_RETURN_NONE;
}