![[image]](http://practical-scheme.net/images/gauche-ao.png)
![[image]](http://practical-scheme.net/images/gauche-pt.png)
Shiro(2009/02/23 01:52:32 PST):
ambient occlusionとかpath tracerをいろいろな言語でやってみるというのが はやっていたらしい： AO bench
少し出遅れた感もあるが、以前から延ばし延ばしにしてきた最適化を仕込むちょうどいい 機会なのでGaucheにポートしてみた。
なお既に(R6RS) Schemeへは .mjt さんの手によってポートされている。
yharaさんによるGaucheへのポートもあるが、今回は.mjtさんのコードを出発点に、Gauche特有のオプティマイズをかけてゆくことにした。 最終的なコードはこちら→aobench.scm。 下にも掲載する。
ambient occlusionとpath tracerとで共通点が多いのでひとつにまとめた。 コマンドラインオプションで切り替えられるようになっている (-tオプション)。 また、レイトレースは並列化がしやすいので、複数スレッドで並列に走らせることもできるようにした (-pオプション)。
今回の目的は計算中心のアプリに向けてのGaucheのランタイムのチューンなので、 かなりGauche特有の最適化がコードにも入っている。一応アルゴリズムは変えないように 気をつけた。
0.8.14 rev6565 speedup
-------------+---------------------+---------------------+--------
AO single real 1m59s (119s) real 1m25s (85s) x1.40
user 1m59s (119s) user 1m25s (85s)
AO parallel real 45s (45s) real 26s (26s) x1.73
(4 threads) user 1m54s (114s) user 1m27s (87s)
PT single real 10m51s (651s) real 8m36s (516s) x1.26
user 10m51s (651s) user 8m35s (515s)
PT parallel real 3m56s (236s) real 2m31s (151s) x1.56
(4 threads) user 10m46s (646s) user 8m42s (522s)
-------------+---------------------+---------------------+--------
parallelでの伸びが大きいのは、trunkではflonum allocation最適化のおかげで ごみの量が大幅に減っているから だと思う。GaucheのGCはまだparallel markをイネーブルしていないので、 GCが走っている間は他のスレッドが止まってしまう。 user/realの比が並列度の目安になると思うが、0.8.14ではまだ4コアを 充分に使いきれない。
0.8.14 rev6565 ------------------------+------------+------------+ AO parallel user/real 2.53 3.36 PT parallel user/real 2.74 3.46 ------------------------+------------+-----------+-
trunkでも並列度に不満はあるが、今後parallel markも試してみよう。
0.8.14以降に入ったGaucheコアの最適化にはこんなものがある。
最後のやつは、f64vector-refとかがVM instructionでサポートされるようになったということ。
0.8.14まで、uniform vectorは完全に拡張モジュールになっていたのだけれど、 他の拡張モジュールでバイナリデータを扱いたい時に必要になることが多く、 拡張モジュールのコンパイル時の依存関係がややこしくなっていた。 なのでかねてから、uniform vectorのCレベルでの構造体定義やコンストラクタを libgauche.so側に取り込んで、gauche.uvector拡張モジュールは拡張Scheme APIを 提供することにしようと思っていたのだった。
今回CレベルでのScmUVectorの定義がコアに移ったことで、uniform vectorの 参照もvector-refと同じようにVMインストラクションに展開できるようになった。 以前は数値計算にuniform vectorを使うとrefの性能が足を引っ張って、 むしろ普通のvectorの方が良いくらい、という問題があったのだけれど、これで解消。 しかも普通のvectorに浮動小数点数を入れちゃうとflonum allocationの最適化の 恩恵を受けられないのだけれど、今回の変更でf64vector-refとかが気軽に 使えるようになるとflonumの出し入れにヒープアロケーションが不要になるので 大変都合が良い。実は上の3項目のうち後ろ2つは、最初のflonum最適化の メリットをより大きく引き出すのにも役に立っているのだ。
aobench.scmにかけたGauche特有のチューニングについてざっとメモしておく。 仕事でも、性能が必要なクリティカルパスではだいたい同じような感じのコードになるんで、 Gauche使っててもう少し性能を出したいという人には役にたつかも。
ただし、こういうチューニングをすると見づらいし保守しにくいコードになるので、 あくまで最後の性能追い込みで、ホットスポットだけに適用するという原則を 忘れずに。これらのチューニングをかけても性能はせいぜい倍になる程度。 もともと実行時間の1%しか占めないようなルーチンにこの手の最適化を 施しても全体で0.5%速くなるだけだ。
基本的な戦略は、
以降、コードに説明をはさむ形で。
;; based on http://lucille.atso-net.jp/aobench/
;; and http://d.hatena.ne.jp/mjt/20090206/p1
;; http://d.hatena.ne.jp/mjt/20090209/p1
(use gauche.uvector)
(use gauche.threads)
(use gauche.parseopt)
(use srfi-1)
(use srfi-27)
(use srfi-42)
(use math.const)
(use util.match)
;;-------------------------------------------------
;; Common utilities
;;
;; vector stuff
(define vec f64vector)
(define-inline (vx v) (f64vector-ref v 0))
(define-inline (vy v) (f64vector-ref v 1))
(define-inline (vz v) (f64vector-ref v 2))
(define-macro (vset! v . args)
(match args
[(x y z) `(begin (f64vector-set! ,v 0 ,x)
(f64vector-set! ,v 1 ,y)
(f64vector-set! ,v 2 ,z))]
[(vv) `(f64vector-copy! ,v ,vv)]))
(define vdot f64vector-dot)
(define vscale f64vector-mul)
(define (vcross a b)
(let ([ax (vx a)] [ay (vy a)] [az (vz a)]
[bx (vx b)] [by (vy b)] [bz (vz b)])
(vec (- (* ay bz) (* by az))
(- (* az bx) (* bz ax))
(- (* ax by) (* bx ay)))))
(define vadd! f64vector-add!)
(define vsub! f64vector-sub!)
(define vmul! f64vector-mul!)
(define vdiv! f64vector-div!)
(define (vnormalize! a)
(let1 len (sqrt (vdot a a))
(when (> len 1.0e-6) (vdiv! a len)))
a)
まずはベクタ関係のユーティリティ。define-inlineは匙加減が難しいのだけれど、 bodyが極めて小さな式の時は効果があることが多い。 body内でletで環境を作ったりしているとむしろインライン展開が逆効果になる。 (VMレベルでは、(固定長引数の)関数呼び出しは極めて軽い。)
vset!をマクロにしているのはどっちかというと見やすさのため。要素ごとに 個別指定する場合と、ベクタ全体を代入する場合とでいちいち関数名を変えないで 済むようにと。結局ほとんどベクタ全体の代入しか使わなくなっちゃったけど。 なおこのマクロは、(vset! v x y z)と使った場合にvが複数回出力に現れるので、 あんまり良くない書き方である。みんなで使いまわすようなマクロならvを一時変数で 受けるように書くのが定石。ただし今回のような、ここでしか使わないことが わかっているような場合はまあokでしょう。letかますと少しオーバヘッドがあるし。
(let1 tmp (gensym)
`(let1 ,tmp ,v
(f64vector-set! ,tmp 0 ,x)
(f64vector-set! ,tmp 1 ,y)
(f64vector-set! ,tmp 2 ,z)))
vcrossについては、先に要素を全部出しちゃうのと、(- (* (vy a) (vz b)) (* (vy b) (vz a))) みたいに計算するのと両方試したんだけどこっちの方が速かった。 状況によって変わると思うのでその都度プロファイルすべし。 今回、vcrossが結構大きな時間を占めていたのでいじったけれど、そうでなければ そんなにこだわる必要はないだろう。
びっくりマークがたくさんついていることからわかるように、 本体は副作用ばりばりである。覚悟されたし。
;; ray
(define org car)
(define dir cdr)
(define new-ray cons)
(define (copy-ray ray)
(new-ray (f64vector-copy (org ray)) (f64vector-copy (dir ray))))
;; clamp for final pixel value
(define (clampu8 f)
(clamp (floor->exact (* f 255.5)) 0 255))
;; calculates basis vectors aligned to N
(define (orthoBasis n)
(let* ([v (cond [(< -0.6 (vx n) 0.6) '#f64(1.0 0.0 0.0)]
[(< -0.6 (vy n) 0.6) '#f64(0.0 1.0 0.0)]
[(< -0.6 (vz n) 0.6) '#f64(0.0 0.0 1.0)]
[else '#f64(vec 1.0 0.0 0.0)])]
[s (vnormalize! (vcross v n))])
(values s (vnormalize! (vcross n s)) n)))
;; returns a unit vector points to a random direction
(define (random-direction)
(let ([r (random-real)] [phi (* 2.0 pi (random-real))])
(let1 rho (sqrt (- 1.0 r))
(values (* (cos phi) rho) (* (sin phi) rho) (sqrt r)))))
レイは起点と方向をペアにするだけ。.mjtさんのコードそのまま。 ただし、副作用でレイを変更することがあるので、copy-rayが必要になった。 あとはちょっとしたユーティリティ。
3項以上の比較演算子は0.8.14では2項にくらべてがくんと速度が落ちてたんだけど trunkではそうでもない。 なお、orthoBasisが返すベクタのうち、最初の2つはfreshにアロケートされたもの、 最後のやつは引数そのまま、というのは後で重要になる。
random-directionは共通部分をくくり出しただけ。(sqrt (- 1.0 r)) を 2回使うんで変数rhoに束縛しているが、普通にコードを書くときはそこまで 気をつかう必要はない。確かにGaucheのコンパイラはcommon subexpression eliminationをまだやってくれないけれど、通常のコンテキストでは(sqrt (- 1.0 r)) くらい何度計算してもたいしたことはない。 但し、今回このルーチンは10^8〜9回のオーダーで呼ばれるので、こういう微妙なところが 効いてくる場合がある。(効果は共通式の内容によるので、常にプロファイルせよ)
;; Macro define-simple-struct creates a bunch of functions and macros
;; to emulate a structure by a vector.
;; (This is a kludge. Future version of Gauche will have records as
;; efficient as this trick.)
;;
;; NAME is a symbol to name the structure type. TAG is some value
;; (usually a symbol or an integer) to indicate the type of the
;; structure.
;;
;; (define-simple-struct <name> <tag> <constructor> [(<slot-spec>*)])
;;
;; <constructor> : <symbol> | #f
;; <slot-spec> : <slot-name> | (<slot-name> [<init-value>])
;;
;; For each <slot-spec>, the following accessor/modifier are automatially
;; generated.
;;
;; NAME-SLOT - accessor (macro)
;; NAME-SLOT-set! - modifier (macro)
(define-macro (define-simple-struct name tag constructor . opts)
(let-optionals* opts ((slot-defs '()))
(define (make-constructor)
(let ((args (gensym))
(num-slots (length slot-defs))
(slot-names (map (lambda (s) (if (symbol? s) s (car s)))
slot-defs))
(init-vals (map (lambda (s) (if (symbol? s) #f (cadr s)))
slot-defs)))
`(define-macro (,constructor . ,args)
(match ,args
,@(let loop ((n 0)
(r '()))
(if (> n num-slots)
r
(let1 carg (take slot-names n)
(loop (+ n 1)
(cons
`(,carg
(list 'vector
,@(if tag `(',tag) '())
,@carg
,@(map (cut list 'quote <>)
(list-tail init-vals n))))
r)))
))))
))
`(begin
,@(if constructor
`(,(make-constructor))
'())
,@(let loop ((s slot-defs) (i (if tag 1 0)) (r '()))
(if (null? s)
(reverse r)
(let* ((slot-name (if (pair? (car s)) (caar s) (car s)))
(acc (string->symbol #`",|name|-,|slot-name|"))
(mod (string->symbol #`",|name|-,|slot-name|-set!")))
(loop (cdr s)
(+ i 1)
(list*
`(define-macro (,acc obj)
`(vector-ref ,obj ,,i))
`(define-macro (,mod obj val)
`(vector-set! ,obj ,,i ,val))
r))))))
))
;; Represents intersection. One instance of intersection is allocated
;; per scanline and reused by all traces.
(define-simple-struct is 'intersection make-intersection
([t 1.0e+30] ; distance from the ray origin
[updater #f] ; a thunk to update intersection
[p (make-f64vector 3 0.0)] ; where the ray intersects
[n (make-f64vector 3 0.0)] ; surface normal
[col #f] ; color of the surface (pathtrace)
[emissiveCol #f] ; emissive color of the surface (pathtrace)
[tmp (make-f64vector 3 0.0)] ; temporary working area
; (ugly, but to reduce allocation)
))
(define (reset-intersection! is)
(is-t-set! is 1.0e+30)
(is-updater-set! is #f)
is)
レイとプリミティブの交点を保持する構造体の定義。
Gaucheで構造体めいたものを扱うにはオブジェクトシステムがあるのだけれど、 残念ながらスロットアクセスがあんまり速くない。ここでは仕方なく、 define-simple-structというマクロで、ベクタを構造体のかわりにしている。 非常にアドホックなマクロなので、あまり一般的に使うことはお勧めしない。 hygienicじゃないので色々落とし穴があるし。
予定としては、このマクロを使うのと同程度に高速なレコードタイプを 組み込みでサポートするようにするつもり。srfi-9とr6rs recordのAPIを その上に載せる。なので特に急ぎでない人はあせってdefine-simple-struct をコピーするよりもそちらを待たれるのがよろしかろう。
;;-------------------------------------------------
;; Primitives.
;; A primitive returns a procedure that calculates intersection.
(define (Sphere center radius col emissiveCol)
(define (update-intersection! isect ray)
(let ([pp (is-p isect)]
[nn (is-n isect)])
(vset! pp (dir ray))
(vmul! pp (is-t isect))
(vadd! pp (org ray))
(vset! nn pp)
(vsub! nn center)
(vnormalize! nn))
(is-col-set! isect col)
(is-emissiveCol-set! isect emissiveCol))
(define radius^2 (* radius radius))
(lambda (isect ray)
(let* ([rs (rlet1 v (is-tmp isect)
(vset! v (org ray)) (vsub! v center))]
[B (vdot rs (dir ray))]
[C (- (vdot rs rs) radius^2)]
[D (- (* B B) C)])
(when (> D 0.0)
(let1 t (- 0.0 B (sqrt D))
(when (< 0.0 t (is-t isect))
(is-t-set! isect t)
(is-updater-set! isect update-intersection!)))))))
(define (Plane p n col emissiveCol)
(define (update-intersection! isect ray)
(vset! (is-n isect) n)
(vset! (is-p isect) (dir ray))
(vmul! (is-p isect) (is-t isect))
(vadd! (is-p isect) (org ray))
(is-col-set! isect col)
(is-emissiveCol-set! isect emissiveCol))
(define p.n (vdot p n))
(lambda (isect ray)
(let1 v (vdot (dir ray) n)
(when (> (abs v) 1.0e-6)
(let1 t (/ (- (- (vdot (org ray) n) p.n)) v)
(when (< 0.0 t (is-t isect))
(is-t-set! isect t)
(is-updater-set! isect update-intersection!)))))))
;; Find closest intersection of ray and objects.
;; Updates isect. Returns true if ray intersects, false otherwise.
(define (find-is! isect ray objs)
(reset-intersection! isect)
(let loop ((objs objs))
(unless (null? objs)
((car objs) isect ray)
(loop (cdr objs))))
(cond [(is-updater isect) => (lambda (u) (u isect ray) #t)]
[else #f]))
ジオメトリプリミティブの定義。.mjtさんの実装に倣い、ジオメトリプリミティブの 実体はレイを受け取って交点を計算するクロージャとした。
交点をアロケートしていると重いので、受け取った交点オブジェクトisectを 副作用で変更する。交点計算のところも副作用をばしばし使って、 f64vectorを一度もアロケートすることなく計算できるようにしている。 読みにくくなるからあまりお勧めしない。
vmul!, vadd!などが並ぶコードはアセンブラのようでもある。 もしこんなコードを大量に書くはめになったら、ベクタ演算式を こういうプリミティブに展開するようなマクロを書いた方が良いだろう。
あと、.mjtさんの実装ではクロージャ内で交点位置の法線を計算していたけれど、 それが必要なのは最も近い交点だけなので、先に一番近い交点を持つオブジェクトを 選択してからそれについてだけ位置と法線を計算するようにしている。ちょっとずるい。
find-is!はプリミティブのリストとレイから一番近い交点を計算して、 isectを破壊的に変更する。ループのところは、普通なら (for-each (cut <> isect ray) objs) の一行で済ませるとこなんだけど、 あいにくこう書くと今のGaucheでは実行時クロージャアロケーションが 起きてしまうのでわざとループにしてる。 これも、普段、99%のケースではfor-eachで書くべき。 ボトルネックとわかった時点でループに書き換える。
;;-------------------------------------------------
;; Main entry
;;
(define-constant IMAGE_WIDTH 256)
(define-constant IMAGE_HEIGHT 256)
(define-constant NSUBSAMPLES 2)
(define-constant NPATH_SAMPLES 128)
(define-constant NAO_SAMPLES 8)
(define-constant MAX_TRACE_DEPTH 16)
(define (main args)
(let-args (cdr args) ([type "t=y" 'ao]
[nprocs "p=i" 1]
[else => (lambda _ (usage))])
(call-with-output-file "out.raw"
(cute (if (> nprocs 1) do-parallel do-serial)
(case type
[(ao) ambientOcclusion]
[(pt) pathTrace]
[else (exit 1 "invalid type (must be either ao or pt): ~a" type)])
nprocs <>))
(print "DONE")
0))
(define (usage)
(print "Usage: aobench [-t ao|pt] [-p nprocs]")
(exit 0))
(define (do-serial sampler _ out)
(dotimes [y IMAGE_HEIGHT]
(let1 scanline (render IMAGE_WIDTH IMAGE_HEIGHT y NSUBSAMPLES sampler)
(print "LINE"y" - " (u8vector-length scanline))
(write-block scanline out))))
;; naive parallelization by static partitioning
(define (do-parallel sampler nproc out)
(let1 lines (make-vector IMAGE_HEIGHT #f)
(define (calc y)
(print "LINE"y)
(set! (ref lines y)
(render IMAGE_WIDTH IMAGE_HEIGHT y NSUBSAMPLES sampler)))
(define (task off)
(do-ec (: y off IMAGE_HEIGHT nproc) (calc y)))
(for-each thread-join!
(list-ec (: off nproc)
(thread-start! (make-thread (cut task off)))))
(dotimes [y IMAGE_HEIGHT] (write-block (ref lines y) out))))
(define (render width height y nsamp sampler)
(define isect (make-intersection))
(define sum (vec 0.0 0.0 0.0))
(rlet1 img (make-u8vector (* width 4))
(dotimes [i width]
(vset! sum 0 0 0)
(dotimes [u nsamp]
(dotimes [v nsamp]
(let ([px (/ (+ i (- (/. u nsamp) (/. width 2))) (/. width 2))]
[py (- (/ (+ y (- (/. v nsamp) (/. height 2))) (/. height 2)))])
(vadd! sum (sampler (new-ray (vec 0.0 0.0 0.0)
(vnormalize! (vec px py -1.0)))
isect)))))
(u8vector-set! img (+ (* i 4) 0) (clampu8 (/ (vx sum) (* nsamp nsamp))))
(u8vector-set! img (+ (* i 4) 1) (clampu8 (/ (vy sum) (* nsamp nsamp))))
(u8vector-set! img (+ (* i 4) 2) (clampu8 (/ (vz sum) (* nsamp nsamp)))))))
main関数と共通部分。このへんは実行される回数がたいしたことないので、 特に工夫する必要はない。 x,yでループするような構造はdotimesとかsrfi-42使った方が見やすいと思う。
;;-------------------------------------------------
;; Ambient occlusion
;;
(define *objects-ao*
(list
(Sphere (vec -2.0 0.0 -3.5) 0.5 #f #f)
(Sphere (vec -0.5 0.0 -3.0) 0.5 #f #f)
(Sphere (vec 1.0 0.0 -2.2) 0.5 #f #f)
(Plane (vec 0.0 -0.5 0.0) (vec 0.0 1.0 0.0) #f #f)))
(define (ambientOcclusion ray isect)
(if (not (find-is! isect ray *objects-ao*))
0.0
(let ([eps 0.00001]
[ntheta NAO_SAMPLES]
[nphi NAO_SAMPLES])
;; set ray origin
(let1 p (org ray)
(vset! p (is-n isect)) (vmul! p eps) (vadd! p (is-p isect)))
(let loop ([j 0] [i 0] [occ 0.0])
(cond
[(= j ntheta) (/ (- (* ntheta nphi) occ) (* ntheta nphi))]
[(= i nphi) (loop (+ j 1) 0 occ)]
[else
(receive (x y z) (random-direction)
(receive (b0 b1 b2) (orthoBasis (is-n isect))
;; update ray direction. be careful not to overwrite b2,
;; which is the same vector as (is-n isect).
(let1 d (dir ray)
(vmul! b0 x) (vmul! b1 y) (vset! d b2) (vmul! d z)
(vadd! d b0) (vadd! d b1))
(let1 hit (find-is! isect ray *objects-ao*)
(loop j (+ i 1) (if hit (+ occ 1.0) occ)))))])))))
ambient occlusionの計算。渡されたrayとisectを使いまわしてる。 ベクタ演算も副作用使いまくり。 チューニング中に、順序を間違えて変なバグを仕込み、取るのに苦労した。 このようなコードは決してお勧めしない。が、背に腹が替えられない時もある。
;;-------------------------------------------------
;; Path tracer
;;
(define *objects-pt*
(list
(Sphere (vec -1.05 0.0 -2.0) 0.5 (vec 0.75 0.0 0.0) (vec 0.0 0.0 0.0))
(Sphere (vec 0.0 0.0 -2.0) 0.5 (vec 1.0 1.0 1.0) (vec 1.0 1.0 1.0))
(Sphere (vec 1.05 0.0 -2.0) 0.5 (vec 0.0 0.0 1.0) (vec 0.0 0.0 0.0))
(Plane (vec 0.0 -0.5 0.0) (vec 0.0 1.0 0.0) (vec 1.0 1.0 1.0) (vec 0.0 0.0 0.0))))
(define (trace! ray isect depth)
(if (not (find-is! isect ray *objects-pt*))
(vec 0.7 0.7 0.7)
(receive (x y z) (random-direction)
(receive (b0 b1 b2) (orthoBasis (is-n isect))
(let ([d (dir ray)] [o (org ray)])
;; update ray direction. be careful not to overwrite b2,
;; which is the same vector as (is-n isect).
(vmul! b0 x) (vmul! b1 y) (vset! d b2) (vmul! d z)
(vadd! d b0) (vadd! d b1)
;; update ray origin
(vset! o d) (vmul! o 0.00001) (vadd! o (is-p isect)))
(let ([fr (is-col isect)]
[Le (is-emissiveCol isect)])
(rlet1 col (if (= depth (- MAX_TRACE_DEPTH 1))
(vec 0.0 0.0 0.0)
(trace! ray isect (+ depth 1)))
(vmul! col fr) (vadd! col Le)))))))
(define (pathTrace ray isect)
(rlet1 col (vec 0.0 0.0 0.0)
(dotimes [i NPATH_SAMPLES]
(vadd! col (trace! (copy-ray ray) isect 0)))
(vmul! col (/. NPATH_SAMPLES))))
パストレーサ特有部分。同様にrayとisectは使いまわし。さらに、 各traceが返す色も破壊的に変更している。 既定ケース(再帰しないケース)で定数ベクタ'#f64(0.7 0.7 0.7)ではなくわざわざ(vec 0.7 0.7 0.7)のように 新たなベクタを作っているのはそのため。
Tag: CG
enami (2009/02/27 04:11:27):
shiro (2009/03/06 00:10:05):
enami (2009/03/08 01:06:22):