pybind11入門(3) NumPy連携その１

その1 https://qiita.com/lucidfrontier45/items/872c38b5eec08e3800d7
その2 https://qiita.com/lucidfrontier45/items/c7cc409a3af962d100a3

pybind11を使ってNumPyと連携する方法を試した。まずは直接ndarrayを触る方法。

## `pybind11::array`型

pybind11には`pybind11::buffer`と`pybind11::array`のに種類のクラスが定義されており、これらを介してndarrayを扱うことができる。前者はndarrayに限らない一般的なPythonのbufferプロトコルと連携できるものであり、ここではndarray専用の`array`型を試してみる。なお、`array_t<T>`というtemplateも用意されており、実際にはこちらを用いる。いつもどおり、cmakeを用いる。

```cmake:CMakeLists.txt
cmake_minimum_required(VERSION 3.2)
project(pybind_test VERSION 0.1.0)
set(CMAKE_POSITION_INDEPENDENT_CODE ON)
add_subdirectory(pybind11)
add_subdirectory(fmt)
pybind11_add_module(mymodule module.cpp)
target_link_libraries(mymodule PRIVATE fmt::fmt)
```

`fmt`という書式フォーマットのライブラリも使用する。
https://github.com/fmtlib/fmt

## `ndarray`へのReadOnlyアクセス

２次元配列の各要素を表示するだけのプログラムで動作を確認してみる。

```c++:module.cpp
#include <iostream>
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <fmt/core.h>

namespace py = pybind11;

template <typename T>
void print_array(py::array_t<T> x)
{
    const auto &buff_info = x.request();
    const auto &shape = buff_info.shape;
    std::cout << "C++" << std::endl;
    for (auto i = 0; i < shape[0]; i++)
    {
        for (auto j = 0; j < shape[1]; j++)
        {
            auto v = *x.data(i, j);
            std::cout << fmt::format("x[{}, {}] = {}", i, j, v) << std::endl;
        }
    }
}

PYBIND11_MODULE(mymodule, m)
{
    m.doc() = "my test module";
    m.def("print_array", &print_array<int32_t>, "");
    m.def("print_array", &print_array<double>, "");
}
```

ndarrayを扱う場合には`pybind11/numpy.h`をincludeする。`request`メソッドで`buffer_info`型の変数を得ることがでる。`buffer_info`は以下のような構造になっており

```c++
struct buffer_info {
    void *ptr;
    ssize_t itemsize;
    std::string format;
    ssize_t ndim;
    std::vector<ssize_t> shape;
    std::vector<ssize_t> strides;
};
```

ndarrayでおなじみの`shape`などが得られるのでこれを用いでfor loopを回せます。データにアクセスする際には`x.data(i, j)`のようにするとその要素へのポインタが得られる。なお、`data`メソッドはReadonlyのアクセスなので得られるポインタにはconstがついている。

最後に、この関数はtemplateを使用しているのでpybind11に登録する際には実体化する必要がる。ここでは整数用の`int32_t`と実数用の`double`の２種類を作る。

Pythonで実行して確認する。

```python
import numpy as np
import mymodule

x = np.arange(2*3).reshape((2, 3)).astype(np.int32)
print("python")
for i in range(x.shape[0]):
    for j in range(x.shape[1]):
        print("x[{}, {}] = {}".format(i, j, x[i, j]))

mymodule.print_array(x)

>>>
python
x[0, 0] = 0
x[0, 1] = 1
x[0, 2] = 2
x[1, 0] = 3
x[1, 1] = 4
x[1, 2] = 5
C++
x[0, 0] = 0
x[0, 1] = 1
x[0, 2] = 2
x[1, 0] = 3
x[1, 1] = 4
x[1, 2] = 5
```

ちゃんと同じ内容が表示できた!

## 配列の変更

`.data(i, j)`の代わりに`.mutable_data(i, j)`を用いることで**constでない**ポインタが得られ、内容を書き換えることができる。以下のプログラムは全配列の要素に定数を足す例である。

