2012-06-15

pyopencl.arrayの使い方

Python OpenCL

PyOpenCLはPythonからOpenCLのAPIにアクセスするためのライブラリですが、それだけではなく素のOpenCLにはない機能もいくつか追加されています。pyopencl.arrayは、numpy.ndarrayライクのインタフェースでOpenCLのデバイスを意識することなく基本的な行列演算を行うことができます。しかしながら、PyOpenCLのドキュメントからは使い方がわかりにくい面もありますので、今回はその基本的な使い方について紹介します。numpyの基本的な知識を要求する気がします。

Arrayオブジェクトの作成

Arrayオブジェクトを作る方法は大きくわけて3つあります。

hostのデータをもとに作成する
device上で生成する
既にあるArrayオブジェクトから作成する

それぞれの方法についてソースコードを例示します。

import pyopencl as cl
from pyopencl import array as clarray
from pyopencl import clmath
from pyopencl import clrandom
import numpy

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

# 1. hostのデータをもとに作成する

host_a = numpy.random.rand(100)
a = clarray.to_device(queue, host_a)

# 2. device上で生成する

# 空の浮動小数点数型Arrayオブジェクト
b = clarray.empty(queue, 1000, numpy.float32)
# 0で初期化された2次元整数型Arrayオブジェクト
c = clarray.zeros(queue, (10, 10), numpy.int32)
# range
d = clarray.arange(queue, 0, 100, 10, dtype=numpy.int32)
# 乱数生成
e = clrandom.rand(queue, 100, numpy.float32)

# 3. 既にあるArrayオブジェクトをもとに作成する

f = clarray.empty_like(a)
g = clarray.zeros_like(a)

なお、emptyやempty_likeで作ったArrayオブジェクトはfillメソッドやsetメソッド、clrandom.fill_randを使って初期化することができます。

b.fill(2.5)

clrandom.fill_rand(f, a=-10, b=10)

h = clarray.empty_like(a)
h.set(host_a)

上記のように、randとfill_rand関数に引数a, bを与えると乱数の生成範囲が[a, b)になります。デフォルト値は[0, 1)です。

Arrayオブジェクトの演算

Arrayオブジェクトはnumpy.ndarrayと同様に演算をすることができます。

スカラー値との演算

a = clrandom.rand(queue, 5, numpy.int32, a=-10, b=10)
print a.get()

print a + 1
print a - 1.5
print a * 3
print a / 2.0
print a ** 2.0

getメソッドはdevice上にあるArrayオブジェクトの内容をhost側のndarrayオブジェクトにコピーします。Arrayオブジェクトの__str__メソッド、__repr__メソッドの中ではgetメソッドが使われており、printなどで表示するときには自動的にhost側に転送されるため、本記事では以降省略するものとします。

実行結果

[ 4 -8 -6 -8  2]  
[ 5 -7 -5 -7  3]  
[ 2.5 -9.5 -7.5 -9.5  0.5] 
[ 12 -24 -18 -24   6] 
[ 2. -4. -3. -4.  1.]
[ 16.  64.  36.  64.   4.]

Arrayオブジェクト同士の演算

b = clrandom.rand(queue, 5, numpy.int32, a=-10, b=10)
print b

print a + b
print a - b
print a * b
print a / b
print a ** b.astype(numpy.float32)

実行結果

[ 9  7 -6  2 -8] 
[ 13  -1 -12  -6  -6]  
[ -5 -15   0 -10  10]
[ 36 -56  36 -16 -16]
[ 0 -1  1 -4  0] 
[  2.62144000e+05  -2.09715200e+06   2.14334705e-05   6.40000000e+01
   3.90625000e-03]

単項演算

print -a
print abs(a)
print len(a)

実行結果

[-4  8  6  8 -2] 
[4 8 6 8 2]
5

numpy.ndarrayでは他にもビット演算などがサポートされていますが、PyOpenCLでは今のところサポートされていません。

数学関数

