CUDA kernel return value shortcut

CUDA の kernel 関数からの戻り値を毎回 cudaMemcpy したくない時につかう

template<class T>
class CudaValue
{
private:
    T* ptr;

public:
    CudaValue()
    {
        cudaMalloc(&ptr, sizeof(T));
    }

    ~CudaValue()
    {
        cudaFree(ptr);
    }

    operator T()
    {
        T ret;
        cudaMemcpy(&ret, ptr, sizeof(T), cudaMemcpyDeviceToHost);
        return ret;
    }

    CudaValue& operator=(const T& rhs)
    {
        cudaMemcpy(ptr, &rhs, sizeof(T), cudaMemcpyHostToDevice);
        return this;
    }

    T* get() { return ptr; }

private:
    CudaValue(const CudaValue&);
    CudaValue operator=(const CudaValue&);    
};

python & CUDA Integration

This template would be greatly helpful to whoever want to integrate python and cpp

test.cu

#include <iostream>

namespace test { namespace cuda {

extern "C" void hello()
{
    std::cout << "Hello" << std::endl;
}

}}

test.cpp

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

using namespace boost::python;

namespace test { namespace cuda {
extern "C" void hello();
}}

BOOST_PYTHON_MODULE(test)
{
    numeric::array::set_module_and_type("numpy", "ndarray");
    test::cuda::hello();
    import_array();
}

and then

setup.py

import os
import sys
from distutils.core import setup, Extension
from subprocess import Popen

modules = []
modules.append(['test', 'test.cpp', 'test.cu'])

for module in modules:
    mod_name = module[0]
    mod_cu_source = module[2]
    proc = Popen(['nvcc', '-Xcompiler', '-fPIC', '-c', mod_cu_source, '-o', '%s.o' % mod_name])
    proc.wait()

extensions = []
for module in modules:
    mod_name = module[0]
    mod_source = module[1]
    ex_module = Extension(
        mod_name,
        sources = [mod_source],
        include_dirs = ['/usr/include/boost'],
        libraries = ['boost_python', 'stdc++'],
        library_dirs = ['.', '/usr/local/cuda-5.5/lib64'],
        extra_compile_args=['-fopenmp', '-O3', '-std=c++0x'],
        extra_link_args=['%s.o' % mod_name, '-lcudart']
    )
    extensions.append(ex_module)

setup(name = "tools", version = "1.0", ext_modules = extensions)

If you place 'cudart' in the libraries section instead of using '-lcudart' in extra_link_args, libcudart.so would not be loaded and the process would stop at 'import test'.
Type 'ldd test.so' and check the results.

This code works fine on my PC (Ubuntu 12.04 LTS 64bit & CUDA 5.5 & Python 2.7.3).

Notes:
If you would omit '-Xcompiler -fPIC', you would get

