Boost.PythonとNumPy Array (C API編)

Boost.PythonでNumPyのarrayを受け取る手段として, NumPyが提供するCのAPIを使うという手がある. 使いづらいが機能的には完全であることと, CとPythonを協調して使う際の知見となるため, 記録を残しておく.

基本的な使い方はこちらのページを見ればわかる. Boost.PythonとNumPy C APIを同時に使っている例は, この投稿に載っている. 今回のサンプルはこの投稿につけられた「こうした方がいいよ」という点を反映したものである.

基本的な操作はBoost.Pythonのobjectクラスを使って行うが, 一部のNumPy APIを使う場合は, object.ptr()でPyObject型のポインタを取り出し, それをNumPy APIに渡してやるというのが基本方針である.

#include <Python.h>

#include <numpy/arrayobject.h>
#include <boost/python.hpp>
#include <boost/python/numeric.hpp>
#include <cmath>

using namespace std;
using namespace boost;
using namespace boost::python;

double trace(boost::python::object& arg)
{
  //配列かチェック
  if(!PyArray_Check(arg.ptr())){
    //例外が発生したことを示すフラグをセット
    PyErr_SetString(PyExc_ValueError,
                    "arg must be nd-array");
    //Boostに例外の発生を伝える
    throw_error_already_set();
    return 0;
  }

  //PyObjectからPyArrayObjectを生成
  PyArrayObject* array = (PyArrayObject*)PyArray_FROM_OTF(arg.ptr(),NPY_NOTYPE,NPY_IN_ARRAY);

  //スマートポインタで, arrayの寿命を管理
  python::handle<> parray((PyObject*)array);

  if (array->nd != 2 || array->descr->type_num != PyArray_DOUBLE) {
    PyErr_SetString(PyExc_ValueError,
                    "array must be two-dimensional and of type float");
    throw_error_already_set();
    return 0;
  }

  int rows = array->dimensions[0];
  int cols = array->dimensions[1];
  int n = min(rows,cols);
  
  double sum = 0.;
  for (int i = 0; i < n; i++)
    sum += *(double *)(array->data + i*array->strides[0] + i*array->strides[1]);
  
  return sum;
}

//BOOST_PYTHON_MODULE()の引数がmodule名
BOOST_PYTHON_MODULE(simple_array){
  import_array();
  numeric::array::set_module_and_type("numpy", "ndarray");
  
  def("trace",trace);
}

まずはarray特有の操作に絞って解説しよう.

if(!PyArray_Check(arg.ptr())){

これは, PyObject*がPyArrayに変換可能かをチェックする関数である.

PyArrayObject* array = (PyArrayObject*)PyArray_FROM_OTF(arg.ptr(),NPY_NOTYPE,NPY_IN_ARRAY);

こっちは, PyObjectからPyArrayObjectに変換する処理である.
PyArray_FROM_OTFはPyArray_FROM_ANYをラッピングしたマクロである. 第二引数はarrayのデータタイプ(doubleとかintとか), 第三引数はArrayが満たさなければいけない条件を渡す.
第二引数がNPY_NOTYPEの場合は, PyObject*を元に型が設定される. この例では第三引数は行列が連続したメモリに配置されることを強制している.
詳細はこちら.
行列の要素を足し合わせる部分は

  for (int i = 0; i < n; i++)
    sum += *(double *)(array->data + i*array->strides[0] + i*array->strides[1]);

である. わざわざstridesを介してアクセスしているのは連続したメモリに配置されているかは, 一般には成り立たないためである. stridesの値はデータタイプの個数では無く, バイトサイズであることに注意.

次に, これまで解説していないBoost部分の機能を説明する.

PyErr_SetString(PyExc_ValueError,
"arg must be nd-array");
//Boostに例外の発生を伝える
throw_error_already_set();
return 0;

ここはPythonのValueErrorを発生させる部分である. このページいわく, PythonAPIの仕組みとして, スレッドごとにglobal indicatorが存在し, 例外を発生させる際は, そのindicatorのフラグを立てた上で関数がNULLなり, 0なりを返せばいいらしい. PyErr_SetStringはindicatorのフラグを立てる関数の一つである.
ところがそのままではBoostがindicatorをクリアしてしまうらしく, throw_error_already_set()をコールしてそれを抑制しなければいけないっぽい.

さて, 残る説明は, オブジェクトの寿命管理に関わる

python::handle<> parray((PyObject*)array);

である. これは少しめんどくさい. しかも自信がないのでここから先間違ってたらごめんなさい.
Pythonは基本的に参照カウントでガベージコレクションを行なっている(参考ページ). よって, C++内でオブジェクトを生成したりした場合, 参照カウントを足したり引いたりしなければいけない.
例えば, おそらくPyArray_FROM_OTFは渡したオブジェクトが変換不要な場合は, そのオブジェクトの参照を返し, そうでない場合は新しい領域を確保した後にPyArrayObject*を生成する.
大事なのは, どちらにしろ参照カウントが1増加する点である. よって参照カウントを減らさなければ後々メモリリークが発生してしまう.

Boostに頼らない場合は, Py_DECREFを使って参照カウントを減らすのだが, 今回はBoost::python::handleを使って参照カウントを管理している. これは, スマートポインタの一種でたぶんスコープを外れた時点で参照カウントを1減らしてくれる.

このhandleにはいくつかバリエーションがある. このページに列挙されているが, NULLポインタを許さないhandleを生成することができる.

重要な概念としてborrowedなポインタとnewなポインタというものがある. NumPy C APIPython APIが返すPyObjectポインタはborrowedとnewに大別される.
newなポインタは, そのポインタを得る際に参照カウントがインクリメントされており, そのポインタが指すオブジェクトを使い終わったなら, Py_DECREFで参照カウントをデクリメントしなければいけない.
一方, borrowedなポインタは, そのポインタを得る際に参照カウントはインクリメントされていないため, Py_DECREFする必要は無い. しかし, そのポインタが指すオブジェクトが余所で参照カウントがデクリメントされた際に, まだ使っているにも関わらず解放されるリスクが存在する.

handleはborrowedなポインタを扱うか, それともnewなポインタを扱うかを,

handle<> y(borrowed(x))
or
handle<> y(x)

という様な記法で区別している. borrowed(x)と書いた場合のみ, この瞬間に予期せぬ解放を防ぐため参照カウントをインクリメントするという違いがある.

今回はnewなポインタを得ているため, borrowed関数は使用していない.

まとめ

まあこれでBoostの使い方と行列アクセスの方法がわかったが, めんどくさい.
なんらかのwrapperを使いたくなる所ではある.