Syoyo Fujita's Blog

raytracing monte carlo

Month: November, 2008

Good introduction and history on Numerical Computation with Guaranteed Accuracy

http://w2.gakkai-web.net/gakkai/ieice/vol6pdf/vol6_09.pdf


This is a must see material on numerical computation.

The document describes a history of Prof. Oishi’s research on numerical computation with guaranteed accuracy.

Prof. Oishi is well known researcher on Numerical Computation with Guaranteed Accuracy(GA is based on interval arithmetic and it includes error free computation, faithful computation, etc)

Theory of Guarantted Accuracy is vital for graphics application for example accurate ray tracing, physics and geometry processing.

[Ja]

lucille のレイトレコアもそろそろしっかりと精度保証して演算誤差をバウンドしたいと思っていて調べていたら見つけました.

大石先生の精度保証に関する研究のイントロダクションと歴史についてかかれたドキュメントです.
短くまとまっていて、また研究がどのような発展を辿ったのかよくわかってとても役に立ちます.

ところで、精度保証の基礎となるアルゴリズムの TwoProduct って FMA(一回だけの丸め) がないと 17 fp 演算も必要なんですよね.
x86 もきちんとした FMA を実装してくれたら(AVX や SSE5 の FMA は二回の丸めなのでダメ), 結構精度保証のテクニックも気軽に使えるものになると思うのです.

Haskell のスコープの仕様でハマる

久しぶりに MUDA を更新していて、
以下のようなコードでハマりました.


test n = muda n

  where

  muda 0 = show 0
  muda n = concat
    [ muda (n-1)
    , ","
    , show bora
    ]

  bora = n

main = putStrLn $ test 3

これを実行すると、 “0, 1, 2, 3” を期待していたのですが、


0,3,3,3

となってしまいます.

C とかと違って、Haskell は関数型だから関数単位でのスコープなので(の解釈で合っているはず)、
bora が参照している n は muda の中の n ではなくて test の n だったのです.

以下のようにして、where でもうひとつ節を作れば期待通りになります.


test n = muda n

  where

  muda 0 = show 0
  muda n = concat
    [ muda (n-1)
    , ","
    , show bora
    ]

    where

      bora = n

main = putStrLn $ test 3

$ runhaskell test.hs
0,1,2,3

じつはこれ、JavaScript を書いているときも同じように変数のスコープでハマりました…
(JavaScript も関数型なスコープ体系なので)

Alchemy too slow…?

LLVM を使って C のプログラムを Flash で動かせるようにした Alchemy がリリースされ、
Alchemy はパフォーマンスが必要なプログラムに向くということで、
ベンチマークを取ってみました.

ベンチマークに使ったのは AmbientOcclusion レンダラの C 版.
(これを今後は AO ベンチマークと呼ぶ事にしよう)

http://lucille.svn.sourceforge.net/viewvc/lucille/angelina/proce55ing/c_reference/


C(gcc4.4 -O3) : 2.6 secs
Alchemy       : 480 secs(8 分)

(参考: 前回のパフォーマンス比較)

う、遅すぎ…
AS3 版 AO の結果よりもさらに 8 倍も遅い…

Alchemy 版は、以下の手順でコンパイルして速度を計りました.


$ alc-on
$ gcc -O3 ao.c
$ time ./a.exe
/Users/syoyo/src/flex_sdk_3.2/bin/adl _sb_19448/app.xml 2> /tmp/adl.trace & echo $!
...

遅さの原因はどこに

Alchemy の gcc を使い、gcc -O2 -c として中間の LLVM ビットコードオブジェクトでコンパイルを止め、
llc -march=c として C に変換して普通の gcc でコンパイル実行したら、
普通のネイティブ(通常の gcc)と同じくらいの早さだったので、
中間言語の時点では大きな速度低下になるような要因はないようです.
(不思議なスタブ関数がいっぱい挿入されるとか)

