Syoyo Fujita's Blog

raytracing monte carlo

Month: January, 2008

Quasi-Monte Carlo light transport simulation by efficient ray tracing

(from ompf.org/forum)

Quite interesting and fresh thesis was published, from the author of BIH.

Quasi-Monte Carlo light transport simulation by efficient ray tracing
http://vts.uni-ulm.de/query/longview.meta.asp?document_id=6265

I didn’t read it yet, I have to read it soon.

# but nowdays, there are a lot of papers I have to read.
# 1 paper per day is absolutely not enough. Dozens of papers per day may be needed to catch up the trend in the literature.

Certifying floating-point implementations using Gappa

New gappa preprint was posted.

Florent De Dinechin, Christoph Quirin Lauter, Guillaume Melquiond
Certifying floating-point implementations using Gappa
http://arxiv.org/abs/0801.0523

I’m planning to combine gappa with my MUDA language for verified and accurate floating point computation.
But gappa is designed for scalar fp code, while MUDA is for vector fp code.
I have to find a way to solve this scalar-vector problem 😉

WaveScript(Regiment) language.

(from Lambda the Ultimate)

Regiment(a.k.a WaveScript) homepage
http://regiment.us/regiment_webpage/index.html

WaveScript is the programming language used to develop WaveScope applications. It is a high-level, functional, stream-processing language that aims to deliver uncompromising performance. WaveScript programs execute in parallel on multiple cores, or distributed across a network. Its compiler uses aggressive partial evaluation techniques to remove abstractions and reduce the source program to a graph of stream operators.

functional, stream-process, partial evaluation and multi-core.
Hmm, it would be worth to incorporate idea of WaveScript language into global illumination language(shader + sampler + raytracing + etc) I’m thinking now. I must investigate it.

I know there is already some advanced language(or shader system) such like RTSL(RayTracing Shading Language) and lightspeed,
but those are still incomplete as a global illumination language.

For example, I think we need below features for a global illumination language, and no one find a way how to do it.
– automatic differentiation computation(not in image-space)
– find & exploit time-spatio coherence
– clustering incoherent rays
– automatic importance sampling for user-defined brdf

See also:

Faust, Signal Processing Language

一人で 7600 億の損失を出した男

いやー、「マネー革命」が放映された頃の 90 年代後半の金融危機よりも
がぜん面白くなってきました。いや、もうすでに面白い段階だったので、さらに加速されているいう感じ。
(時は加速するというのは本当だったのか!)

仏政府は、同国金融史上最大といわれる巨額不正の全容解明を求めている
http://www.afpbb.com/middle/790

一人のトレーダーが出した損失としては、
日経デリバティブで 1600億円 の損失を出してベアリングス銀行を破綻に追い込んだニック・リーソン、
住友商事のトレーダーで、銅取引で 3000 億の損失をだした浜中泰夫、
が有名どころですが、
この記録をぶっちぎって、一人のトレーダーによる損失では過去最大の記録が打ち立てられました。

さすがソジェン!
俺たちにできないことを平然とやってのける!
(ベアリングスの教訓が生きていなかったとは!)

次はどこかな〜

任天堂、おまえもか!

任天堂、3Q 決算発表で、創業以来初の売り上げ 1 兆円越え、経常利益が前期比で 2 倍と過去最高になるも、
為替差損による見通し据え置きが失望されて(?)、1/28日はストップ安。

http://quote.yahoo.co.jp/q?s=7974.o&d=c&k=c3&z=m&h=on

29 日の一般教書演説で、大統領は「米経済は不透明な時期にある」と、景気減速を認めた。
http://www.yomiuri.co.jp/editorial/news/20080129-OYT1T00854.htm

やっとかよ…

Fast Barrier for x86 Platforms

自動最適化で Spiral の資料をあさることが多くなっていますが、
ひとつ面白そうなネタを見つけました。

x86 での最速なバリア同期処理は何か?というもの。

Fast Barrier for x86 Platforms
http://www.spiral.net/software/barrier.html

まず、驚きだったのが、pthread のバリア同期関数の遅さ。
2×2 cores では 27000 cycles もかかっている。
これは CPU の周波数が 2.7 GHz なら、0.1 msec かかるということになる。ひょえ〜〜

lock + アトミック加算でのバリア同期だと、数千サイクル程度に減る。

どれもそうだけど、2×2 cores(2 dies) では、FSB をまたぐ必要があるわけで 2 cores(1 die)
より遅くなるのは当然。

