🎃

pybind11を使ってC++/FortranプログラムをPythonで動かす <その3>

に公開
C++
Fortran
pybind11
tech

はじめに

Python ➔ C++ ➔ Fortranの順にバインディングを行うことで、PythonからFortranプログラムを動かすテストの第3回です。入力パラメータがPython ➔ C++ ➔ Fortranの順に流れ、計算結果がFortran ➔ C++ ➔ Pythonの順に返ってくるような処理をイメージしています。

プログラム例 4

プログラム例3では関数の引数はコピーセマンティクスでしたが、NumPyのバッファプロトコルを使えばコピー処理を無くして、関数内から直接データにアクセスするようにできます。

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

// Using Fortran code within C++
extern "C" {
    void statistics(const int n, const double *x, double *mean, double *var);
}

namespace py = pybind11;

// Mean and variance
std::tuple<double, double> c_statistics_pyarray(const py::array_t<double> &x) {
    // Protocol buffer
    py::buffer_info buf = x.request();

    const int n = buf.shape[0];
    const double *x_ptr = static_cast<const double *>(buf.ptr);
    double mean, var;

    statistics(n, x_ptr, &mean, &var);

    return std::forward_as_tuple(mean, var);  // C++11
    // return {mean, var};                    // C++17
}

// Bindings with Python
PYBIND11_MODULE(example, m) {
    m.def("statistics", &c_statistics_pyarray);
}

ポイントは、

  • 関数に与える引数の型を py::array_t<double> にしたこと

にあります（pyはnamespaceで定義した名前です）。この型に実装されているrequestメソッドを使用すればプロトコルバッファ（buffer_info型）を取得することができ、これによって配列データの形状や先頭データを指すポインタなどの情報を引けるようになります。なお、pybind11/numpy.hヘッダのインクルードが必要です。倍精度浮動小数点数以外のデータでもキャストが働くので動かすことができますが（NumPy配列に限らずリストデータも扱えます）、後で述べるようにムーブセマンティクスでなくなるので注意が必要です。

>>> import example
>>> import numpy as np
>>> x = [2, 4, 6, 8, 10]
>>> print('(mean, var) =', example.statistics(x))
(mean, var) = (6.0, 10.0)

プログラム例 5

もうひとつ、データに直接アクセスする方法があります。unchecked、または mutable_unchecked メソッドを使用します。

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

extern "C" {
    void statistics(const int n, const double *x, double *mean, double *var);
}

namespace py = pybind11;

std::tuple<double, double> c_statistics_pyarray_unchecked(const py::array_t<double> &x) {
    // Unchecked proxy object
    auto r = x.unchecked<1>();

    const int n = r.shape(0);
    const double *x_ptr = r.data(0);
    double mean, var;

    statistics(n, x_ptr, &mean, &var);

    return std::forward_as_tuple(mean, var);  // C++11
    // return {mean, var};                    // C++17
}

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

先ほどと違うのは、

  • 関数内で unchecked<N> メソッドを使用している

点です。Nは配列の次元を表します。こちらを使っても直接的なデータアクセスを実現できますが、メソッドの名前が示すように、uncheckedで得られる参照には配列の次元やインデックス範囲などのチェックが働きません。

プログラム例 6

mutable_unchecked メソッドを使用するとデータを書き換えることもできます。要素をインクリメントするプログラムでこれを確かめます。

Fortranソースコード

increment.f90
module increment_mod
    use, intrinsic :: iso_c_binding, only: c_int, c_double
    implicit none
contains
    ! Increment all element by 1
    subroutine increment(n, x) bind(C)
        integer(c_int), intent(in), value :: n
        integer(c_int), intent(inout) :: x(n)
        x(:) = x(:) + 1
    end subroutine
end module

C++ソースコード

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

extern "C" {
    void increment(const int n, int *x);
}

namespace py = pybind11;

void c_increment_unchecked(py::array_t<int> &x) {
    // Unchecked proxy object
    auto r = x.mutable_unchecked<1>();

    const int n = r.shape(0);
    int *x_ptr = r.mutable_data(0);

    increment(n, x_ptr);
}

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

結果。

>>> import example
>>> import numpy as np
>>> x = np.arange(1, 10, 2).astype(np.int32)
>>> print('x =', x)
x = [1 3 5 7 9]
>>> example.increment_unchecked(x)
>>> print('x =', x)
x = [ 2  4  6  8 10]

もちろん、バッファプロトコルによっても実現できます。

example.cpp
void c_increment(py::array_t<int> &x) {
    // Protocol buffer
    py::buffer_info buf = x.request();

    const int n = buf.shape[0];
    int *x_ptr = static_cast<int *>(buf.ptr);
    
    increment(n, x_ptr);
}

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

コピーセマンティクスとムーブセマンティクス

requestメソッドによるバッファプロトコルの利用（プログラム例4）も、uncheckedメソッドによるダイレクトアクセス（プログラム例5,6）も、データがNumPyの32ビット整数型か倍精度浮動小数点数型の配列でなければコピーセマンティクスになることに注意が必要です。

例えば次のようなプログラムでそれが確認できます。

example.cpp
py::array_t<int>& c_increment(py::array_t<int> &x) {
    // Protocol buffer
    py::buffer_info buf = x.request();
    const int n = buf.shape[0];
    int *x_ptr = static_cast<int *>(buf.ptr);
    increment(n, x_ptr);
    return x;
}

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

Pythonで32ビット整数型と64ビット整数型（64ビット計算機のデフォルトの整数型）を与えたときの結果を比較します。

>>> import numpy as np
>>> import example
>>> x_int32 = np.arange(1, 10, 2).astype(np.int32)
>>> x_int64 = np.arange(1, 10, 2)
>>> print(x_int32, example.increment(x_int32))
[ 2  4  6  8 10] [ 2  4  6  8 10]
>>> print(x_int64, example.increment(x_int64))
[1 3 5 7 9] [ 2  4  6  8 10]

結果が異なっています。これは、64ビット整数型を与えても32ビット整数型への変換が内部的に発生し、インクリメントは変換後のコピーに作用するからです。倍精度浮動小数点数型の変数に対して32ビット整数を与えても同様のことが起こります。なまじプログラムが動くだけに注意が必要です。

PYBIND11_MODULEマクロのモジュール定義に py::arg().noconvert() アノテーションを追加して型の変換を抑止することで、意図しない動作を防ぐことができます。こうしておくと、変換が必要な型を与えたときにはエラーで終了してくれるので安心です。

example.cpp
PYBIND11_MODULE(example, m) {
    m.def("increment_int", &c_increment_int, py::arg().noconvert());
}

>>> import numpy as np
>>> import example
>>> x_int32 = np.arange(1, 10, 2).astype(np.int32)
>>> x_int64 = np.arange(1, 10, 2)
>>> print(x_int32, example.test_increment_int_noconvert(x_int32))
[ 2  4  6  8 10] [ 2  4  6  8 10]
>>> print(x_int64, example.test_increment_int_noconvert(x_int64))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: test_increment_int_noconvert(): incompatible function arguments. The following argument types are supported:
    1. (arg0: numpy.typing.NDArray[numpy.int32]) -> numpy.typing.NDArray[numpy.int32]

Invoked with: array([1, 3, 5, 7, 9])

64ビット整数型（64ビット計算機のデフォルトの整数型）を与えたときには、32ビット整数型しか対応していないよとメッセージを出してエラー終了してくれます。

参考

https://pybind11.readthedocs.io/en/stable/

Discussion