なので、やはり AVM の性能や、LLVM 中間言語 -> AVM トランスレータあたりに速度低下の原因があると思います.
関数呼び出しのオーバヘッドが大きいとか、JIT が利いていないとか、サンドボクシングでポインタ経由のメモリアクセスにいろいろチェックが入っているとか、なのかもしれません.

Shark で Alchemy で生成したバイナリをプロファイリングしてみましたが、
シンボル情報がないのか JIT 生成された関数なのか、
時間を占める上位関数はアドレス情報しかありませんでしたので、
なにをやっているところが時間がかかっているのかよくわかりません.

ただ、時間を占めている関数のアセンブラを見ると分岐と call ばかりの命令列でしたので、
たぶんやはりスタブ関数呼び出しとかサンドボクシングのチェックとか、インタプリタ実行とか、そこらへんなのかなぁと.

iPhone 3G = 4.8 GFlops?

genki さんと一緒に,
iPhone や iPod touch に載っている ARM プロセッサの VFP(ベクトル浮動小数点ユニット) について調べていました.

とりあえず現時点で理解したこと.

VFP は一度に演算するベクトル長を指定できる.

浮動小数点ステータスレジスタの LEN ビットにベクトル長を書き込むことで、
FP 命令 (fmuls など)はそのベクトル長ぶんの演算を一度に行うという、ちょっと変態なアーキになっている.
通常の FPU としての利用だと、ベクトル長 = 1 ということ.
ベクトル長は単精度だと最大 8 個、倍精度だと 4 個まで設定できる.

ただし、ベクトル長ぶんのレジスタバンク(レジスタペア)を消費する.
単精度だと全部で 32 レジスタ(s0 ~ s31)、倍精度だと 16 レジスタ(d0 ~ d16)あるけど、
たとえばベクトル長 4 なら、各命令のオペランドは [s8, s11] などという風に 4 個レジスタを消費する.
つまり、こんな感じの動作モデルになる.


# ベクトル長  = 1 の場合:
fmuls s8, s9, s10    ; s8 = s9 * s10

# ベクトル長 = 4 の場合:
fmuls s8, s12, s16  ; s8  = s12 * s16
                    ; s9  = s13 * s17
                    ; s10 = s14 * s18
                    ; s11 = s15 * s19

また、レジスタすべてをベクトル演算用に使うことはできなくて、
単精度の場合、s0 ~ s7 のレジスタバンクはスカラ演算のみに使うことができる.
つまり、単精度でベクトル長=4 のときだと、実質的には 6 SIMD レジスタ(4×6=24) しか扱えない.
倍精度に至っては、ベクトル長=4 だと 3 SIMD レジスタしか使えない.
パイプライン実行で命令レイテンシを隠してパフォーマンスを出すことを考えると、
6 レジスタでぎりぎりなんとかストールなしになるかも..、というところか.

VFP で SIMD 演算のプログラミング

通常の浮動小数点演算コードは、iPhone SDK の gcc-arm により、
ベクトル長=1 の条件で浮動小数点命令にマップされるけれども、
ベクトル長 = 4 などの SIMD 演算命令出力はサポートされていない.

llvm の ARM バックエンドも、 型の演算は、単純に 4 回の fp 演算に展開されてしまう.
また gcc で C 言語から VFP 命令を直接扱える intrinsic 関数のインターフェイスもない.
(VFP の後発の NEON だと, intrinsic 関数がある).

なので、現状 VFP で SIMD 演算をしたかったら、アセンブラでプログラミングするほかないようです.
ここはやはりオレコンパイラバックエンドを書ければ理想ですね.

VFP 演算ピークはいくら?

VFP には FMAC(積和命令) があり、レイテンシ/スループットは 1/8 である.
(タイミングは、ベクトル長には関係ないのかな?)

http://infocenter.arm.com/help/index.jsp?topic=/com.arm.doc.ddi0344c/ch16s07s01.html