以外と健闘しているのが、Icc + OpenMP によるバリア同期。
ただこれはコンパイラ次第によって変わるだろうから、あまり参考にはならなそう。

で、最後に一番早いと提案されているのが、なんか独自の?配列を使った方法。
lock + アトミック加算より 2 倍くらい早くて、アセンブラを使わなくてよいという利点がある。
2 cores の時で 200 サイクルくらいなので、一回メインメモリにアクセスするのと同じくらいの
コストで同期がとれるということになる。

いずれにせよ、ここらへんで最適化するよりは、
まずは同期の回数が少ない並列プログラムを書くようにするほうがよいでしょうね。

floating point. part 1 : éžæ­£è¦åŒ–数の恐怖

浮動小数点演算(IEEE 754)のことを調べていると、
いろいろと今まで lucille などでの浮動小数点計算コードが実はダメダメで、
一から考え直す必要があることが分かってきました。
(MUDA 用の libm ももっと fp の仕組みを知ってから書き始める必要がありそう)

そんなわけで、グラフィクス野郎が知っておくべき浮動小数点演算シリーズ.

第一回は、非正規化数の恐怖です。

非正規化数(denormal number, subnormal number)は、演算がアンダーフローしたときに、
いきなりアンダーフローの結果がゼロになると困るような計算のために IEEE 754 に導入されたようです。
しかし、非正規化数は、基本的にはソフトウェアで処理されるので、実はとてもパフォーマンスに影響があることが知られています。
(x86 OS では、基本的に非正規化数の扱いはデフォルトで有効になっている)

アンダーフローはそんなに滅多に起こるものではない。
しかし一度起きてしまうと、アンダーフローは一気にやってくる。
一回でもアンダーフローを生じさせるような計算は、通常何度もアンダーフローを発生させるようになるからだ。
そのため、プログラムは時々アンダーフローの嵐に見舞われる。これがプログラムの実行を遅くするのだ。

W. Kahan
IEEE 754R minutes from September 19, 2002

アンダーフロー(非正規化数)は通常あんまり起こらないが、
起こってしまうとずっとそれは生じてしまう、ということになる。

つまりは、今の金融市場と同じですね。
大きく指数が下落すると、ヘッジ(or 損切り)のための売りがでて、
その売りがさらに他の売りを誘発するという流動性の罠みたいな感じ。

では実際どれくらい、演算に非正規化数があると効率が悪くなるのだろう?
[1] によると、非正規化数の処理のコストは 1000~ サイクルである!

世の中は広い。なんとすでに非正規化数の演算がどれくらい効率が悪いかのベンチマークプログラムが存在する [2]。
(私が一日がかりで x86 の denormal の振る舞いや fp control register の仕様と
格闘しながら作ったプログラム [5] よりも遥かに出来がいいなぁ)

私の環境(gcc 4.0.1, core 2 intel mac) では約 25 くらいのベンチ結果となった。
(これは [5] を実行したときのパフォーマンス比ともほぼ一致する)

[2] のベンチは [3,4] で解説されているように、至極単純である。


slow:
a <- [1, 0, 0, 0, 0...]
for i <- 1 to ITER
  for n <- 3 to SIZE
    a[n] <- (a[n] + a[n-1] + a[n-2]) / 3

fast:
tiny <- small normalized value
a <- [1, tiny, tiny, tiny, tiny...]
for i <- 1 to ITER
  for n <- 3 to SIZE
    a[n] <- (a[n] + a[n-1] + a[n-2]) / 3

time = slow / fast

という感じで、どんどん平均を取っていくコードを繰り返して実行する。
最初の slow ループは、a[0] だけ 1 で、残りの a[] の配列の要素の値は 0 である。
これはつまり、n が増えるにつれて、どんどんと a[] の配列の値は小さくなっていっていき、ある時点で
計算された a[n] は非正規化数の領域へと突入する。

それに対して、 fast のループでは 0 を小さな正規化数である tiny に置き換えることで、
計算結果が非正規化数にならないように防いでいる。

パフォーマンスの差は、この slow と fast のループの処理時間から割り出される。

[2, 3, 4] によれば、大きな傾向として、最新の CPU アーキティクチャほど
非正規化数の処理パフォーマンスが悪いという結果が読み取れる。

対応策

プログラムがアンダーフローを引き起こすか、演算結果が非正規化数になるかどうかは、
Prof. Kahn の引用のように、これはアルゴリズムに起因する部分が大きいと思う
(たとえば行列計算や, 物理計算での time integration などでは生じやすいと思う)。