```c++
template <typename T>
void modify_array_inplace(py::array_t<T> x, T a)
{
    const auto &buff_info = x.request();
    const auto &shape = buff_info.shape;
    for (auto i = 0; i < shape[0]; i++)
    {
        for (auto j = 0; j < shape[1]; j++)
        {
            *x.mutable_data(i, j) += a;
        }
    }
}

PYBIND11_MODULE(mymodule, m)
{
    m.doc() = "my test module";
    m.def("modify_array_inplace", &modify_array_inplace<int32_t>, "");
    m.def("modify_array_inplace", &modify_array_inplace<double>, "");
}
```

Pythonで試してみる。

```py
x = np.random.randn(2, 3).astype(np.float64)
print("before")
for i in range(x.shape[0]):
    for j in range(x.shape[1]):
        print("x[{}, {}] = {:.6f}".format(i, j, x[i, j]))

mymodule.modify_array_inplace(x, 10.0)
print("after")
for i in range(x.shape[0]):
    for j in range(x.shape[1]):
        print("x[{}, {}] = {:.6f}".format(i, j, x[i, j]))

>> 
before
x[0, 0] = 1.307738
x[0, 1] = -1.625780
x[0, 2] = 0.196201
x[1, 0] = 0.614617
x[1, 1] = -0.175023
x[1, 2] = 1.490781
after
x[0, 0] = 11.307738
x[0, 1] = 8.374220
x[0, 2] = 10.196201
x[1, 0] = 10.614617
x[1, 1] = 9.824977
x[1, 2] = 11.490781
```

配列の各要素に一律に10.0が足されたのが分かる。

なお、

```py
mymodule.modify_array_inplace(x, 10)
```

のようにすると`x`がdoubleの配列なのにintを足していて`modify_array_inplace<int32_t>`のほうが呼ばれたらしく、`x`は型が異なるのでコピーが渡されて配列の中身は変更されなかった。渡す引数がintなのか、floatなのかは注意したほうが良さそう。

## array_tを返す関数

`array_t`型の変数をを新規に作成し、戻り値にすることでPython側に新しいndarrayを返すことができる。

```c++
template <typename T>
auto modify_array(py::array_t<T> x, T a)
{   
    const auto &buff_info = x.request();
    const auto &shape = buff_info.shape;
    py::array_t<T> y{shape};
    for (auto i = 0; i < shape[0]; i++)
    {
        for (auto j = 0; j < shape[1]; j++)
        {
            *y.mutable_data(i, j) = *x.data(i, j) + a;
        }
    }
    return y;
}
```

`array_t`型のコンストラクタは色々あるが、`shape`を渡すものが最も簡単に使用できると思う。Pythonで利用する場合は以下の通り。

```py
x = np.arange(2*3).reshape((2, 3)).astype(np.int32)
x2 = mymodule.modify_array(x, 10)
```

## direct access

余計なチェックを無くしてよりオーバーヘッドが少ないアクセスの仕方もある。`@cython.boundscheck(False)`的なものだと思う[^unchecked]。

```c++
template <typename T>
void modify_array_direct(py::array_t<T> x, T a)
{
    auto r = x.template mutable_unchecked<2>();
    for (auto i = 0; i < r.shape(0); i++)
    {
        for (auto j = 0; j < r.shape(1); j++)
        {
            r(i, j) += a;
        }
    }
}
```

公式docでは`x.mutable_unchecked<2>();`となっているがこれだとコンパイルエラー。githubのissueに[同様のもの](https://github.com/pybind/pybind11/issues/1412)があって`x.template mutable_unchecked<2>();`とすれば良いらしい[^template]。ちなみにtemplateを使用せず、

```c++
void modify_array_direct(py::array_t<double> x, double a){
    auto r = x.mutable_unchecked<2>();
    ...
}
```
のような場合は問題なかった。

## 参考

いつもの公式doc: https://pybind11.readthedocs.io/en/stable/advanced/pycpp/numpy.html
uncheckedとtemplateを組み合わせる黒魔術: https://github.com/pybind/pybind11/issues/1412

[^unchecked]: `unchecked`を使用しない方法だとindexの範囲外にアクセスするとPythonの`IndexError`が出たが、`unchecked`の場合は特にエラーが出ずにそのまま上書きできた。おそろしやー。
[^template]: 自分のC++力が低すぎてなぜそれで良いのかよくわからない。。。