ベクトル長 = 4, 単精度のときを考える.
iPhone 3G だと 600 MHz くらいのクロックのようなので、 VFP もこの速度で動いているとすると、

0.6(GHz) x 2(FMAC) x 4(SIMD LEN) = 4.8 GFlops

となります.
仮説が正しければ、ポテンシャルとしては、
iPhone 3G は Core2 2GHz の 1/4 くらいの浮動小数点パフォーマンスがあることになります.

けっこう、iPhone でメガデモというのも可能性としてありそうかな?… 🙂

CPU performance of Amazon EC2 instances.

ちょっとそろそろ lucille を使ってクラウドレンダリング(cloud rendering) を
やってみようかと考えているので、まずは Amazon EC2 のインスタンスの CPU パワーはいかほどか、調べてみました.

測定をした CPU インスタンスは, Standard instances から small, それから High-CPU instances から medium と xlarge です.
http://aws.amazon.com/ec2/instance-types/

ベンチマークプログラムは himeno bench,
コンパイラは gcc 4.2.3, コンパイルオプションは
-mfpmath=sse -msse2 -mtune=native -O3 -DSMALL
として、SSE 命令を明示的に使うようにしています(No more x87!).

結果.


Instance          MFLOPS(larger is faster)
========          ======
small              454.6
c1.medium         1016.9
c1.xlarge(64bit)  1416.6

クライアントマシンとの比較として、
Intel Mac Core2 2.16 GHz, gcc 4.2.1 での結果をとってみました.
こちらは、 908.9(32bit), 1056(64bit) でした.

High-CPU medium instance がだいたい Core2 2GHz 相当(1 core あたり)と言えそうです.
small はちょっと非力ですね.

medium は仮想 2 コアサポートなので、うまくマルチスレッドプログラムを作れば small に比べて、
(カタログスペック通りに) 4 倍のパフォーマンスを得られる可能性があります.

なので、数値計算が主な利用用途ですと、
medium は $0.2/hour, small は $0.1/hour なので、
medium を CPU インスタンスとして使うのがお得と言えそうです.

もちろん、計算をガンガンしたい場合は xlarge という手もありますが、 $0.8/hour でちょっと高いのが難点ですね.

RenderMan BOF in Japan.

http://info.rendersan.org/2008/11/rendersan3.html

RednderMan BOF in Japan is scheduled on Nov 22.
I’ll do a talk on lucille in this BOF:

TITLE:
luiclle global illumination renderer
RenderMan capable renderer and Beyond.

[Ja]

11/22 に、RenderSan という RenderMan 関係野郎の集まりにて、
lucille について話します.

RenderMan サポートしている lucille レンダラの紹介と、
RenderMan をサポートするのを越えて, これからの lucille の野望などについて話ます.

集まりの詳細はリンク先をご参照ください.

iRender : Global Illumination with iPhone.

irender.png

I’ve ported ambient occlusion renderer to iPhone as a OpenGL ES application.

http://lucille.svn.sourceforge.net/svnroot/lucille/angelina/irender

irender.tar.gz (project files)

The renderer core itself is simple and I already ported it to many platforms.

Proce55ing
Flash10
JavaScript
Performance comparison

Unfortunately I haven’t real iPhone, so I just tested that the renderer runs correctly in the simulator.
I don’t know how much the rendering time is in real iPhone device.

[Ja]

iPhone にアンビエントオクルージョンレンダラをポーティングしてみた.
とりあえず OpenGL ES が使えるということで、iPhone の開発環境(Cocoa?)も ObjC もまったくわからない状態で,
レンダラ部分を無理矢理 C + GL だけ使って組み込むことで実現してます.
2 時間くらいで作ったのでやっつけです. ソースもいいかげん.

残念ながら実機がないので速度がどれくらいかは分かりません.
シミュレータだと結構早いです. Core2 2.15GHz で、ネイティブ C が 2.5 秒なのに, iPhone simulator だと 5 秒ほど.