Arrayオブジェクトの各要素に数学関数などOpenCLでサポートされた関数を適用するためにはpyopencl.clmath以下の関数を使います。OpenCLの仕様上Arrayオブジェクトのdtypeは浮動小数点数型である必要があります。利用可能な関数についてはPyOpenCLのドキュメントやOpenCLの仕様書もご覧ください。

a = clrandom.rand(queue, 5, numpy.float32)
print a
     
print clmath.sqrt(a)
print clmath.sin(a)
print clmath.exp(a)

実行結果

[ 0.45427519  0.19955009  0.48750299  0.80684304  0.14642352]
[ 0.67399943  0.44671029  0.69821417  0.89824444  0.38265327]
[ 0.43881115  0.19822837  0.46842125  0.72210687  0.14590086]
[ 1.5750314   1.22085333  1.62824535  2.24082255  1.15768635]

Reduction

配列の総和を求めたり、最小値、最大値を求めるような処理はリダクションと呼ばれます。GPUでこのようなリダクション処理を効率よく行わせる場合には数々のチューニングテクニックが必要で、結構な量のコードを書く必要がありますが、PyOpenCLではこのようなショートカットが用意されているので詳細を気にせずに高速な処理を行わせることができます。

a = clrandom.rand(queue, 1000, numpy.float32)
b = clrandom.rand(queue, 1000, numpy.float32)

print clarray.sum(a)
print clarray.min(a)
print clarray.max(a)
print clarray.dot(a, b)

実行結果

510.130096436          
0.00353157520294
0.998835802078         
251.52911377

Subset

Arrayオブジェクトから指定したインデックスの要素を取り出した新しいArrayオブジェクトを作ったり、そのような部分配列からリダクションを行う機能が用意されています。

a = clrandom.rand(queue, 1000, numpy.float32)
b = clrandom.rand(queue, 1000, numpy.float32)
indices = clarray.arange(queue, 0, 1000, 2, dtype=numpy.int32)

c = clarray.take(a, indices)
print clarray.subset_dot(indices, a, b)
print clarray.subset_min(indices, a)
print clarray.subset_max(indices, a)

実行結果

120.567108154
0.00146424770355
0.999368607998

Conditionals

以下の関数は2つのArrayオブジェクトの各インデックスにおいて条件にあう要素を取り出します。if_positiveでは第1引数の要素が0以下の場合第3引数の要素が、そうでない場合第2引数の要素が選ばれます。

a = clrandom.rand(queue, 5, numpy.float32)
b = clrandom.rand(queue, 5, numpy.float32)
criterion = clrandom.rand(queue, 5, numpy.int32, b=2)
print a
print b
print criterion

print clarray.if_positive(criterion, a, b)
print clarray.maximum(a, b)
print clarray.minimum(a, b)

実行結果

[ 0.45082718  0.26803267  0.16043031  0.3467738   0.0948478 ]
[ 0.35381329  0.95943987  0.29650748  0.53794491  0.68245429]
[0 1 0 0 0]
[ 0.35381329  0.26803267  0.29650748  0.53794491  0.68245429]
[ 0.45082718  0.95943987  0.29650748  0.53794491  0.68245429]
[ 0.35381329  0.26803267  0.16043031  0.3467738   0.0948478 ]

カスタムKernel

先述のリダクションのように効率よく並列できるパターンというのがいくつかあります。PyOpenCLではそういったパターンの内Elementwise(transform)、Reduction、Scanについてメタプログラミングの仕組みが提供されています。

from pyopencl.elementwise import ElementwiseKernel
from pyopencl.reduction import ReductionKernel
from pyopencl.scan import ExclusiveScanKernel
from pyopencl.scan import InclusiveScanKernel
    
a = clrandom.rand(queue, 5, numpy.float32)
print a

k1 = ElementwiseKernel(ctx, 'float* a', 'a[i] = sin(a[i]);')
b = a * 1
k1(b)                        
print b
                             