対応としては、演算結果がアンダーフロー(非正規化数)になるかどうか演算のコードをチェックして、
回避策を探すのがよい。

– 非正規化数による漸進的アンダーフローがアルゴリズムに
それほど影響を与えないのであれば、非正規化数を無効にする
(x86 では DAZ(Denormals are zero), FZ(Flash to zero) flag を設定)
基本的にグラフィックス野郎にとっては非正規化数の精度は必要ないだろうから、
とりあえず非正規化数は off にしておくのがよいだろう。

– 計算コードがアンダーフローしないように入力の値の定義域を制限したり、
アルゴリズムの見直しを行う
(値域、定義域のチェックには gappa などが使えるだろう)

ちなみに、ゲーム機などのプロセッサでは非正規化数は通常サポートされていません。

おわりに

というわけで、知らずのうちにあなたのプログラムが非正規化数の演算を引き起こしていて、
そのせいで実は本来の 1/10 のパフォーマンスしか出ていなかったかもしれない、
という身の震えるような(?) 浮動小数点計算の落とし穴を紹介しました。

[1] SSE Performance Programming
http://developer.apple.com/hardwaredrivers/ve/sse.html

[2] Subnormal Floating-Point Degradation in Parallel Applications
http://charm.cs.uiuc.edu/subnormal/

[3] Performance Degradation in the Presence of Subnormal Floating-Point Values
Hari Govind, Isaac Dooley, Michael Breitenfeld, Orion Lawlor and Laxmikant Kale
http://research.ihost.com/osihpa/
[PPT]

[4] Quantifying the Interference Caused by Subnormal Floating-Point Values.
Isaac Dooley, Laxmikant Kale
http://osihpa.eng.utep.edu/2006/DooleySubnormal06.pdf

[5] http://lucille.svn.sourceforge.net/viewvc/lucille/angelina/fp/denormal/

[6] x86 では、演算結果が非正規化数になる時点でパフォーマンスダウンが発生する模様。
非正規化数になったかどうかの判定に、浮動小数点例外が使える。
しかし denormal 例外は、オペランドに denormal が与えられたときでないと発生しない。
つまり非正規化数になったかどうかを例外で捉えるには denormal 例外だけではだめで、
underflow 例外を使うことになる。
しかし、underflow したら非正規化数の領域をすっ飛ばしてすぐに 0 になるケースでは、
パフォーマンスダウンは発生しない。つまり underflow 例外を補足したからといって、
それが denormal が発生したとは言えないことになる。ややこしい。

Recent advances on motion blur rendering

Last week, I’ve talked/discussed a lot with renderists over the world about the offline renderer.

After the discussion, I think we’ve confirmed that there is still a lot of things to do in (global illumination) rendering.

One of it is fast & efficient computation of motion blur, this still has not been solved fully in the context of RenderMan style renderer, as well as ray tracing renderer.

After I’ve researched recent advances on motion blurring, I found interesting paper.

Eulerian Motion Blur
Doyub Kim and Hyeong-Seok Ko,
In Eurographics Workshop on Natural Phenomena, 2007,

http://www.doyub.com/research/

eulerian_mb.png
(image from Eulerian Motion Blur)

Almost all previous motion blur technique focuses on Lagrangian motion(e.g. moving objects).
They propose new Eulerian motion blur techniques.
(Some people may imagine this approach is similar to image-based motion blur technique)

Natural phenomenon & fine geometries(fluids, clouds, gases, furs, etc) become dominant part in the (production) graphics, thus I think such a Eulerian motion blur techniques will be developed & used intensively near the future.

Here is another recent papers about motion blurring.

Zheng, Yuanfang ; Köstler, Harald ; Thürey, Nils ; Rüde, Ulrich:
Enhanced Motion Blur Calculation with Optical Flow .
In: Aka GmbH (Hrsg.) : RWTH Aachen (Veranst.) : Proceedings of Vision, Modeling and Visualization 2006 (Vision, Modeling and Visualization Aachen 22. – 24. Nov 2006). Aachen : IOS Press, 2006, S. 253–260. – ISBN Erscheinungsjahr.
(see my previous post on this paper)

Frank Dachille, Arie Kaufman,
High-Degree Temporal Antialiasing
Proceedings of Computer Animation 2000, pp. 49-54, 2000.
http://citeseer.ist.psu.edu/455657.html

highdegree_mb.png
(image from High-Degree Temporal Antialiasing)

(my old post on this paper)

Pixmotor: A Pixel Motion Integrator
SIGGRAPH ’07 Technical Sketch, August 2007
Ivan Neulander
http://www.rhythm.com/~ivan/pubs.html