あ、ちなみにレンダリングは一回やるだけで、あとはテクスチャにして回転させているだけです.
毎回回転するたびにレンダリングを行っているわけではありません.

iPhone は 1,000 万台も普及しているとのことなので、
分散レンダリングを iPhone でやるのも面白いんじゃないかなぁと思っています.

iPhone で画像テクスチャを貼付ける

GLES では glDrawPixels が無いということが分かったので、
今回レンダリング結果の画像表示はテクスチャ経由で行っています.

参考までに、 iPhone で OpenGL ES を使ってテクスチャを貼付けて表示する部分はこうなってます.
(タブの設定のせいでインデントがずれているところがあります)

テクスチャを使う手順はこんな感じ.

– テクスチャになるデータを配列で用意
– (最初だけ)テクスチャオブジェクトを作成 (glGenTextures)
– テクスチャオブジェクトをバインド (glBindTexture)
– (最初だけ)テクスチャ拡縮フィルタを設定 (glTexParameterf)
– (更新するときだけ)テクスチャになるデータをテクスチャメモリに転送 (glTexImage2D)
– テクスチャマッピングを有効にする (glEnable(GL_TEXTURE_2D))
– テクスチャ座標アレイを有効にする (glEnableClientState(GL_TEXTURE_COORD_ARRAY))
– テクスチャ座標ポインタを設定 (glTexCoordPointer)
– ポリゴン描画

矩形ポリゴンに矩形テクスチャを貼付けています. テクスチャ座標は texCoords で与えています.
テクスチャデータは width*height*4 bytes の GLubyte 型の配列で用意すれば OK.
よって配列をあとは自由にいじればピクセル処理を行い、それを GLES を通して表示することができます.

GLES だとテクスチャデータの更新は glTexImage2D でやるしかないのかな?
それだと更新コストが高いと思うので、あまり毎フレーム画像をアニメーションさせるというのには向かないかもしれません.

それにしても、しばらくぶりに OpenGL を使いました(SW raytracing に移行しつつあるので OpenGL を使う理由がなくなってきている).
以外と API やら中身の動作様式を覚えているものです.

[追記]

早速実機を持っているひとに試してもらいました.

iPhone3G -> 130 秒 

う、遅すぎ…

でも実は、こんなのが原因のようです.

浮動小数点演算の高速化
http://d.hatena.ne.jp/kstn/20081031/1225458244

というわけで、-mthumb オプションを無効にして再度計測してもらいました.

130 秒(-mthumb あり. デフォルト) -> 38 秒(-mthumb オフ)

結構早くなりました.

ただ、-mthumb を付けると(デフォルトの動作)、コードサイズは小さくなるという効果があるとのこと.
多くのアプリは浮動小数点演算は多用しないので、このオプションがデフォルトで有効になっているようです.

やはりパフォーマンスについて調べるときにはアセンブラレベルでちゃんと確認したいですね.
残念ながら実機がないと、実機向けコンパイル(armv6 ターゲット)の段階でライセンスが何たらと出てコンパイルがまず走らないので、アセンブラ出力が確認できません… 実機を手に入れるしかないのか…

ちなみに CPU(armv6) 自体には VFP と呼ばれる SIMD 命令があるので、これを使えるようになるともっと早くなるかも.
でも iPhone SDK 経由で使えるのかな?(たとえば、インラインアセンブラを使ってとか).

VFP 向けの SIMD 数学ライブラリはありました.
http://code.google.com/p/vfpmathlibrary

インラインアセンブラを使っています.なのでたぶんユーザ側でもきっと VFP を使えるはず.

MUDA の VFP バックエンドをちょっと書きたくなってきた.
でも今は MUDA は C コードを吐くので、 gcc 側でイントリンジックスなりで VFP に対応してくれていないとダメなんですよね.
AVX の時も思ったけど、やはり gcc とかのコンパイラの対応度に依存しないように、
いずれは直接最適化されたアセンブラを吐く、自作の SIMD 命令むけコンパイラを書いたほうがいいんだろうな.