k2 = ReductionKernel(ctx, numpy.float32, neutral='0', reduce_expr='a + b', map_expr='x[i]', arguments='__global float* x')
print k2(a)

k3 = ExclusiveScanKernel(ctx, numpy.float32, scan_expr='a + b', neutral='0')
print k3(a * 1)
           
k4 = InclusiveScanKernel(ctx, numpy.float32, scan_expr='a + b')
c = clarray.empty_like(a)
print k4(a, c)

実行結果

[ 0.48551953  0.27677017  0.63273305  0.30237836  0.28056026]
[ 0.46666792  0.27325016  0.59135091  0.29779151  0.27689403]
1.9779613018
[ 0.          0.48551953  0.7622897   1.39502275  1.69740105]
[ 0.48551953  0.7622897   1.39502275  1.69740105  1.9779613 ]

ここではPyOpenCLについて解説しましたが、PyCUDAでもgpu_arrayという同等の機能が提供されていて、ほぼ同じように利用することが可能です。

2012-06-04

OpenGL連携を有効にしつつpipからPyCUDA/PyOpenCLをインストールする

Python OpenCL CUDA

PyCUDA / PyOpenCLをソースからインストールする場合には、setup.pyを実行する前にconfigure.pyを実行します。configure.pyが何をしているかというと、以下のようなオプション情報を記述したsiteconf.pyというファイルを生成します。

BOOST_INC_DIR = []
BOOST_LIB_DIR = []
BOOST_COMPILER = 'gcc43'
BOOST_PYTHON_LIBNAME = ['boost_python']
USE_SHIPPED_BOOST = True
CL_TRACE = False
CL_ENABLE_GL = False
CL_ENABLE_DEVICE_FISSION = True
CL_INC_DIR = []
CL_LIB_DIR = []
CL_LIBNAME = []
CXXFLAGS = ['-arch', 'i386', '-arch', 'x86_64', '-isysroot', '/Developer/SDKs/MacOSX10.7.sdk']
LDFLAGS = ['-arch', 'i386', '-arch', 'x86_64', '-isysroot', '/Developer/SDKs/MacOSX10.7.sdk', '-Wl', '-framework', 'OpenCL']

私はよく、PyCUDA / PyOpenCLをインストールする際にはpipを使うことを勧めているのですが、OpenGL連携を有効にしたりCUDA / OpenCLのインストール場所が通常とは異なる場合にはconfigure.pyでそれらを指定する必要があり、pipコマンドを実行するだけではそれらを反映させることが出来ません。今回はその場合の対処方法を紹介します。

ホームディレクトリに.aksetup-defaults.pyというファイルを作成し、上述のsiteconf.pyと同じフォーマットで必要なオプションを記述します。

CL_ENABLE_GL = True

これでpipによるPyCUDA / PyOpenCLのインストールをカスタマイズすることができます。これらはPyCUDA / PyOpenCLのsetup.pyの動作なので、pipだけではなくeasy_installでも同様に動くはずです。

2012-05-30

VimのPyOpenCL用シンタックスハイライト

Python Vim OpenCL

PyOpenCLでプログラムを作るときには、以下のようにPythonプログラムの中にOpenCL C言語のプログラムを埋め込むことが多いですが、Pythonの文字列としか認識されないので埋め込まれたプログラムは見辛くなってしまいます。

prg = cl.Program(ctx, """
__kernel void sum(
    __global const float *a,
    __global const float *b,
    __global float *c
)
{
    int gid = get_global_id(0);
    c[gid] = a[gid] + b[gid];
}
""").build()

PyOpenCLのリポジトリにはvim用のsyntaxファイルが含まれています。これを使うことでPythonの文字列の中に埋め込んだOpenCL C言語を認識してシンタックスハイライトをかけてくれるようになります。実際に使うと以下のようになります。