pixmotor.png
(image from Pixmotor)

Vlastimil Havran, Cyrille Damez, Karol Myszkowski and Hans-Peter Seidel
An Efficient Spatio-Temporal Architecture for Animation Rendering
Proceedings of the 13th Eurographics workshop on Rendering, Leuven, Belgium, Pages: 106 – 117, 2003
http://www.mpi-inf.mpg.de/resources/anim/EGSR03/

spatio_temporal_arch.png
(image from An Efficient Spatio-Temporal Architecture for Animation Rendering)

(my old post on this paper)

[Ja]

まだまだ、グラフィックスの研究はやるべきことは一杯あるなー、と実感。
まあモーションブラーについては前にも指摘したと思いますが、今一度。

Parsing floating point format in Haskell + Parsec

When I tried to port simdmath library code to MUDA language,
I’ve faced a problem on parsing floating point format in Haskell.

Basic read function provided by Haskell Prelude cannot parse
C style floating point format fully. For example,


$ ghci
Prelude> read "1.0" :: Double
1.0
Prelude> read "1.0f" :: Double
*** Exception: Prelude.read: no parse
Prelude> read "1.0e-3f" :: Double
*** Exception: Prelude.read: no parse
Prelude> read "-1.0e-3" :: Double
-1.0e-3
Prelude> read "+1.0e-3" :: Double
*** Exception: Prelude.read: no parse
Prelude> read "1.0e+3" :: Double
1000.0

Hmm, explicitly specifying sign and “f” seems not supported by read function.

So I’ve searched a lot(really a lot, it was a waste of my time…)
and found Parsec provides better support for parsing floating point format.


Text.ParserCombinators.Parsec.Token
provides floating point parser float.


$ cat Parser.hs
import Text.ParserCombinators.Parsec
import qualified Text.ParserCombinators.Parsec.Token as P
import Text.ParserCombinators.Parsec.Language (emptyDef)

parseFloat :: GenParser Char st Double
parseFloat = do x <- P.float $ P.makeTokenParser emptyDef
                return x


main :: IO ()
main = do
            parseTest parseFloat "1.0f"
            parseTest parseFloat "1.0e-2f"
            parseTest parseFloat "-1.0e-2f"
            parseTest parseFloat "+1.0e-2f"
            parseTest parseFloat "12.3E+3"
            parseTest parseFloat ".1"

$ runhaskell Parser.hs

1.0
1.0e-2
parse error at (line 1, column 1):
unexpected "-"
expecting float
parse error at (line 1, column 1):
unexpected "+"
expecting float
12300.0
parse error at (line 1, column 1):
unexpected "."
expecting float

float parser still lacks support for handling sign character.
Thus, I’ve added extra code which parses sign.

Final solution is,



$ cat Parser.hs

import Text.ParserCombinators.Parsec
import qualified Text.ParserCombinators.Parsec.Token as P
import Text.ParserCombinators.Parsec.Language (emptyDef)

parseFloat :: GenParser Char st Double
parseFloat = do sign <- option 1 (do s <- oneOf "+-"
                                     return $ if s == '-' then-1.0 else 1.0)
                x    <- P.float $ P.makeTokenParser emptyDef

                return $ sign * x


main :: IO ()
main = do
            parseTest parseFloat "1.0f"
            parseTest parseFloat "1.0e-2f"
            parseTest parseFloat "-1.0e-2f"
            parseTest parseFloat "+1.0e-2f"
            parseTest parseFloat "12.3E+3"
            parseTest parseFloat ".1"

$ runhaskell Parser.hs

1.0
1.0e-2
-1.0e-2
1.0e-2
12300.0
parse error at (line 1, column 1):
unexpected "."
expecting float


Parsing “.1” still fails,
but MUDA doesn’t support such a format for floating point string
since it is something a confused format, by the low of my sense of beauty 😉

AAPL, お前もか!

Apple, 四半期決算で過去最高の売り上げと利益を達成するも、
景気減速による来期見通しが悲観され、時間外で株価 12% down とボロクソに売られた.

http://itpro.nikkeibp.co.jp/article/NEWS/20080123/291756/

現在 -17% まで売り込まれている.

http://finance.google.com/finance?q=NASDAQ:AAPL

USD/JPY は一時 105 円を割った… おそろしや…
円は世界最弱通貨のはずなのだが、他通貨(他国家)の健康状態があまりにもひどいので、消去法で買われている ^^)

usdjpy105.png