GPUは楽し！
トゥットゥルー♪ まゆしぃ☆だいすきなかえるのクーです。
GPUを買ったので遊びます。
今日はchainer もしくはCuPyの「elementwise」に入門します。
PythonからCUDAが簡単に使えるとのことで、とてもうれしいです。
その前に、CUDAとは何なのか少しだけ勉強しました。
CUDAプログラミングの基礎 パート１(PDF)スタックとメモリ
グリッド、ブロック、次元、インデックスなどが重要なことがわかりました。
諸先輩方に教えてもらいます。感謝。
elementwiseを使ってみます
デバイス(GPU)側で並列処理を行う部分の 計算/ロジックを記載した関数のことを「カーネル関数」というそうです。このカーネル関数の部分を記載するのが、elementwiseです。
chainer.cudaのelementwiseはこちら
chainerの中ではCuPyではなく chainer.cudaのelementwiseを使用しています。
CuPyのElementwizeKernelはこちら
カーネルの中で使える関数は本家CUDAを参照します。
まずはMath APIを使うと思うので貼っておきます。
CUDA Math API :: CUDA Toolkit Documentation
メモリ割り当てとか、グリッド、ブロックとかを意識せずに、ほぼ計算だけを書けばよいのは本当にすごいです！！！
スライドも読みました
本家chainer(5.0.0)でelementwiseを使用している箇所は(chainerのソースの中で)ざっとの検索で128か所ありました。ifやforなどの分岐のないものがほとんどのようですが、中には結構複雑なものもあり、勉強中です。
とりあえず、CuPy本家のチュートリアルをやってみました
User-Defined Kernels — CuPy 5.0.0 documentation
import cupy as cp
import numpy as np
if __name__ == '__main__':
# カーネルオブジェクトを変数に入れます。機能の名前でよいと思います。
squared_diff = cp.ElementwiseKernel(
'float32 x, float32 y', #１番目の引数:in_params:入力を定義します。(複数あればカンマ区切りのリストで)
'float32 z', #２番目の引数:out_params:出力を定義します。
'z = (x - y) * (x - y);', #３番目の引数:operation: CUDA-C/C++の表記で計算内容を記載します。
'squared_diff') #４番目の引数:name:内部で使用する関数の識別です。
x = cp.arange(10, dtype=np.float32).reshape(2, 5) # x:{ndarray}[[0,1,2,3,4]\n[5,6,7,8,9]]
y = cp.arange(5, dtype=np.float32) # y::{ndarray}[0,1,2,3,4]
z = squared_diff(x, y)
print(z)
[[ 0. 0. 0. 0. 0.]
[25. 25. 25. 25. 25.]]
C言語/CUDAのコーディングと比較しても夢のようなシンプルさです。素晴らしい！！！
float32の部分をTにするとジェネリック(なんでもあり？）になります。
カーネル関数の中では 3項演算、if や for や while などの制御や、CUDAで準備されている各種の数値演算関数が使えます。
if を使ってみました
本家チュートリアルはあまりバリエーションがないので、諸先輩の説明を参考にしてやってみました。
import cupy as cp
import numpy as np
if __name__ == '__main__':
add_evenzero = cp.ElementwiseKernel(
'T x, T y',
'T z',
'''
if (int(x + y) % 2 == 0) {
z = 0;
} else {
z = x + y;
}
''',
'add_evenzero')
x = cp.arange(10, dtype=np.float32).reshape(2, 5)
y = cp.arange(5, dtype=np.float32) * 2
z = add_evenzero(x, y)
print(z)
[[ 0. 3. 0. 9. 0.]
[ 5. 0. 11. 0. 17.]]
ブロードキャスト((配列に足りない何かを補って完全なる（求められる）配列にしてから計算すること。人類補完計画の配列版みたいなもの(ウソ)))も普通に使えています。
複数行の文字列はシングルクォート３つ(''')で囲むのがPythonのお約束のようです。
3項演算を使ってみました
すいません。これ以降、こんなのが続きます。まぁ、入門なので。。。
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.arange(9, dtype=np.float32).reshape(3, 3)
# ElementwiseKernel のオブジェクトを変数にいれず、そのまま使用することもできます。
y = cp.ElementwiseKernel(
in_params='T x', # パラメータ名を指定しても書けます。
out_params='T y',
operation='y = int(x) % 2 ? x : 0;',
name='add_conditional')(x)
print(y)
[[0. 1. 0.]
[3. 0. 5.]
[0. 7. 0.]]
パラメータ名を指定して書けますが、 in_params, out_params, operation, name は(たぶん)必須だと思うので、この４つはパラメータ名を書かないほうがすっきりすると思います。
カーネルオブジェクトは単発使用の場合は変数に入れないほうがすっきりするかもしれません。（このあたりは好みがでるかもです）
ちなみにchainer本家はカーネルオブジェクトを変数に入れずに使っているケースがほとんどでした。
raw指定をするとインデックスで他の要素を参照できます
ここまでは、CuPyが要素固有の番号などを隠して意識せずに記述できるようにしてくれていましたが、配列の１要素が「他の要素」を参照したい場合、はインデックスを使用します。raw 指定をすることにより、それが可能になります。
自分がカーネルを実行する１つのCUDAコアになったつもりで考えればイメージがわきやすいかもしれません。
raw指定をすると指定した引数に対して、iという特殊変数が使えるようになります。
また、_ind.size()
がインデックスの総数を表します。該当する配列のインデックス指定内で使用するところがポイントだと思います。
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.arange(9, dtype=np.float32)
y = cp.ElementwiseKernel(
'raw T x, int32 shift', # x に raw を指定
'T y',
'y = x[(i + shift) % _ind.size()];', #rawを指定したパラメータで _ind.size() と i が使える
name='array_shift')(x, 2)
print(y)
実行するとエラーがでました！
ValueError: Loop size is Undecided
とのことです。先輩方の説明ではインデックス処理時には、出力する配列の大きさはあらかじめ確定いる必要があり、それを伝えるために出力オブジェクトを引数としてあたえてやる必要があるそうです。
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.arange(9, dtype=np.float32)
y = cp.zeros_like(x, dtype=np.float32) # yにあらかじめオブジェクト(インスタンス)を設定
y = cp.ElementwiseKernel(
'raw T x, int32 shift', # x に raw を指定
'T y', # 入力にraw指定があると出力も引数指定が求められます。
'y = x[(i + shift) % _ind.size()];', #rawを指定したパラメータで _ind.size() と i が使える
name='array_shift')(x, 2, y) # 引数に y を追加
print(y)
[2. 3. 4. 5. 6. 7. 8. 0. 1.]
うまくいきました。
x = cp.arange(9, dtype=np.float32).reshape(3,3)
にすると
[[2. 3. 4.]
[5. 6. 7.]
[8. 0. 1.]]
が出力になります。
i や _ind.size()は、配列の全体数を示していることがわかります。
あれ？yはインスタンスで渡してGPU計算側のデータを詰めてもらっているので、yに改めて代入する必要はないのでは？その通りで、
y = cp.ElementwiseKernel(
は
cp.ElementwiseKernel(
とだけ書けばOKでした。
じゃぁ出力をrawにすると？？？
出力をraw指定してみました。
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.arange(9, dtype=np.float32).reshape(3,3) + 2 # x:{ndarray}[[2. 3. 4. ]\n[5. 6. 7.]\n[8. 9. 10.]]
y = cp.zeros_like(x, dtype=np.float32)
cp.ElementwiseKernel(
'S ind', 'raw T r', # Sは要素をインデックスとして使用する場合のお約束？
'r[ind] = i;',
'shift_indices')(x, y)
print(y)/pre>
[[0. 0. 0.]
[1. 2. 3.]
[4. 5. 6.]]
配列のインデックス範囲を外れてもエラーにならない！！！（マイナスにしても同様にエラーにならず）これはありがたいです。
同時に同じ出力要素に対して計算しようとすると
各CUDAコアは並列に計算しているので、以下のようなコードは計算結果が正しくありません。
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.arange(9, dtype=np.float32).reshape(3,3) + 2 # x:{ndarray}[[2. 3. 4. ]\n[5. 6. 7.]\n[8. 9. 10.]]
y = cp.zeros_like(x, dtype=np.float32)
cp.ElementwiseKernel(
'S ind', 'raw T r', #Sは要素をインデックスとして使用する場合のお約束？
'r[0] += ind ;',
'shift_indices')(x, y)
print(y)
[[9. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
明らかに不正解です。各コアがよってたかってr[0]を更新しようとして、競合が発生していると思うのですが、エラーにならないんですね。バグがわからないので、これは困ります。注意が必要です。
そんなときのためのatomic
アトミック侍。。。じゃなかった、atomicAddとかatomicXorとかのatomicシリーズを使えば、競合を防いで処理してくれるそうです。競合した場合スレッドがどちらか待たされることと、競合のチェックをするので、実行コストはそれなりに高くなるはずです。可能ならば使わないでおきたい機能です。
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.arange(9, dtype=np.float32) + 2 # x:{ndarray}[[2. 3. 4. ]\n[5. 6. 7.]\n[8. 9. 10.]]
y = cp.zeros_like(x, dtype=np.float32)
y = cp.ElementwiseKernel(
'T x', 'raw T r',
'atomicAdd(&r[0], x) ;',
'atomic_add0')(x, y)
print(y)
[54. 0. 0. 0. 0. 0. 0. 0. 0.]
加算したい配列のアドレス(&で指定)を渡すのがポイントのようです。
atomicAddで値の数を数えてみます。
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.array([1,2,3,2,3,4,2,3,5,3,3,7,5])
y = cp.zeros(int(x.max()) + 1, dtype=np.int32)
# 以下ElementwiseKernelはCuPyの説明スライドから
cp.ElementwiseKernel(
'S x',
'raw U bin',
'atomicAdd(&bin[x], 1) ;',
'bincount_kernel')(x, y)
print(y)
[0 1 3 5 1 2 0 1]
ちゃんと数えてくれています。
preambleを使ってみます
共通で使うようなCUDA-C/C++ の定義を別に記述しておけます。
import cupy as cp
import numpy as np
if __name__ == '__main__':
_preamble = '''
template <typename T> __device__ T sigmoid(T x) {
const T half = 0.5;
return tanh(x * half) * half + half;
}
'''
x = cp.arange(9, dtype=np.float32).reshape(3,3)
beta = cp.ones_like(x, dtype=np.float32) * 0.1
y = cp.ElementwiseKernel(
'T x, T beta', 'T y',
'''
T bx = beta * x;
y = x * sigmoid(bx);
''',
'swish_fwd', preamble=_preamble #必ず使うパラメータ以外はパラメータ名指定をするがよいと思います。
)(x, beta)
print(y)
[[0. 0.5249792 1.099668 ]
[1.7233275 2.3947506 3.1122968]
[3.8739376 4.6773143 5.519796 ]]
forを使ってみます
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.arange(9, dtype=np.float32)
y = cp.zeros_like(x, dtype=np.float32)
cp.ElementwiseKernel(
'raw T x', 'T y',
'''
for(int j = 0; j <= i; j++) {
y += x[j];
}
''',
'cumulative_sum')(x, y)
print(y)
[ 0. 1. 3. 6. 10. 15. 21. 28. 36.]
スレッド毎にまわる回数が異なる、あえてひどい例で書きました。
forは どうしても使わなければならない時だけ使うものだと思います。
CUDAで準備されている関数を使ってみます
minとmax を使った例です。いろいろそろっています。上にも張りましたが、あらためてCUDAで使える算術APIをはっておきます。
CUDA Math API :: CUDA Toolkit Documentation
import cupy as cp
import numpy as np
if __name__ == '__main__':
x = cp.arange(9, dtype=np.float32) * 0.4 - 0.5
cap = cp.ones_like(x, dtype=np.float32) * 2.0
y = cp.ElementwiseKernel(
'T x, T cap', 'T y',
'y = min(max(x, (T)0), cap);',
'clipped_relu_fwd')(x, cap)
print(y)
[0. 0. 0.3 0.70000005 1.1 1.5 1.9000001 2. 2. ]
（これぐらいかなぁ。。。ほぼ、諸先輩のトレースなのですが、自由研究は自分でやってみることに意義があります！）
reductionは現在全くわからず！
ReductionKernel についてはまだよくわかっていません。並列計算したものをまとめあげる処理のイメージはあるのですが。。。資料を読んで学習します。
(追記)資料を何度もみているとだんだんイメージがわいてきました。時間ができたら入門してみます。
cudaについてしっかりと説明してくれていそうな資料がありました。時間がかかりそうですが、読んでみます。
*1:平気に貼っていますが、英語を読めるわけではありません。 図や表やわかる単語をみていれば、なんとなく想像できるレベルです。