導入方法はpyopencl.vimをダウンロードして~/.vim/syntax/pyopencl.vimに保存するだけです。これは、Unix系OSの場合ですので、Windowsで使用する場合は適当に読み替えてください。また、OpenCL C言語自体のシンタックスハイライトがopencl.vimにあるので、同様に~/.vim/syntax/opencl.vimに保存しておくといいでしょう。
これで、ファイルタイプがpyopenclのとき、上の画像のように//CL//から始まるdocstringの中がOpenCL C言語として認識されます。ファイルタイプをpyopenclに設定するには、編集時に

:set filetype=pyopencl

とコマンドを打つか、適当な行に

# vim: filetype=pyopencl.python

と書いておくことで自動的に認識させることができます。

2011-10-17

PythonからOCamlにアクセスするための低レベルインターフェースを作ってる

Python OCaml

Python Developers Festaに参加しました。OCamlでPythonの拡張モジュールを作っているという話をしていたら、何人かの人に興味を持ってもらって、色々と意見やアドバイスをいただきました。その時点ではPython/C APIとOCaml/C Interfaceの橋渡しをCで実装していたのですが、ctypesでやってみるのはどうかという話があったので早速やってみました。
OCaml/C Interfaceの全機能の実装はまだ行っていませんが，ソースコードをBitbucketで公開しました。
https://bitbucket.org/likr/otypes

otypes

OCaml/C InterfaceのC関数部分はctypesから呼び出せるので問題はないのですが、他にも数々のマクロが定義されているのでそれらをctypesで地道に実装しています。otypes.pyの一部を抜粋して以下に掲載します。

#!/usr/bin/env python

import ctypes
import sys

class MLDLL(object):
    def __init__(self, name):
        self.lib = ctypes.CDLL(name)
        ArgArray = ctypes.c_char_p * (len(sys.argv) + 1)
        args = ArgArray()
        for i,arg in enumerate(sys.argv):
            args[i] = ctypes.c_char_p(arg)
        args[len(sys.argv)] = None
        self.lib.caml_startup(args)
        self.lib.caml_named_value.restype = ctypes.POINTER(ctypes.c_long)

    def caml_named_value(self, name):
        return self.lib.caml_named_value(ctypes.c_char_p(name)).contents
    
    def val_long(self, x):
        return (x << 1) + 1
    
    def long_val(self, x):
        return x.value >> 1
    
    def val_int(self, x):
        return self.val_long(x)
    
    def int_val(self, x):
        return self.long_val(x)

OCamlプログラムのビルド

以下のようなOCamlプログラムを用意してsample.mlとします。CやPythonなど外部からアクセスしたい値をCallback.registerで設定します。

let f x = x + 3

let sort_list list =
  List.sort (fun a  b -> b - a) list

let sort_array array =
  Array.iter (fun x -> Printf.printf "%d " x) array;print_newline ();
  Array.sort (fun a b -> b - a) array

let _ =
  Callback.register "x" 6;
  Callback.register "y" 7.;
  Callback.register "f" f;
  Callback.register "sort_list" sort_list;
  Callback.register "sort_array" sort_array;
  Callback.register "tuple" ('a','b');
  Callback.register "max3" (fun (a,b,c) -> max a (max b c))

ビルドには-output-objオプションを使ってCからリンク可能なsoファイルを作成します。

ocamlopt -c sample.ml
ocamlopt -output-obj -o sample.so sample.cmx

OCamlMakefileを使う場合は以下のようにMakefileを記述します。

SOURCES = sample.ml
RESULT  = sample.so

OCAMLLDFLAGS = -output-obj

all: native-code

include OCamlMakefile

OCamlプログラムのビルド（Mac）

当方のMac環境では上記の方法ではctypesからのロード時に以下のようなエラーが発生して動かすことができませんでした。

importing mllibrary
Traceback (most recent call last):
  File "./test.py", line 7, in <module>
    lib = MLDLL('sample.so')
  File "/Users/likr/workspace/otypes/otypes.py", line 10, in __init__
    self.lib = ctypes.CDLL(name)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/ctypes/__init__.py", line 353, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: dlopen(sample.so, 6): Symbol not found: _caml_atom_table
  Referenced from: /Users/likr/workspace/otypes/test/sample.so
  Expected in: flat namespace
 in /Users/likr/workspace/otypes/test/sample.so

