pybind11とNumCppによるPythonとC++の連携

Pythonで書いたコードの速度が遅いときに一部の関数だけC++で記述したい場合があります。pybind11を使えばC++とPythonを比較的容易に連携させることができるのですが、NumPyオブジェクトをC++に渡したときの処理の記述が非常に面倒です。そこでNumPyのC++版として開発されているライブラリのNumCppを使って、NumPyオブジェクトの受け渡しを行う方法を検討しました。

TL; DR 現時点ではNumCppよりもEigenを使うほうが良さそうです

pybind11だけを使った場合

まずはpybind11だけでNumPyを使うとどうなるかについて説明します。pybind11ドキュメントのNumPyのページにある足し算を行う例を紹介します。

py::array_t<double> add_arrays(py::array_t<double> input1, py::array_t<double> input2) {
    py::buffer_info buf1 = input1.request(), buf2 = input2.request();

    /* shapeのチェックは省略 */

    auto result = py::array_t<double>(buf1.size);

    py::buffer_info buf3 = result.request();

    double *ptr1 = static_cast<double *>(buf1.ptr);
    double *ptr2 = static_cast<double *>(buf2.ptr);
    double *ptr3 = static_cast<double *>(buf3.ptr);

    for (size_t idx = 0; idx < buf1.shape[0]; idx++)
        ptr3[idx] = ptr1[idx] + ptr2[idx];

    return result;
}

行列の要素にアクセスするためにわざわざポインタを経由する必要がある上に、単なる行列の足し算をわざわざループで書き直す必要があります。この例は1次元配列でしたが、2次元以上になるとさらに処理が面倒になります。そこでC++の行列演算ライブラリを組み合わせて簡単に記述できるようにします。

pybind11は公式でEigenとの連携をサポートしているのですが、Eigenを使う方法は既に日本語解説が存在しているのでここではNumCppを使ってみます。

pybind11入門(4) NumPy連携その2 - Qiita

NumCppを使った場合

インストール

pybind11のインストールはpipを使う方法が一番簡単です。その他の方法はドキュメントのInstalling the libraryを参照してください。

pip install pybind11

NumCppはheader-onlyライブラリなのでファイルをダウンロードするだけです。ここでは作業ディレクトリにcloneすることにします。

git clone https://github.com/dpilger26/NumCpp.git

簡単な例

NumCppを使って二つの行列の和を求める関数を書いてみます。

以下のC++ファイルを example.cpp という名前で保存します。

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

#include <NumCpp.hpp>

namespace py=pybind11;

py::array_t<double> add_arrays(py::array_t<double, py::array::c_style> input1, py::array_t<double, py::array::c_style> input2) {
    auto array1 = nc::pybindInterface::pybind2nc(input1);
    auto array2 = nc::pybindInterface::pybind2nc(input2);
    auto sum = array1 + array2;
    return nc::pybindInterface::nc2pybind(sum);
}

PYBIND11_MODULE(example, m) {
    m.def("add_arrays", &add_arrays);
}

Linuxの場合コンパイルは以下のように行います。

c++ -O3 -Wall -shared -std=c++14 -fPIC $(python3 -m pybind11 --includes) -I NumCpp/include/ -DNUMCPP_INCLUDE_PYBIND_PYTHON_INTERFACE example.cpp -o example$(python3-config --extension-suffix)

これはpybind11のコンパイルコマンドに NumCppのヘッダファイルの場所指定 (-I NumCpp/include/) と NumCppのpybind11オプション (-DNUMCPP_INCLUDE_PYBIND_PYTHON_INTERFACE) を追加してC++14を指定 (-std=c++14) したものです。CMakeを使う方法などもありますが、複雑なのでここでは解説しません。

コンパイルに成功すると example.cpython-***.so のようなファイルが生成されます。これを読み込むPythonプログラムは以下のようになります。

import numpy as np
import example

a = np.asarray([1, 2, 3, 4]) 
b = np.asarray([2, 3, 4, 5])
print(example.add_arrays(a, b).shape) # [[3. 5. 7. 9.]]

正しく実行できているように見えるのですが、1次元配列同士の和なのに2次元配列が返ってきています。ドキュメントによると、NumCppのNdArrayクラスは2次元配列として実装されているため、1次元配列は1xNの2次元配列として扱う仕様になっているそうです。

2次元配列は問題なく処理できます。

import numpy as np
import example

a = np.asarray([[1, 2], [3, 4]]) 
b = np.asarray([[2, 3], [4, 5]])
print(example.add_arrays(a, b)) 
# [[3. 5.]
#  [7. 9.]]

しかし、3次元配列を渡すとエラーが起こります。

import numpy as np
import example

a = np.asarray([[[1], [2]], [[3], [4]]]) 
b = np.asarray([[[2], [3]], [[4], [5]]]) 
print(example.add_arrays(a, b))

エラーメッセージ

File: NumCpp/include/NumCpp/PythonInterface/PybindInterface.hpp
    Function: pybind2nc
    Line: 88
    Error: input array must be no more than 2 dimensional.Traceback (most recent call last):
  File "test.py", line 6, in <module>
    print(example.add_arrays(a, b)) 
ValueError: File: NumCpp/include/NumCpp/PythonInterface/PybindInterface.hpp
    Function: pybind2nc
    Line: 88
    Error: input array must be no more than 2 dimensional.

多次元配列についてのissueを見る限り3次元以上の配列はまだ実装されてなさそうなので、現時点でNumPy互換と呼ぶのはちょっと難しいかもしれません。

おまけ: Eigenを使った場合

参考としてEigenを使った場合を示しておきます。MatrixをRowMajorにすることと、VectorとMatrixを区別することだけ気をつければそんなに難しくないので、よりメジャーなEigenを使うほうが良い場合が多いと思います。

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>

#include <Eigen/Dense>

namespace py = pybind11;

using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;

template <typename T>
T add_arrays(Eigen::Ref<const T> input1, Eigen::Ref<const T> input2) {
    return input1 + input2;
}

PYBIND11_MODULE(example, m)
{
    m.def("add_vectors", &add_arrays<Eigen::VectorXd>);
    m.def("add_matrices", &add_arrays<RowMatrixXd>);
}

Pythonはこんな感じです。

import numpy as np
import example

a = np.asarray([1, 2, 3, 4]) 
b = np.asarray([2, 3, 4, 5])
print(example.add_vectors(a, b))  # array([3., 5., 7., 9.])

a = np.asarray([[1, 2], [3, 4]])
b = np.asarray([[2, 3], [4, 5]])
print(example.add_matrices(a, b))
# array([[3., 5.],
#        [7., 9.]])