/usr/bin/ld: test.o: relocation R_X86_64_32 against `.rodata' can not be used when making a shared object; recompile with -fPIC
test.o: could not read symbols: Bad value

and if you would place 'cudart' instead of '-lcudart', you would get

>>> import test
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ImportError: ./test.so: undefined symbol: __cudaUnregisterFatBinary

For more information, please see:
https://devtalk.nvidia.com/default/topic/468304/linking-c-and-cuda-files-with-nvcc-and-gcc/
http://daily.belltail.jp/?p=791

boost-python ではじめる大規模機械学習（６）- boost::multi_array のエクスポート・前編

あらすじ

boost-python を使用して、Python と C 言語両方を活用する方法を説明しています。

本記事では、boost::multi_array を Python 側にエクスポートする方法を解説します。

こんなものがつくりたい

Python でデータを読み込み、C で計算を行い、Python で途中経過を出力し、Python で結果を出力するプログラムを目標とします。
これは、たとえば次のようなコードになります。

import mylib, wrap
from matplotlib import pyplot as pl
from scipy import array, float64

# データを読み込む
x = array(..., dtype = float64)

# boost::multi_array 配列を確保
# x は入力、y は出力とします
wrap_x = wrap.make_vector(len(x))
wrap_y = wrap.make_vector(5)

# データを C 言語側に転送
wrap.copy_vector(x, wrap_x)

# (Python では不可能なほど大変な処理を) 実行します
for i in xrange(100):
    # 処理をおこなう
    mylib.execute(wrap_x, wrap_y)

    # 途中結果を表示
    y = wrap.scipy_vector(wrap_y)
    pl.plot(y)

# 結果をとりだす
y = wrap.scipy_vector(wrap_y)

mylib の中身は C 言語なので並列化などの最適化が容易です。
しかも、mylib の execute 関数にあたえる引数は boost::multi_array なので、実装もとても簡単です。
はじめに、wrap モジュールの概要を説明し、これを実装する方法を順番に紹介していきます。

おもな機能

wrap モジュールは大きく分けて、make (配列の確保)、reset (配列の初期化)、copy (SciPy -> C のデータ転送)、scipy (C -> SciPy のデータ転送) の 4 種類の機能を持ちます。
各関数はさらに、データ型と配列タイプをサフィックスとして持ちます。たとえば、reset_bool_matrix であれば、boolean 型の 2 次元配列を初期化する関数です。
さらに簡単のため、データ型が double の場合のみ型を省略することにします。
今回・次回合わせて、以下の関数の実装方法を説明します。

• make_vector
• make_matrix
• make_bool_vector
• make_bool_matrix
• reset_vector
• reset_matrix
• reset_bool_vector
• reset_bool_matrix

お急ぎの方は、ページ最後のソースコードをコピペして貼り付けるでも OK です。

共通モジュール

C 言語側でひろく利用する機能を、共通モジュールとして作り込みます。
ソースコードの名前は cmath.h とします。

インクルード、using 宣言

#include <boost/multi_array.hpp>
#include <boost/shared_ptr.hpp>
using boost::extents;
using boost::multi_array;
using boost::multi_array_ref;
using boost::shared_ptr;

型宣言

以下の型は、C 言語モジュールで繰り返し使うので、短縮名を与えておきます。
Python とのやりとりには multi_array に対するスマートポインタを使用するため、こちらも同時に定義しておきます。

typedef multi_array<double, 1> _wrap_vector;
typedef multi_array<double, 2> _wrap_matrix;
typedef multi_array<bool, 1> _wrap_bool_vector;
typedef multi_array<bool, 2> _wrap_bool_matrix;

typedef shared_ptr<_wrap_vector> wrap_vector;
typedef shared_ptr<_wrap_matrix> wrap_matrix;
typedef shared_ptr<_wrap_bool_vector> wrap_bool_vector;
typedef shared_ptr<_wrap_bool_matrix> wrap_bool_matrix;

dim, rdim

multi_array, multi_array_ref の次元サイズを取得するための補助関数を定義します。

template<class T>
static inline size_t dim1(T x)
{
    return x->size();
}

template<class T>
static inline size_t dim2(T x)
{
    return (*x)[0].size();
}

template<class T>
static inline size_t rdim1(T x)
{
    return x.size();
}

template <class T>
static inline size_t rdim2(T x)
{
    return x[0].size();
}

rnelem

multi_array_ref が何要素のデータ領域を指しているかを取得する補助関数を定義します。

template<class T>
static inline size_t rnelem(multi_array_ref<T, 1> x)
{
    return rdim1(x);
}

template<class T>
static inline size_t rnelem(multi_array_ref<T, 2> x)
{
    return rdim1(x) * rdim2(x);
}

ref 関数

C 言語側で、スマートポインタをはがすための関数を実装します。

template<class T>
multi_array_ref<T, 1> ref(shared_ptr<multi_array<T, 1> > x)
{
    multi_array_ref<T, 1> r(x->data(), extents[dim1(x)]);
    return r;
}

template<class T>
multi_array_ref<T, 2> ref(shared_ptr<multi_array<T, 2> > x)
{
    multi_array_ref<T, 2> r(x->data(), extents[dim1(x)][dim2(x)]);
    return r;
}

ソースコード

今回はここまで。
cmath.h はこれで完成です。次回は wrap モジュールを組み立てます。

cmath.h

#include <boost/multi_array.hpp>
#include <boost/shared_ptr.hpp>

#ifndef PYLIB_CMATH_INCLUDED
#define PYLIB_CMATH_INCLUDED
namespace cmath {
using boost::extents;
using boost::multi_array;
using boost::multi_array_ref;
using boost::shared_ptr;

typedef multi_array<double, 1> _wrap_vector;
typedef multi_array<double, 2> _wrap_matrix;
typedef multi_array<bool, 1> _wrap_bool_vector;
typedef multi_array<bool, 2> _wrap_bool_matrix;

typedef shared_ptr<_wrap_vector> wrap_vector;
typedef shared_ptr<_wrap_matrix> wrap_matrix;
typedef shared_ptr<_wrap_bool_vector> wrap_bool_vector;
typedef shared_ptr<_wrap_bool_matrix> wrap_bool_matrix;

// ========== ========== ========== ==========
//
//          dim
//
// ========== ========== ========== ==========
template<class T>
static inline size_t dim1(T x)
{
    return x->size();
}

template<class T>
static inline size_t dim2(T x)
{
    return (*x)[0].size();
}

// ========== ========== ========== ==========
//
//          rdim
//
// ========== ========== ========== ==========
template<class T>
static inline size_t rdim1(T x)
{
    return x.size();
}

template <class T>
static inline size_t rdim2(T x)
{
    return x[0].size();
}

// ========== ========== ========== ==========
//
//          rnelem
//
// ========== ========== ========== ==========
template<class T>
static inline size_t rnelem(multi_array_ref<T, 1> x)
{
    return rdim1(x);
}

template<class T>
static inline size_t rnelem(multi_array_ref<T, 2> x)
{
    return rdim1(x) * rdim2(x);
}

// ========== ========== ========== ==========
//
//          ref
//
// ========== ========== ========== ==========
template<class T>
multi_array_ref<T, 1> ref(shared_ptr<multi_array<T, 1> > x)
{
    multi_array_ref<T, 1> r(x->data(), extents[dim1(x)]);
    return r;
}

template<class T>
multi_array_ref<T, 2> ref(shared_ptr<multi_array<T, 2> > x)
{
    multi_array_ref<T, 2> r(x->data(), extents[dim1(x)][dim2(x)]);
    return r;
}
}
#endif

boost-python ではじめる大規模機械学習（５）- 小休止

あらすじ

boost-python を使用して、Python と C 言語両方を活用する方法を説明しています。
前記事では、SciPy の配列に C 言語からアクセスする方法を解説しました。
今回は、実装の方向性について考えます。

PyArrayObject は面倒

前回見た通り、SciPy の配列にアクセスするには何段もの手順を踏む必要があり、面倒です。
PyUblas を経由して BLAS を使用することもできますが、BLAS の最大の問題は 3 次元以上のテンソルをサポートしないことです。

C 言語の配列を Python にエクスポートする

問題解決のひとつの方法は、C 言語の配列オブジェクトを Python にエクスポートすることです。
多次元 vector などのコンテナを Python 側から読み取る、という方法です。
この方法の欠点は、C 言語で定義されたコンテナのラッパーを一通り定義する必要があることと、SciPy の便利関数の恩恵を得られないことです。

boost::multi_array を使おう

本記事では、上記 2 つのアプローチとは異なる方法、Python と C 言語それぞれで配列を確保し、必要に応じて相互にコピーするという戦略を取ります。
メモリ使用量が 2 倍になるため、当然処理は非効率です。しかし、この実装は純 Python の実装よりは十分に速く、かつ C 言語、Python それぞれで自然なコードが記述できるため、プロトタイピングという Python の長所を最大限に活かしているとも言えます。

次回からは、この方針のもと、機械学習アルゴリズム実装の一連の流れを紹介していきます。

boost-python ではじめる大規模機械学習（４）- 配列アクセス

あらすじ

boost-python を使用して、Python と C 言語両方を活用する方法を説明しています。
前記事では、任意の Python モジュールを C 言語から呼び出す方法を解説しました。
本記事では、SciPy の配列に C 言語からアクセスする方法を解説します。

配列オブジェクトの取得

色々な方法がありますが、メモリに直接アクセスするのが最も高速と考えられます。
はじめに、PyObject* object.ptr() を使用して PyObject のポインタを取得します。
次に、PyArray_FromObject で配列の構造体へのポインタを得ます。

static inline PyArrayObject* py_array(PyObject* ptr, int ndim)
{
    return (PyArrayObject*)PyArray_FromObject(ptr, NPY_FLOAT64, ndim, ndim);
}

void pick_double(object src, int index)
{
    PyArrayObject* psrc = py_array(src.ptr(), 1);
    Py_DECREF(psrc);
}

PyArray_FromObject の引数は、順にポインタ（PyObject*）、データ型（NPY_FLOAT64, NPY_INT32 など）、最小次元数、最大次元数です。最小次元数に 2、最大次元数に 3 を指定して 1 次元配列を入力するとエラーになる、という仕掛けです。
本シリーズでは、最小次元数・最大次元数には常に同じ値を指定することにします。

PyArray_FromObject 関数は Python 内部の参照カウンタの値を 1 増やすので、処理が終わったら Py_DECREF 関数を呼び出し、参照カウンタの値を減らす必要があります。

値の読み出し

SciPy の配列はメモリ上に連続に並んでいるとは限らないので、添字ごとのアドレスの移動量（ストライド）を取得します。
次に、char* PyArrayObject.data を使用して値を読み出せば完成です。

static inline int str1(PyArrayObject* p)
{
    return p->strides[0];
}

template<class T>
static inline T* valof(char* x, size_t i, int stride)
{
    return (T*)&x[i * stride];
}

double pick_double(object src, int index)
{
    PyArrayObject* psrc = py_array(src.ptr(), 1);
    char* pbuf = psrc->data;
    double ret = *valof<double>(pbuf, index, str1(psrc));
    Py_DECREF(psrc);
    return ret;
}

動作テスト

hello.cpp

#include <boost/python.hpp>
#include <numpy/arrayobject.h>
using namespace boost::python;

//
//  Helpers
//
static inline PyArrayObject* py_array(PyObject* ptr, int ndim)
{
    return (PyArrayObject*)PyArray_FromObject(ptr, NPY_FLOAT64, ndim, ndim);
}

static inline PyArrayObject* py_int_array(PyObject* ptr, int ndim)
{
    return (PyArrayObject*)PyArray_FromObject(ptr, NPY_INT32, ndim, ndim);
}

static inline int str1(PyArrayObject* p)
{
    return p->strides[0];
}

template<class T>
static inline T* valof(char* x, size_t i, int stride)
{
    return (T*)&x[i * stride];
}

//
//  Implementations
//
double pick_double(object src, int index)
{
    PyArrayObject* psrc = py_array(src.ptr(), 1);
    char* pbuf = psrc->data;
    double ret = *valof<double>(pbuf, index, str1(psrc));
    Py_DECREF(psrc);
    return ret;
}

int pick_int(object src, int index)
{
    PyArrayObject* psrc = py_int_array(src.ptr(), 1);
    char* pbuf = psrc->data;
    int ret = *valof<int>(pbuf, index, str1(psrc));
    Py_DECREF(psrc);
    return ret;
}

BOOST_PYTHON_MODULE(hello)
{
    numeric::array::set_module_and_type("numpy", "ndarray");
    def("pick_double", pick_double);
    def("pick_int", pick_int);
    import_array();
}

実行結果

>>> import hello
>>> from scipy import float64, int32, rand
>>> x = rand(10)
>>> x
array([ 0.38193373,  0.24639094,  0.02470218,  0.25281427,  0.41856296,
        0.52292722,  0.7646326 ,  0.04793535,  0.41331719,  0.88596331])
>>> y = int32(x * 10)
>>> y
array([3, 2, 0, 2, 4, 5, 7, 0, 4, 8], dtype=int32)
>>> hello.pick_double(x, 2)
0.02470218244761957
>>> hello.pick_int(y, 3)
2
>>>

boost-python ではじめる大規模機械学習（３）- Python モジュールの呼び出し

あらすじ

boost-python を使用して、Python と C 言語両方を活用する方法を説明しています。
前記事では、Python・C 言語間で簡単なオブジェクトを受け渡す方法を解説しました。
本記事では、任意の Python モジュールを C 言語から呼び出す方法を解説します。

準備

以降、色々と SciPy の機能を使うことがあります。このとき、以下のコードを拡張モジュール側に書き込む必要があります。

#include <numpy/arrayobject.h>

BOOST_PYTHON_MODULE(hello)
{
    numeric::array::set_module_and_type("numpy", "ndarray");
    import_array();
}

理由は SciPy のリファレンスに書いてあるハズ。

Python モジュールの読み込み

import 文を使用することで、Python モジュールを読み込むことができます。
各モジュールの属性は attr 関数で呼び出すことができます。
おどろくほどシンプルな動作です。

#include <boost/python.hpp>
#include <numpy/arrayobject.h>
using namespace boost::python;
static object scipy = import("scipy");
static object float64 = scipy.attr("float64");
static object uint8 = scipy.attr("uint8");
static object zeros = scipy.attr("zeros");

Python 関数の呼び出し

非常にシンプルです。迷うところがありません。

object get_float()
{
    return float64;
}

object make_float_zeros(int nx, int ny)
{
    return zeros(make_tuple(nx, ny), float64);
}

object make_byte_zeros(int nx, int ny)
{
    return zeros(make_tuple(nx, ny), uint8);
}

動作テスト

hello.cpp

#include <boost/python.hpp>
#include <numpy/arrayobject.h>
using namespace boost::python;
static object scipy = import("scipy");
static object float64 = scipy.attr("float64");
static object uint8 = scipy.attr("uint8");
static object zeros = scipy.attr("zeros");

object get_float()
{
    return float64;
}

object make_float_zeros(int nx, int ny)
{
    return zeros(make_tuple(nx, ny), float64);
}

object make_byte_zeros(int nx, int ny)
{
    return zeros(make_tuple(nx, ny), uint8);
}

BOOST_PYTHON_MODULE(hello)
{
    numeric::array::set_module_and_type("numpy", "ndarray");
    def("get_float", get_float);
    def("make_float_zeros", make_float_zeros);
    def("make_byte_zeros", make_byte_zeros);
    import_array();
}

実行結果

$ python
>>> import hello
>>> hello.get_float()
<type 'numpy.float64'>
>>> hello.make_float_zeros(2, 3)
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
>>> hello.make_byte_zeros(2, 4)
array([[0, 0, 0, 0],
       [0, 0, 0, 0]], dtype=uint8)
>>> 

次回は配列の内容にアクセスする方法を解説します。

動く隠れマルコフモデル（導出編・前編）- 動く PRML シリーズ（３）

やりたいこと

動く PRML シリーズ、第３回は隠れマルコフモデル (hidden Markov model, HMM) です。混合ガウス分布 (GMM)、変分混合ガウス分布 (VB-GMM) に続き第三回です。

反復繰り返し型の機械学習アルゴリズムを理解するには大きく分けて二つの方法があります。一つ目はもちろん、

• 更新式を導出すること。

反復アルゴリズムの理論的性質はすべて更新式の形に反映されています。従って、この更新式を自力で導出することはとても勉強になります。

そして、二つ目は、

• 各イテレーションの内容をグラフにプロットし、実際の挙動を確認すること。

本シリーズでは、更新式の導出・各イテレーションのプロット、の二本立てで機械学習アルゴリズムを深く理解することをめざします。
以下で早速、HMM の生成モデルについて解説していきたいと思います。

生成モデル

はじめに、隠れマルコフ系列の尤度を与えます。時系列のデータ数を N とし、各時刻でとりうる状態の数を K とします。観測データが時刻 n で k 番目の状態から生成されているとき、これを であらわします。
初期確率を 、状態 j から状態 k への遷移確率を であらわすと、 すべての同時分布は次のように書くことができます。

　　

観測変数の尤度は GMM とします。ここで重要なことは、隠れマルコフ系列・観測変数の尤度は独立に与えられる、ということです。このことを覚えておくと、後々の Forward-Backward アルゴリズムの実装が非常にすっきり見通しの良いものになります。

　　

HMM と EM アルゴリズム

さて…
GMM の E-Step では を計算しました。それは、この値が M-Step で必要になるからです。
HMM の場合は事態がやや複雑です。そこで、今回は M-Step を先に導出することにします。

M-Step

完全データの対数尤度関数は

　　
　　

と書けます。これを各変数で偏微分することで更新式を得ます。
以下、, と書きます。

初期確率の更新

　　

をもとにラグランジュの未定乗数法を使用します。すると

　　

を得ます。

遷移確率の更新

　　

をもとにラグランジュの未定乗数法を使用します。すると

　　

を得ます。

平均・分散（精度）の更新

GMM の場合と同じなので省略します。
基本的に、HMM の観測分布パラメータの更新式は、GMM などの単純な混合分布のパラメータ更新則と一致します。

動く隠れマルコフモデル（導出編２）執筆予定