試行錯誤してみたところ、ccオプションを指定してCコンパイラのオプションを外してやることで正しく動くようになりました。原因がしっかりとわかっていないので要調査です。

ocamlopt -cc="gcc" -output-obj -o sample.so sample.cmx

Makefileの場合はこうなります。

make OCAMLLDFLAGS='-output-obj -cc "gcc"'

サンプル

otypesを使ってPythonからsample.mlの値にアクセスしてみます。

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

import otypes
from otypes import MLDLL

# ライブラリのロード
print 'importing mllibrary'
lib = MLDLL('sample.so')
print

# int値の取得
x = lib.caml_named_value('x')
print lib.int_val(x)
print

# int -> int 関数の呼び出し
f = lib.caml_named_value('f')
res = lib.caml_callback(f, x)
print lib.int_val(res)
print

# float値の取得
y = lib.caml_named_value('y')
print lib.double_val(y)
print

# float値の作成
d = lib.caml_alloc(1, otypes.DOUBLE_TAG)
lib.store_double_val(d, 0.3)
print lib.double_val(d)
print

# int * int の取得
ab = lib.caml_named_value('tuple')
print chr(lib.int_val(lib.field(ab, 0)))
print chr(lib.int_val(lib.field(ab, 1)))
print

# int * int * int -> int 関数の呼び出し
max3 = lib.caml_named_value('max3')
arg = lib.caml_alloc(3, 0)
lib.store_field(arg, 0, lib.val_int(3))
lib.store_field(arg, 1, lib.val_int(8))
lib.store_field(arg, 2, lib.val_int(5))
ret = lib.caml_callback(max3, arg)
print lib.int_val(ret)
print

# Arrayのソート
l = [1,4,2,35,64,3,12,14,46,57]
sort_array = lib.caml_named_value('sort_array')
ml_array = lib.caml_alloc(len(l), 0)
for i,item in enumerate(l):
    lib.store_field(ml_array, i, lib.val_int(item))
lib.caml_callback(sort_array, ml_array)
for i in range(len(l)):
    print lib.int_val(lib.field(ml_array, i))
print

# Listのソート
sort_list = lib.caml_named_value('sort_list')
def cons(hd, tl):
    l = lib.caml_alloc(2, 1)
    lib.store_field(l, 0, hd)
    lib.store_field(l, 1, tl)
    return l
tl = lib.val_emptylist()
tl = cons(lib.val_int(57), tl)
tl = cons(lib.val_int(46), tl)
l = cons(lib.val_int(58), tl)
l = lib.caml_callback(sort_list, l)
print lib.int_val(lib.field(l, 0))
l = lib.field(l, 1)
print lib.int_val(lib.field(l, 0))
l = lib.field(l, 1)
print lib.int_val(lib.field(l, 0))

TODO

OCamlのGC制御
PythonフレンドリーなAPIの提供
OCaml/C Interfaceの未実装部分の開発

2011-09-25

PyOpenCLハンズオン報告

Python OpenCL

kyoto.py PythonハンズオンでPyOpenCLハンズオンのチューターを担当させていただきました。ハンズオン資料を公開しています。

PyOpenCLとはOpenCLという並列計算フレームワークをPythonから扱うためのライブラリで、Andreas Klöckner氏らによって開発されています。OpenCLのメリットや特徴については株式会社フィックスターズさんが公開している資料にわかりやすくまとめられています（OpenCL策定当初なので少々古いですが）。

OpenCLはGPUコンピューティングができるとして注目を集めているのですが、PyOpenCLを使ってPythonからOpenCLを扱えるようになるとまた別のメリットが現れてきます。

