C++ Advent Calender 10日目の記事です。

Boost 1.63でBoost.NumPyがBoost.Pythonにマージされます。
それに伴い以下の記述では不具合があるかもしれません。修正記事を書く予定です。

動機

Pythonまじ便利です。
シミュレーションとその結果の解析がメインとなるアカデミア(見習い)である私には、
対話的にデータの加工、解析、可視化がシームレスに実行できるIPython Notebookは必須です。
残念ならがC++だけで同等の機能を提供してくれる環境はありません(たぶん)。
CERNが作ってるROOTのclingはC++を対話的に実行できるそうですが情報が少なくて使ったことありません。
だれかC++で文芸的プログラミングできるIC++ Notebook作ってください（切実）。

しかしながらシミュレーション自体は数日から数週間実行するもので、
速度が必須でありPythonなんて遅い言語で書けません。
なのでシミュレーションはC++で書く事になります（Fortranなにそれおいしいの）。

一方でscikit-learnやPyMCのようなプロジェクトがPython用のフロントエンドを用意している所を見ると、
今後暫くは他の言語で実装してPythonから使うという形が一般的になっていくような気がします。
ということで、Pythonで使える関数をC++で書くための方法のうち、
Boost.NumPyを使う方法についてまとめます。

何故Boost.NumPy?

私は以前一度Boost.Pythonを使おうとしてコンバータが理解できずに挫折しましたが、
Boost.NumPyは非常に簡単に使えました(大事)。

簡単に他の選択肢を述べておくと、
今回紹介するBoost.NumPyはシンプルにPythonのnumpy.ndarrayのC++用ラッピングを構築する方法を取ります。
C++には既にEigenのような線形代数演算を実行するライブラリが(無数に)あるので、
numpy.ndarrayをそのままEigenの多次元ベクトルに変換する選択肢もあります。
もっと単純にBoost.Pythonのみを使用してC++のvectorなどをPythonのlistに変換し、
その後numpy.ndarrayに変換する選択肢もあります。
他にもSWIGやCython等沢山あります。

目的

前回の記事でインストールとコンパイルの仕方をまとめたので今回は実際に使う方法をまとめます。
通常、C++のコードのインターフェイス部分だけをPython側から使いたいの思いますので、
難しい事は避けてなるべく簡単に使う方法をまとめます。

Boost.Pythonに関する説明はしないので、各自で調べてください。
以下にBoost.Pythonの解説サイトを挙げます：

Tutorial for Boost.NumPy

以下はBoost.NumPyのチュートリアルを元に私が書いたものです。
表記の複雑性の排除のため以下の様にnamespaceを略記します。

namespace p = boost::python;
namespace np = boost::numpy;

`np::ndarray`の使い方(一次元の場合)

前置きが長かったので、さっさと動くコードを載せましょう：

mymod1.cpp

# include "boost/numpy.hpp"
# include <stdexcept>
# include <algorithm>

namespace p = boost::python;
namespace np = boost::numpy;

/* 2倍にする */
void mult_two(np::ndarray a) {
  int nd = a.get_nd();
  if (nd != 1)
    throw std::runtime_error("a must be 1-dimensional");
  size_t N = a.shape(0);
  if (a.get_dtype() != np::dtype::get_builtin<double>())
    throw std::runtime_error("a must be float64 array");
  double *p = reinterpret_cast<double *>(a.get_data());
  std::transform(p, p + N, p, [](double x) { return 2 * x; });
}

BOOST_PYTHON_MODULE(mymod1) {
  Py_Initialize();
  np::initialize();
  p::def("mult_two", mult_two);
}

これはそのままコンパイルできます。
一次元配列を2倍にするだけの簡単な関数ですが、重要な要素を含んでいます。

配列の次元はget_nd()を使う
配列の大きさはshape(n)を用いる
配列の要素の型は動的に決定されget_dtype()で取得する
生ポインタでデータにアクセスできる
C++のstd::runtime_errorがPython側のRuntimeErrorに変換される

メモリはnp::ndarray側で管理しているので特に意識する必要はありません。
残念ながら型が動的に決定されるので関数のオーバーロードができません。
実行時にif文で判定する必要があります。
例外を変換してくれるのは地味に便利です。

Pythonから使ってみましょう。

mymod1.py

# !/usr/bin/env python
# -*- coding: utf-8 -*-

import mymod1
import numpy as np

if __name__ == '__main__':
    a = np.array([1,2,3], dtype=np.float64)
    mymod1.mult_two(a)
    print(a)

    b = np.array([1,2,3], dtype=np.int64)
    mymod1.mult_two(b) # raise error
    print(b)

bの部分は上述の通り型がdoubleでなくlong longなのでエラーになります：

[ 2.  4.  6.]
Traceback (most recent call last):
  File "/path/of/mymod1.py", line 13, in <module>
    mymod1.mult_two(b)
RuntimeError: a must be float64 array

なおnumpy上で整数を使用する場合、特に指定しない限りint64(C++だとlong long)が使用されます。
intでないので注意しましょう。

オーバーロードができないなど不満は多少ありますが、
簡単にPythonで使える関数をC++で実装できました。

`np::ndarray`の使い方(多次元の場合)

多次元の場合も基本は同じですが、メモリの並び順に気をつける必要があります。
numpy.ndarray.strideの説明が分かり易いですが、以下で簡単に説明します。

2次元の配列を一つの連続したメモリ領域で管理したい事があります。

double *a_as1 = new double[N*M];
double **a_as2 = new double*[N];
for(int i=0;i<N;++i){
  a_as2[i] = &a_as1[i*M];
}

こうすればa_as2[i][j]で呼び出されるメモリはa_as1[i*M+j]と同じになります。
このi*M+jのような事をするための方法がndarray.strideです。
今の場合j方向に1進めるのにメモリ上では8Byte進みます(doubleは8Byte)。
一方i方向に1進めるのにメモリ上では8*MByte進みます((i+1)*M+j = i*M+j + M)。
この8, 8*Mをstride(ひとまたぎ、歩幅)と呼びます。
この考え方はより高次元になっても使えます。

これに気をつければあとは簡単です。

# include "boost/numpy.hpp"
# include <iostream>
# include <stdexcept>
# include <algorithm>

namespace p = boost::python;
namespace np = boost::numpy;

void print(np::ndarray a) {
  int nd = a.get_nd();
  if (nd != 2)
    throw std::runtime_error("a must be two-dimensional");
  if (a.get_dtype() != np::dtype::get_builtin<double>())
    throw std::runtime_error("a must be float64 array");

  auto shape = a.get_shape();
  auto strides = a.get_strides();

  std::cout << "i j val\n";
  for (int i = 0; i < shape[0]; ++i) {
    for (int j = 0; j < shape[1]; ++j) {
      std::cout << i << " " << j << " "
                << *reinterpret_cast<double *>(a.get_data() + i * strides[0] +
                                               j * strides[1]) << std::endl;
    }
  }
}

BOOST_PYTHON_MODULE(mymod2) {
  Py_Initialize();
  np::initialize();
  p::def("print", print);
}

strideを見て小さい順に回す方がメモリアクセスの点で高速ですが、
面倒なので省略しました。各自がんばってください。

最後に

PythonというかNumPyまじ便利です。
SciPyが標準的な数値計算ライブラリのインターフェイスを丁寧なドキュメント付きで用意してくれている点が大きいのでしょう。

以上でnp::ndarrayのデータにアクセスする方法をまとめたので、
C++で実装したアルゴリズムをPythonから使えるようになったと思います。
次はPyPIへの投稿や、scikitsへ投稿について調べたいと思います。