[Problem] Vector swizzle in AVX.

Swizzling vector element(e.g. v.xwvz) in AVX requires a bit sophisticated solution because there’s no single instruction to do it.

We have to swizzle vector element with combination of shuffle, permute, blendv instructions.

It is like a puzzle solving, thus I need time to solve it optimally.

To derive my own solution, I’ve searched and collected some existing works.

Franz Franchetti and Markus Püschel
Generating SIMD Vectorized Permutations
Proc. International Conference on Compiler Construction (CC), Lecture Notes in Computer Science, Springer, Vol. 4959, pp. 116-131, 2008
http://spiral.ece.cmu.edu:8080/pub-spiral/abstract.jsp?id=95

PerfectShuffle from LLVM
http://llvm.org/svn/llvm-project/llvm/trunk/utils/PerfectShuffle/

[Ja]

v.xxxx を AVX でやる方法は作ったけど、もうちょっと一般化して、
AVX(double x 4)でベクトルのスウィズルをやるのは実はちょっと頭がいる.

スウィズルとはこんな感じのベクトル要素入れ替えのことです.


vec v = (1, 2, 3, 4);
vec w = v.xzwy;              // w = (1, 3, 4, 2)

AVX の場合、128bit の制約があるのでデータを自由に動かせないなどいくつかの制約があるので、
その制約の中で命令の組み合わせで最適なものを探しだす必要がある.

レジスタの内容をメモリにいったん書き戻して要素を並び替えて読み戻しとかいうのはここでは無しね.
それにこの方法だとパフォーマンスもよくないだろうし.

生真面目に手で求めるとすると, 4^4 = 256 通りあるわけで、 4 要素でなくて 1, 2, 3 要素のスウィズル(v.x, v.yz, v.xzw など)まで考えると 340 通りもある. さらに float x 8 や、将来の double x 8 などを考えると組み合わせパターンは爆発して手に負えなくなる(か、サルになってコードをずっと書き下すかだ).

そんなわけで手持ちの命令セット(shuffle, perm, blendv あたり)からうまく最適な組み合わせでスウィズルを行うコードを作るのはちょっとした最適化問題やパズル解きのようなことをしないといけない.

自分でその解法を一から考えるのも時間の無駄だし(考えるの止めたわけではありません)、
ちょっと解法のアイデアはあるけど最適な気がしないので、
より最適な解法はないかちょっくら web をあさってみたら、やはり当たり前のように既存研究がある.

しかも spiral プロジェクトと LLVM から!
解き方は、どちらも動的計画法を使うようですね.

コンパイラバックエンドでレジスタ割り当てや命令スケジューリングの一手法である
ILP(Iterated Linear Programming)とちょっと似ているかな.

Eagle eye is rendererd with Arnold

映画イーグルアイで GI レンダラ Arnold が使われたということで、観てきました.

http://www.eagleeyemovie.com/intl/jp/


ほーほー、なるほどこいつに使われているわけね、
これだとたしかに prman じゃやりにくそうだから Arnold にしたのかな?
なるほどなるほど…

より詳しいメイキングは、ここらへんで手に入ります.

http://vfxworld.com/?atype=articles&id=3779

http://www.fxguide.com/article501.html

http://digitalcontentproducer.com/mil/features/step_eagleeye/

メイキングを見ると、あのシーンでも使われていたようです.

ちなみに、とあるレンダラ野郎に言わせれば、Arnold は漢なレンダラとのこと.
まあたしかにそういわれればそんな気もします.
下手な小細工しないで、一直線なレンダラが長期的にはいいのかもしれません.

Introduction to type inference

JS 実装をしようかな、と思い立ち、
まずは型推論の知識をしっかり取得するところから初めようと考えています.

そもそも、型推論とはなんぞや、というところからおさらいしてみます.