タスク並列、データ並列プログラムどちらにも対応するため、従来マルチスレッド、マルチプロセスによる並列化やC言語拡張でPythonプログラムの高速化していた部分の代用になる
PILやPyOpenGLなど、Pythonの豊富な外部ライブラリと連携して比較的簡単にアプリケーションを構築可能
デバイス側で実行されるカーネル関数のコードはOpenCL C言語そのものなので、プロトタイピングに有効

OpenCL、PyOpenCLについてなかなか情報発信ができていない状態なのですが、わからないことなどがあれば気軽に質問でもしてください。

2011-06-12

Pyomoの紹介

Python 最適化勉強会

Python京都勉強会で発表してきましたという報告とその補足です．関西ではPythonの勉強会はまだ多くはないのでいろいろ貴重でマニアックな話を聞くことができました．プログラムはこんな感じでした．
30分のセッション

Python3（http://www.slideshare.net/naoina/python-kyoto-study）
パッケージング
Django
数理最適化とPython（自分：http://www.slideshare.net/likr/python-8277902）
DotCloud

Vim勉強会（http://atnd.org/events/16581）
itertools（自分）
Sphinxとblockdiag（http://www.slideshare.net/lab1092/0611-sphinx-pykyoto）
Vim プラグイン（http://www.slideshare.net/naoina/python-kyoto-study-lt-8283397）

さて，自分の発表したPyomoという数理最適化モデリングのライブラリについてですが，プロジェクトのWebサイトはhttps://projects.coin-or.org/Cooprです．PyomoはCOIN-OR傘下のプロジェクトで，オープンな最適化環境を目指して開発が進められています．最適化問題を解くためには通常ソルバーという種類のソフトウェアを使います．商用ではCPLEXやGurobi，フリーのものではGLPKやSYMPHONYなど様々なソルバーが公開されています．問題を解くときには問題の情報をファイルに記述してソルバーに渡してやったりするのですが，そのようなファイルはソルバーが処理しやすいようなフォーマットをしていて人間が直接読んだり書いたりするには適していません．そのため問題データを出力するときにはモデリングツールという別のソフトウェアを使って問題の情報を記述するのが一般的です．モデリングツールを使うことによってソルバー間のファイル形式の違いを吸収したりもしてくれます．この辺りはプログラム言語の高級言語とマシン語の関係に似ています．

モデリングツールにも様々な種類があり，単独のソフトウェアとして提供されているもの，CやJavaなどの言語のライブラリとして提供されているものなどがあります．Pyomoは後者で，Python用のライブラリになっています．PyomoがPythonを採用した理由の1つとして，Pythonはインタプリタ言語だからコンパイルすることなしにモデルを手探りで修正していけることが挙げられています．Pythonのライブラリとして使用できるため，アプリケーションと統合しやすいこともPyomoの特徴です．

以上，昨日の発表でこのあたりのことを喋りきれていなかったので補足でした．

2011-05-19

SL4A(Python)でQRコード作成

Android Python

SL4Aで簡単なテキストなどをQRコード化する方法を紹介します．こちらのサイト（http://libpanda.s18.xrea.com/qrcode.html）でPythonでQRコードを作成する方法が解説されていますが，Android上ではPILを使ったりswfに出力したりするのは現実的ではありません．
WebViewとCanvasで表示という方法なども考えられるのですが，Android端末にZXing TeamのQRコードスキャナーがインストールされている場合，Intentを送ることでQRコードを表示させることができるのでこれを利用します．（参考：http://awalkingcity.com/blog/2008/10/23/qr-codes-made-even-easier-with-android/）

import android

droid = android.Android()
extras = {}
extras['ENCODE_TYPE'] = 'TEXT_TYPE'
extras['ENCODE_DATA'] = 'test'
intent = droid.makeIntent('com.google.zxing.client.android.ENCODE', None, None, extras).result
droid.startActivityIntent(intent)

これだけでOKです．外部アプリとIntentで連携することで簡単に高機能なスクリプトが書けそうです．