Python/C++混合编程利器Pybind11实践
混合编程简介
对于Python这种脚本语言，编程方便，便于快速验证算法，但是在实际开发过程中，免不了涉及到混合编程，比如在一些计算密集型的情景下面，还是需要使用C/C++来完成。
通常我们在Python平台上进行算法编程工作，测试正确后将性能有瓶颈的模块或者需要并行的函数部分用C++重新实现，并通过python import调用动态链接库（.so/.pyd等）实现本地代码加速。
pybind11 是一个轻量级的只包含头文件的库，用于 Python 和 C++ 之间接口转换，可以为现有的 C++ 代码创建 Python 接口绑定。其目的和语法类似于 Boost.Python 库，但是Boost库十分庞大而复杂，而pybind11十分轻量级。pybind11 通过 C++ 编译时的检查来推断类型信息，来最大程度地减少传统拓展 Python 模块时繁杂的样板代码。已经实现了所有主要的 Python 类型到C++的转换，包括handle, object, bool_, int_, float_, str, bytes, tuple, list, dict, slice, none, capsule, iterable, iterator, function, buffer, array, array_t。C++到Python的转换，包括STL 数据结构、智能指针、类、操作符&函数重载、实例方法等，并且实现了对Eigen中的矩阵和向量的接口支持，对数组的操作十分方便。
Windows环境配置
查看网上windows系统下的环境配置的文档，编译环节主要在visual studio 2019上，环境配置过程比较复杂，这里贴出来作为记录参考。运行环境如下
Win10，Microsoft Visaul Studio 2019 x64
Miniconda3 , with python 3.7 base
pybind11是 header-only的，因此不需要编译动态链接库，直接解压使用即可。
下载地址： pybind/pybind11 （原地址） pybind/pybind11 （镜像）
官方文档： Intro - pybind11 documentation
下载后解压到你通常存放代码库的地方，并新建环境变量PYBIND11指向该路径，如：C:\VS_Lib\pybind11，内部文件夹结构如下
实际上只用到了include文件夹，test文件夹下是pybind11各模块测试的cpp和py代码，可以作为学习模块方法的参考。
系统Path环境变量添加（默认安装Miniconda3时自动添加）：
C:\ProgramData\Miniconda3
C:\ProgramData\Miniconda3\Library\mingw-w64\bin
C:\ProgramData\Miniconda3\Library\usr\bin
C:\ProgramData\Miniconda3\Library\bin
C:\ProgramData\Miniconda3\Scripts
C:\ProgramData\Anaconda3\Lib\site-packages
C:\ProgramData\Anaconda3\Lib\site-packages\numpy\core\include
如果需要在C++调用Python解释器执行代码(py::scoped_interpreter guard{};)，在Windows环境中需要添加额外的两个系统环境变量让pybind11能够找到解释器链接：
PYTHONHOME：C:\ProgramData\Miniconda3
PYTHONPATH：C:\ProgramData\Miniconda3\Lib\site-packages; C:\ProgramData\Miniconda3\DLLs; C:\ProgramData\Miniconda3\Lib
如果缺少这两个环境变量，会出现运行时错误，这是windows下的一个bug，相关讨论见此贴： pybind11 python embed on Windows 10, Fatal Python error: initfsencoding: unable to load the file system codec #1930
Fatal Python error: initfsencoding: unable to load the file system codec
接下来就是配置Visaul Studio项目属性，调用pybind11了。使用VS2019新建一个VC++项目和x64解决方案，添加新建属性页pybind11.props，方便以后复用。
1）设置编译输出类型
通用属性--常规--常规--目标文件扩展名：.pyd
通用属性--常规--项目默认值-配置类型：动态库.dll
2）添加include包含：
通用属性--VC++目录--常规--包含目录：
$(PYBIND11)\include
C:\ProgramData\Miniconda3\include
3）链接器配置：
链接器-常规-附加库目录：C:\ProgramData\Miniconda3\libs
链接器-输入-附加依赖项：python3.lib; python37.lib
4）其他配置(可选)：
VC++目录--常规--包含目录：$(EIGEN3) （使用Eigen3，环境变量EIGEN3需要事先配置）
C/C++-语言-符合模式：否 (/permissive)
C/C++-语言-OpenMP支持：是 (/openmp) （使用OPENMP）
编译文件，在输出目录下生成$(ProjectName).pyd文件，拷贝到python运行同级目录下，即可通过 import 导入模块。
pybind11绑定代码测试
通过编写PYBIND11_MODULE(ModuleName,m)模块，可以绑定C++编写的函数，详细的入门教程及语法可参考官方文档，这里我们简单演示一些基本功能，如
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <omp.h>
#include <iostream>
namespace py = pybind11;
using namespace py::literals;
// define func
int add(int i = 1, int j = 2) {
return i + j;
py::array_t<double> add_c(py::array_t<double> arr1, py::array_t<double> arr2) {
py::buffer_info buf1 = arr1.request(), buf2 = arr2.request();
if (buf1.shape != buf2.shape)
throw std::runtime_error("Input shapes must match");
/* No pointer is passed, so NumPy will allocate the buffer */
auto result = py::array_t<double>(buf1);
py::buffer_info buf3 = result.request();
double* ptr1 = (double*)buf1.ptr,
* ptr2 = (double*)buf2.ptr,
* ptr3 = (double*)buf3.ptr;
#pragma omp parallel for
for (ssize_t idx = 0; idx < buf1.size; idx++)
ptr3[idx] = ptr1[idx] + ptr2[idx];
return result;
PYBIND11_MODULE(example1, m) {
m.doc() = "pybind11 example-1 plugin"; // optional module docstring
m.def("add", &add, "A function which adds two numbers",
py::arg("i") = 1, py::arg("j") = 2);
m.def("add_c", &add_c, "A function which adds two arrays with c type");
这里m.doc即是python中要显示的模块信息，通过m.def绑定了两个函数“add”和“add_c”，分别实现整数相加与数组相加。
编译得到example1.pyd文件，注意文件名与PYBIND11_MODULE(ModuleName,m)中的ModuleName要一致，否则无法导入模块。把.pyd文件放在运行目录下启动python，接下来在python中使用如下
In [1]: import example1
In [2]: dir(example1)
Out[2]:
['__doc__',
'__file__',
'__loader__',
'__name__',
'__package__',
'__spec__',
'add',
'add_c']
In [3]: example1.add(10,20)
Out[3]: 30
In [4]: import numpy as np
In [5]: example1.add_c(np.ones((3,4,5)),np.ones((3,4,5)))
Out[5]:
array([[[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.]],
[[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.]],
[[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.],
[2., 2., 2., 2., 2.]]])
绑定Numpy数组
操作Numpy数组，需要引入头文件<pybind11/numpy.h>，通过pybind11::array_t<T>类可以接收numpy.ndarray数组。数组本质上在底层是一块一维的连续内存区，通过pybind11中的request()函数可以把数组解析成py::buffer_info结构体，buffer_info类型可以公开一个缓冲区视图，它提供对内部数据的快速直接访问。
struct buffer_info {
void *ptr; //指向数组(缓冲区)数据的指针
py::ssize_t itemsize; //数组元素总数
std::string format; //数组元素格式(python表示的类型)
py::ssize_t ndim; //数组维度信息
std::vector<py::ssize_t> shape; //数组形状
std::vector<py::ssize_t> strides;//每个维度相邻元素的间隔(字节数表示)
别的类型都比较直观清楚，需要特别注意strides这个元素代表的是每个维度相邻元素的字节间隔，numpy和C++的数组默认都是行优先存储的，那么对于一个row行cols列的二维float数组(矩阵)来说，每行“相邻”元素（如array(0,N)和array(1,N)）在行优先的内存上就相隔了一行元素大小的数量，即 sizeof (float) * cols，每列“相邻”元素（如array(N,0)和array(N,1)）在内存中是紧挨着的，即元素开头的位置相隔了一个元素类型的字节大小 sizeof (float)，那么这个数组和缓冲区的strides就是{ sizeof (float) * array.cols, sizeof (float) }，这在解析和重包装数组对象时经常用到。
py::array_t<double>add_c(py::array_t<double> arr1, py::array_t<double> arr2) 重新实现数组相加，arr.request()解析获得buffer_info对象，通过py::array_t<double>(buffer_info)传入整个buffer_info可以重新构造一个完全相同的数组对象，即实现深拷贝，也可以传入一个buffer_info.shape构造形状相同的新数组，都是开辟了新的内存空间。获取缓冲对象元素的指针buffer.ptr就可以操作元素完成运算，这里buffer.size是元素的总数，不管数组是多少维度的，其底层表示均是一维数组，可以用一个大循环直接遍历所有元素，实现数组元素相加。
double* ptr1 = (double*)buf1.ptr,
* ptr2 = (double*)buf2.ptr,
* ptr3 = (double*)buf3.ptr;
#pragma omp parallel for
for (ssize_t idx = 0; idx < buf1.size; idx++)
ptr3[idx] = ptr1[idx] + ptr2[idx];
这里采用openMP进行简单并行加速，实测运算速度可以达到numpy的90%以上。
Eigen数组接口
py::array_t的类函数通常比较有限（仅有访问元素\dims\shape\flag\dtype等基础功能），对标numpy中丰富的线性代数操作难以满足，通常需要转换成C++的数组和矩阵类型来完成相应功能，其中最广泛使用的的矩阵库就是Eigen，pybind11也实现了对Eigen一维和二维数组的直接转换支持，需要包含头文件<pybind11/eigen.h>。
下面是一个简单例子，函数参数和返回值都可以直接使用Eigen::Matrix<T>和Eigen::Array<T>的类型，pybind11会自动帮你转换。
using namespace Eigen;
MatrixXd add_mat(MatrixXd matA, MatrixXd matB)
return matA + matB;
namespace py = pybind11;
PYBIND11_MODULE(add_mat_moudle, m)
m.doc() = "Matrix add";
m.def("mat_add_py", &add_mat);
不过这样使用会发现函数运行速度比直接传入py::array_t<T>要慢很多，尤其是大型矩阵。实际上，当我们普通的Eigen::Matrix对象作为参数和返回值时，为了保证内存安全，pybind11接受 numpy.ndarray 的输入值，将其值复制到适当的临时数组变量，然后用这个临时变量调用C++函数。即 默认情况下是多进行了一次数组内存拷贝 的，对于计算量很小的矩阵四则运算等操作，这会显著增加函数运行的总时间！
那么有没有不进行复制而是引用传递的方法呢？答案是 使用 Eigen::Ref<MatrixType> 和 Eigen::Map<MatrixType> ，这样pybind11 默认会简单地引用返回的数据，但是须注意确保这些数据保持有效(不能在返回前被销毁)。特别注意，由于 numpy 和 Eigen 对数据的默认存储顺序不同(Numpy行优先，Eigen列优先) ，有时会遇到限制，需要在创建Eigen::Matrix对象使用Eigen::Rowmajor参数指定为行优先数组，否则转换时有可能会发生内存泄漏导致函数崩溃。
如果自定义函数中没有使用Eigen::Ref 和 Eigen::Map接收和返回参数，为了避免数组被复制，可以在绑定函数中使用pybind11的返回值策略 py::return_value_policy::reference_internal 来返回引用值
PYBIND11_MODULE(add_mat_moudle, m)
m.doc() = "Matrix add";
m.def("mat_add_ref", &add_mat, py::return_value_policy::reference_internal);
稀疏矩阵类型 scipy.sparse.csr_matrix/scipy.sparse.csc_matrix 不支持按引用传递，它们总是被复制按值传递的。
高维数组
无论是Eigen::Matrix还是Eigen::Array都是2D数组，其他C++矩阵库如numcpp也不支持多维数组类型，想要在C++中延续Numpy高维数组(ndim ≥ 3)的诸多特性并不容易。Eigen其实在附加模块中实现了一个全新的数组类型Eigen::Tensor，可以支持多维数组的部分运算和功能，但是到目前为止并没有被官方加入到加入到Eigen的主干当中，而是放在<unsupported/Eigen/CXX11/Tensor>模块下。其内置的方法和Eigen::Matrix类型差异很大，也没有直接好用的方法在Tensor和Matrix之间相互转换，官方文档也没有Eigen::Tensor的介绍（坑），可以从以下博客中快速学习Eigen::Tensor的常用特性：
在pybind11中需没有直接转换Eigen::Tensor类型的接口，需要通过py::buffer_info进行重构转换，下面是一个3维数组Tensor的例子
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <unsupported/Eigen/CXX11/Tensor>
namespace py = pybind11;
// return 3 dimensional ndarray with reversed order of input ndarray
template<class T>
py::array_t<T> eigenTensor(py::array_t<T> inArray) {
// request a buffer descriptor from Python
py::buffer_info buffer_info = inArray.request();
// extract data an shape of input array
T *data = static_cast<T *>(buffer_info.ptr);
std::vector<ssize_t> shape = buffer_info.shape;
// wrap ndarray in Eigen::Map:
// the second template argument is the rank of the tensor and has to be known at compile time
Eigen::TensorMap<Eigen::Tensor<T, 3>> in_tensor(data, shape[0], shape[1], shape[2]);
// Operate Eigen::Tensor Here
// build result tensor with reverse ordering
Eigen::Tensor<T, 3> out_tensor(2, 2, 2);
for (int i=0; i < shape[0]; i++)
for (int j=0; j < shape[1]; j++)
for (int k=0; k < shape[2]; k++)
out_tensor(k, j, i) = in_tensor(i, j, k);
// return numpy array wrapping eigen tensor's pointer
return py::array_t<T>(shape, // shape
{shape[0] * shape[1] * sizeof(T), shape[1] * sizeof(T), sizeof(T)}, // strides
out_tensor.data()); // data pointer
PYBIND11_MODULE(eigen, m) {
m.def("eigenTensor", &eigenTensor<double>, py::return_value_policy::move,py::arg("inArray"));
这里同样可以使用 Eigen::TensorMap 来包装得到Tensor数组的引用视图，从而避免接收参数时的复制操作，之后对 in_tensor 的操作都会等价映射到 inArray 上。同样返回值也必须转换成py::array_t<T>类型，通过构造函数py::array_t<T>(shape, strides , data pointer)来返回Tensor的引用视图。
直接使用Numpy和Python功能
通过Eigen::Tensor可以部分处理高维数组数据，但是Eigen::Tensor的特征方法较少且文档不健全，很多Numpy功能没有对应的函数，实际使用来说还是很不方便。为了保证Numpy的完整数组功能，我们希望可以在C++中直接使用Numpy的函数功能。
通过py::moudle::attr()就可以实现，py::moudle::attr()可以链接到当前激活的python环境，直接调用python中相应类型的函数，需要 #include <pybind11/embed.h>。py::module::import()可以将python标准库或当前python环境中定义的对象到C++环境下共同使用，这真正意义上的“混合编程”。
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <pybind11/embed.h>
#include <omp.h>
#include <iostream>
namespace py = pybind11;
using namespace py::literals;
py::object np = py::module::import("numpy");
void showarray(py::array_t<double> arr1)
auto local = py::dict();
py::object amax = arr1.attr("max")();
py::object shape = arr1.attr("shape");
py::array_t<double> a1 = np.attr("ones")(std::make_tuple(3, 4, 5),"dtype"_a="double");
py::print(a1);
local["arr1"] = arr1;
auto a2 = py::eval("(arr1==5)",local);
py::tuple a3 = np.attr("where")(a2);
int index = a3[0].cast<py::array_t<int>>().at(0);
py::print(a2);
py::print(index);
py::print(shape);
py::object np = py::module::import("numpy") 等价于python中的 import numpy as np，之后就可以使用np的函数和属性了，使用格式为 (py::object变量).attr("python中的函数名")(参数)。
py::object amax = arr1.attr("max")();
py::object shape = arr1.attr("shape");
py::array_t<double> a1 = np.attr("ones")(std::make_tuple(3, 4, 5),"dtype"_a="double");
等价于python中的
arr1.max()
arr1.shape
a1 = np.ones((3,4,5),dtype=double)
_a是py::literals中的迭代器别名，用来输入python中的Keyword参数，如"dtype"_a = "double"。
仍然有一些python特性无法通过py::moudle::attr()和C++方法给出（attr()中只能输入函数名，非函数特性则不行），比如数组的切片和列表索引特性、布尔数组等。pybind11 提供
eval
,
exec
和
eval_file
函数来直接运行 Python 表达式和语句，如下所示。
// At beginning of file
#include <pybind11/eval.h>
auto local = py::dict();
local["arr1"] = arr1;
auto a2 = py::eval("(arr1==5)",local);
// Evaluate a sequence of statements
py::exec(
"print('Hello')\n"