私が型推論について知りたいことの一番の目的は、
「型推論を行うことで動的言語の事前コンパイルを可能にしプログラムの実行を早くすることができるのではないか」,
ということになりますので、それに特化した内容になっています.
もちろん型推論には、プログラムを早くするという以外の目的もありますが、
ここではそれらは取り上げないことにします.

型推論(type inference)とは、簡単に言うと、


var muda = 7

というコードがあったら、

「変数 muda って Int 型(整数型)じゃね?」

とコンパイラがよろしく変数の型を決定してくれる(推測してくれる)機能になります.

静的な関数型言語(Ocaml, Haskell)などでは、変数(関数)には型を通常プログラマが指定しませんが、コンパイラが型推論を用いて型を割り付け、文法にエラーが無いかの判断などに使ったり、最適なアセンブラコードの選択に使ったりします.

動的言語でも、この型推論を行うことで、実行時にエラーが起きるようなスクリプトを、あらかじめ実行前にエラー判定したり、実行前に型が不変であることが分かれば事前にコードパスを特定の型に特化してコンパイルすることで、実行時の処理速度を上げることができます.

たとえば,


def muda(a):
  print a

muda(10)

というようなスクリプトがあった場合、muda() 関数の引数 a の型は int しかとらないことが分かるので、int に特化したネイティブコードをあらかじめスクリプト実行前(JIT コンパイルでもよし)に生成することができます.

型推論のアルゴリズム

では実際に、型推論ってどうやればいいのでしょうか?
推論っていうくらいだから、なにか AI みたいなことをするのでしょうか?
(実はしません)

ひとまず型推論のアルゴリズムについて調べてみました.

Hindley-Milner

関数型言語向けでは、Hindley-Milner のアルゴリズムが基本になるそうです.

Damas, Luis & Milner, Robin (1982), “Principal type-schemes for functional programs”, POPL ’82: Proceedings of the 9th ACM SIGPLAN-SIGACT symposium on Principles of programming languages, ACM, pp. 207–212.
http://groups.csail.mit.edu/pag/6.883/readings/p207-damas.pdf

イメージとしては再帰的に式をたどって型を割り当てていく感じ?
py2llvm で実装した型推論の方法と似ているかな?
(こういう車輪の再発明なことをして時間の無駄をしたくないから、既存研究を調べるわけです.
先人たちのほうがきちんと sophisticate された体系としてまとめてくれているわけですから)

CPA

オブジェクト指向言語向けには、CPA と呼ばれる Cartesian Product Algorithm から入っていくのがよさそうです.

http://research.sun.com/self/papers/cpa.html

イメージとしては、取りうる型の組み合わせをどんどん狭めていってひとつの型を導きだすという感じでしょうか.
たとえば最初は変数 muda は [String, Number, Boolean] のいずれか、とセットアップしておき、
いろいろ周りの状況をみて候補を絞っていって、 muda : Number とするような感じ?

Java や JavaScript(?) の基礎になった Self 言語の型推論エンジンとして使われたという実績があります.

また、Starkiller と呼ばれる Python 向けの型推論にも使われているアルゴリズムです.

http://salib.com/writings/thesis/thesis.pdf

DDP

オブジェクト指向言語向けにもう少し発展した型推論として、DPP というアルゴリズムがあります.

DDP: Demand-Driven Analysis with Goal Pruning
http://www.lexspoon.org/ti/

SmallTalk の型推論向けに開発されました.
DDP は CPA よりも、大規模なプログラムに対してスケールしやすいことを挙げています.

DDP ベースの型推論を Ruby に組み込もうとした SoC があり、
http://soc.jayunit.net/
成果は Ruby 向けコード補完として Aptana に取り込まれたそうです.
(ところで、この SoC プロジェクトのメンターだったのは Chris Williams なのですね.
JDK6 の JIT を書いたひと. いろんなところでつながっているんだなぁ.)

まとめ

JS 実装という観点からは、まずは CPA からいろいろ調べていくとよさそう.