Haskellのプログラムを高速化するための７つの方法

以前のエントリでHaskellのプログラムが遅いというのが不満だったので、Haskell-cafeにメールして聞いてみました。

いくつか高速化するためのテクニックを教えてもらったので、まとめてみます。

-O2オプション

ghcでコンパイルするときに-Oを付けると最適化してコンパイルされますが、それを-O2オプションに変えるとさらに最適化されます。

-fexcess-precision

ghcでコンパイルするときに、-fexcess-precisionオプションを付けてnumber-crunchingをすると、処理が速くなります。
★number-crunchingって何だろう？

-fvia-C,-optc-O3

これらのオプションを付けると、HaskellのコードをいったんCに変換し、それを-O3オプション付きでコンパイルします。２回コンパイルするのでコンパイルにかかる時間は増えますが、さらなる高速化を期待できます。

-fvia-Cと-optc-O3を付けることで、約10%高速化できるようです。

UNPACKプラグマと正格フラグ

データ構造の要素にUNPACKプラグマと正格フラグ(!)を付けることで、計算の高速化を期待できます。

たとえば、Particleの定義を、

data Particle = Particle { pos, vel :: {-# UNPACK #-} !Vector3,
mass, rho, prs :: {-# UNPACK #-} !Scalar }

のようにします。

-fexcess-precisionオプションとUNPACKプラグマ、正格フラグを使うことで、約30%高速化できるようです。

★!フラグを付けることで正格評価されるという理解でいいのかな？
★UNPACKプラグマを付けるとどうなるんだろう？非ボックス化というやつ？抽象化せず直で要素を扱うのかな？

ByteString

HaskellのStringは処理が遅いらしいので、代わりにByteStringを使うことで高速化を望めます。

コンパイルオプションの最適化

Haskellに限ったテクニックではありませんが、Acoveaのようなコンパイルオプション最適化ツールを使うことで、コンパイルオプションを最適化して高速化できます。Acoveaは遺伝的アルゴリズムでコンパイルオプションを最適化するようです。

DPH

DPHを使うことでも高速化できるようです。
★DPHってなんだろう？

参考になれば。
-

遅延ストリームによってシミュレーションと出力処理を分離

このエントリで紹介したHaskellのコードでは、遅延ストリームという手法を使っています。

このエントリでは、遅延ストリームとそれを使うメリットについて書きます。

遅延ストリームとは？

遅延ストリームとは、遅延評価される無限リストです。

使っているのは、

output $ take (read (n!!0)) $ simulation ps0

の部分です。

simulation関数の型は

[Particles]->[[Particles]]

となっていて、粒子の初期状態を引数にとり、各ステップでの粒子の状態をリストで返します。

そのリストから、コマンドライン引数で与えたステップ数だけtake関数でとりだし、output関数（アクション）でファイルに出力しています。

このように書くと、取り出すステップ数だけ先に計算して、すべて計算し終わってからまとめて出力されることになりそうですが、Haskellの場合は遅延評価されるので、１ステップずつ評価されてファイルに出力されるのです。

遅延ストリームを使うメリットは？

このように遅延ストリームを使うメリットは、関数のモジュール性を高めて、シミュレーションと出力を完全に分離できることです。

simulation関数とoutput関数はまったく独立していて、たとえば、output関数のみをdisplay関数に差し替えて画面に出力する、ということができます。

display $ simulation ps0

逆に、この処理を遅延評価でない言語、たとえばC++で書くと、

for ( int i = 0 ; i < n; i++ )
  {
    output_particles( p_ps, i );
    simulation( p_ps );
  }

のように、forループの中にシミュレーションと出力を書くことになり、お互いが強く結びついてしまいます。遅延ストリームをなら、このように強く結びつくことを避けられます。

以上のように、遅延ストリームを使うと、関数のモジュール性を高めることができるのです。

OCamlで遅延ストリームは使える？

OCamlも言語に遅延評価の仕組みが含まれているので、同様に遅延ストリームを使うことができます。

OCamlでの遅延ストリームは、このエントリにあるコードから抜粋すると、

open Lazy;;

type 'a stream = Nils
               | Cons of 'a * 'a stream lazy_t;;

let rec stream_take (n : int) (s : 'a stream) : 'a stream =
  match (n,s) with
  | (_,Nils)       -> Nils
  | (0,Cons(x,_ )) -> Nils
  | (n,Cons(x,xs)) -> Cons(x, lazy(stream_take (n-1) (force xs)));;

let rec simulation ps = 
  let ps'  = calc_amount ps (mk_neighbor_map ps) in
  let ps'' = advance ps' (mk_neighbor_map ps')
  in Cons(ps, lazy(simulation ps''));;

let main =
  let n = int_of_string (Sys.argv.(1))
  in print_number_of ps0;
     stream_output (stream_take n (simulation ps0));;

のように書くことができます。

'a stream型のcarを'a型、cdrを'a stream lazy_t型とし、cdrを順番にforceしていくことで遅延ストリームを実現しています。

Haskell版と同様に、粒子の状態を１ステップずつ評価してファイルに出力します。

C++では？

C++では言語に遅延評価の仕組みがありません。なので、遅延ストリームをそのまま使うことはできません。

ただ、シミュレーションと出力の分離という意味では、関数オブジェクトを使うことで、同等のことは一応可能です。
-

Haskellのコードを速くするにはどうしたらいいんだろう？

前のエントリでSPHのコードをC++とHaskellとOCamlで書いて比較したんだけど、あまりにHaskell版が遅いなぁ。

C++に対して約15倍も時間がかかってしまう。プロファイラを使ってチューニングしてもこの速度差。せめてC++の３倍くらいで済んで欲しい。

2009/11/5追記：計りなおすと、Haskell版はC++版の約５倍ですんでました。。

どうしたらもっと速くなるんだろう？

正格評価したらいいの？どうやるの？Data.Vectorモジュールには!フラグが付いているんだけど、それじゃだめなの？

ステップの更新を非破壊的にやっているから遅い？けどOCamlはそれでもC++の２倍で済んでるのに。

純粋な科学技術計算だから、どのみちすべての式を計算する必要があるし、その部分では遅延評価のメリットはないよなぁ。だから遅いのかなぁ？けど、ただ計算するだけだから遅延評価も正格評価も関係ないようにも思うんだけど。（遅延評価の方がメモリは使いそうな気がする）

やりたいこととしては、

思考をそのまま落とせるHaskellでプロトタイプを作る→速度の速いC++で書き直す

というフローを組みたいんだけど、それにしても遅すぎる。せめてC++の３倍くらいで済んで欲しいなぁ。あと、Haskellの書きやすさを損なうような方法はとりたくない。（けど、レンダリングの時間も含めて考えたら、シミュレーションは速くても遅くても関係ない、というケースはあるなぁ）

速くするにはどうしたらいいんだろう？何を調べたらいいんだろう？
-

【Emacs】行ジャンプするには

Emacsで行ジャンプするには、

M-x goto-line 行番号

とする。

よく使う機能なので、.emacsに

(global-set-key "\C-x\C-g" 'goto-line)

と書いておくと\C-x \C-gで行ジャンプできるのですごく便利。
-

Haskell、OCamlでSPH法

このエントリに載せていたC++のコードを、HaskellとOCamlで書いてみました。（あと、C++のコードも整理しました）

コードを置いておきます。

• Haskell版(sph.hs)
• OCaml版(sph.ml)
• C++版(sph.cpp)

Haskell版(sph.hs)の実行方法は、

$ ghc --make -O sph.hs
$ ./sph 300

です。./sphの引数は計算するステップ数です。上記のsph.zipにはData/Vector.hsというベクトル演算のためのモジュールを入れてあります。

同様に、OCaml版(sph.ml)の実行方法は、

$ ocamlopt -pp 'camlp4o -parser Camlp4ListComprehension' -o sph sph.ml
$ ./sph 300

です。-ppオプションは、プリプロセッサにCamlp4を使ってOCamlでリスト内包表記を使えるようにするためのものです。

もう１つC++版(sph.cpp)の実行方法は、

$ g++ -O2 -o sph sph.cpp
$ ./sph 300

です。以前はboostを使ってみていたのですが、boostを用意する手間やコンパイルの遅さを考えるとあまり意味がないので、使わないようにしました。（あと、このコードを書く前にSTLの総称プログラミングのスタイルでも書いてみたのですが、かえってわかりづらくなったのでやめました）

この３つの実行時間ですが、

▼C++
$ time ./sph 30
real 0m6.367s
user 0m4.916s
sys 0m0.157s

▼OCaml
$ time ./sph 30
real 0m11.718s
user 0m9.461s
sys 0m0.241s
→ cppの約２倍

▼Haskell
$ time ./sph 30
real 1m23.961s
user 1m2.786s
sys 0m2.169s
→ cppの約15倍

On PowerBook G4 1.5GHz (Mac OSX 10.4.11)

のようになりました。

OCamlがC++の約２倍、HaskellがC++の約15倍です。

Haskellがあまりに遅いのでOCamlで書いてみたというのがあるのですが、それにしても遅い。ここまで遅いものなのだろうか？

2009/11/5追記：計りなおすと、Haskell版はC++版の約５倍ですんでました…

それに対して、OCamlは本当に速いですね。いまのコード(sph.ml)は粒子を非破壊的に更新しているのですが、それでもC++の約2倍で実行できています。

次は、Haskell版とOCaml版で使っている遅延ストリームについて書きます。


■流体シミュレーションに関するエントリ
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
SPHによる巻き波のシミュレーション２
Haskell、OCamlでSPH法
このあとやりたいこと
カメラ位置を変えて流体をレンダリング
固液連成シミュレーション

■剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション

粒子法(SPH)のプログラム

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

-

ただ、君を愛してる

宮崎あおいつながりで、友達から薦められて「ただ、君を愛してる」をみた。

あー、もう泣きそう。胸が締めつけられる。

死とひきかえに誠人に愛される女性になろうとする、静流が健気すぎる。

SPHによる巻き波のシミュレーション第３回

前回のエントリで、巻き波のシミュレーションができた。今回は、それを綺麗にレンダリングしてみる。

まず、粒子でレンダリングしたものがこれ。




これを、マーチングキューブで陰関数曲面を抽出してレンダリングするとこうなる。





さらに、屈折や反射を考慮してレンダリングしたものがこれ。





いい感じ！次は何をしようかな？


■関連するエントリ
SPHによる巻き波のシミュレーション第１回
SPHによる巻き波のシミュレーション第２回

-

SPHによる巻き波のシミュレーション第２回

前回のエントリでは、可動壁によって波を起こすことで巻き波のシミュレーションをしようとしたが、うまく巻き波が起きなかった。

そこで、巻き波を起こすにはどうしたらいいのかいろいろと調べてみたところ、RadovitzkyとOritzの"Lagrangian finite element analytics of newtonian fluid flows."という論文を見つけた。

この論文は、有限要素法によってラグランジュ的に自由表面をシミュレーションする手法を提示している。その適用例の１つとして砕波を扱っており、そこに孤立波の理論解が示されていた。これを適用すれば波の岸に打ち寄せる波を再現することができそう。

孤立波の初期条件は、水面の高さをη、x方向の初速をu、y方向の初速をvとすると、



のようになる（x方向のオフセットを省略）。ここで、gは重力加速度で9.8[m/s^2]、Hとdはそれぞれ基準となる水面の高さと波の高さである。

また、tanhとsechは双曲線関数で、その定義は

である。

この初期条件を適用して、またシミュレーションしてみると、



お！巻き波が起きた！

シミュレーション結果を見ると、左から進んできた孤立波が、斜面にさしかかったところでせり上がって、巻き波になってくずれる様子が見てとれる。また、巻き波が打ち寄せたあとにさーっと引いていく様子も、実際の砂浜とよく似ている。

いい感じに巻き波が起きるようになったので、次は、このエントリのように綺麗にレンダリングしてみます。

■関連するエントリ
SPHによる巻き波のシミュレーション第１回
SPHによる巻き波のシミュレーション第３回
-

SPHによる巻き波のシミュレーション第１回

波の一種に巻き波というものがある。浅海域に進入した波が対称性を失って、前面が切り立った壁のようになり、頂部の水が前に飛び出す形で巻くように砕けるものをいう。

粒子法では、しぶきがあがるような激しい変形を表現することができる。そこで、巻き波をシミュレーションしたいと思った。



左側の境界条件を周期的に動かすことで波を起こしている。また、真ん中から右は斜面になっている。

ところが、結果をみてみるとうまくいかない…

波が減衰して斜面まで届かないのと、斜面の底部分で不自然に左向きの流れが起きているのが気になる。

波が減衰してしまう点については、容器の形や可動壁の動きを調整しようと思う。

斜面の底部分で不自然な流れが起きている点については、理由がよくわからない…境界条件にはペナルティ法を使っているが、そのせいで斜面方向に変に力がかかっているんだろうか？

もうすこし調整が必要そうだ。

■関連するエントリ
SPHによる巻き波のシミュレーション第２回
SPHによる巻き波のシミュレーション第３回
-

【論文】p = k ( ρ - ρ0 )でのρ0の意味

Mullerらの"Particle-Based Fluid Simulation for Interactive Applications"では、圧力pを密度ρから求める際に、定常密度ρ0を引いて

p = k(ρ - ρ0)

としている（kは定数）。ρ0を引く理由としてMullerらは、

... Since pressure forces depend on the gradient of the pressure field, the offset mathematically has not effect on pressure forces. However, the offset does influence the gradient of a field smoothed by SPH and makes the simulation numerically more stable. ...

と述べている。

なんのこっちゃ？

オフセット（ρ0）が数学的に圧力に影響を与えないといいつつ、それによってシミュレーションがより安定化するといっている。しかし、数学的に影響を与えないのであれば、安定化するもなにもないではないか。

そこで、引用元である、Desbrunの"Smoothed Particels: A new paradigm for animating highly deformable bodies"の中で、なぜρ0が導入されているのかを見てみた。

Desbrunらによると、

... In astrophysics applications, pressure forces were often combined with gravitational forces balancing the expansion phenomenon.
  In contrast, we would like to animate materials with constant density at rest. Consequently, the material sould exhibit some internal cohesion, resulting in attraction-repulsino forces as in the Lennard-Jones model. ...

とある。

なるほど。

もともとSPHは天体のシミュレーションで用いられていた手法なので、真空中での粒子の運動を対象としている。一方で、Desbrunらが行おうとしている大変形体のシミュレーションは、大気中での運動が対象となる。

大気中での運動なので、本来であれば、大変形体の周りに空気があって、その気圧と大変形体の圧力が釣り合っていないといけない。

それを論文でどう扱っているかというと、周りの空気も粒子としてシミュレーションに含めるのではなく、その代わりに、ρからρ0を引いて気圧の影響を補正した圧力を用いるという方法をとっている。つまり、定常状態での圧力を仮想的に0気圧となるようにオフセットをとって計算している。

つまり、定常状態での圧力を仮想的に0気圧とすることで、気圧と大変形体の圧力を釣り合わせているということになる。

実際に、ρ0を変えてシミュレーションをしてみた。

【ρ0=600】


【ρ0=300】


【ρ0=0】


確かに、ρ0が小さくなるにつれて、粒子が空気中に広がっていっている様子がわかる。ρ0=0に至っては、まさに真空中に拡散しているように見える結果が得られた。

結論
圧力pを密度ρから求める際に定常密度ρ0を引いてp = k (ρ - ρ0)とするのは、気圧の影響を加味して粒子が拡散してしまわないようにするためだった。SPHを宇宙から地球上に持ってくるわけですな。

■関連する記事
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…

ClickTaleの無料利用分を使い切ってしまった

おととい入れてみたClickTale。

２日で1週間分のレコーディング数を使い切ってしまった…。これじゃ使えないなぁ。

１人１人のユーザの動きを追える素晴らしいツールなんだけど、月$99は払えないや。

自分で作ろうか？？

ClickTaleの再生がAnalyticsに計測されてしまわないか？

このブログにClickTaleを入れてみた。

ClickTaleはユーザのマウスの動きをレコーディングしてくれるアクセス解析ツールで、レコーディングされたマウスの動きを再生できる。

再生機能で１つ気になったのが、再生したときのページ読み込みが解析対象ページに埋まっているGoogle Analyticsで計測されてしまわないかという点。

Analyticsで確認してみると、ClickTaleからのアクセスは計測されないみたいでよかった。

OneApp

MicrosoftがOneAppという携帯電話向けのアプリケーションフレームワークを発表した。

OneAppは、途上国の人々がTwitterやFacebookのようなサービスに携帯電話からアクセスできるようにするサービスで、プロセッサが非力でメモリが少なくてもJavaが動きさえすれば使うことができる。

大型Webアプリケーションを動かすために必要な処理の一部をMicrosoftのサーバがクラウドとして代行する仕組みで、デベロッパはJavaScriptとXMLによる比較的簡単な作業でアプリをOneAppに対応させることができる。

Microsoftが途上国向けにサービスを提供する理由には、人道的理由もあるが、途上国が潜在的に巨大な市場であるというビジネス上のものもあるだろう。

iTunesでムービーのサムネイルを変えるには

iTunesのムービーのサムネイルは好きなタイムラインのものに変えることができる。

▼変え方

1. サムネイルにしたいタイミングでムービーを停止する
2. ムービーを右クリック
3. 「ポスターフレームを設定」を選ぶ

これで好きなカットをサムネイルにできる♪

粒子法のシーンを２倍のサイズにしてみた

「粒子法のシーンを２倍のサイズにしてみたが…」で失敗したものをもう一度計算した。

【２倍のサイズ】


【元のサイズ】


初期配置や境界の形はまったく同じにして、xyz各方向に２倍している。一見同じように見えるけど、床のチェック模様の細かさが違うところに注目して欲しい。

２倍サイズにすると、より高いところからさらに多くの流体が落下するため、盛大なしぶきがあがった。上のレンダリング結果だと黒くなってしまってわかりづらいが、陰関数曲面をそのままレンダリングするとその様子がよくわかる。

粒子法(SPH)のプログラム

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--

【論文】Particle-Based Fluid Simulation for Interactive Applications

Matthias Mullerらの「Particle-Based Fluid Simulation for Interactive Applications」を読んだ。

この論文では、SPHをベースに、自由表面を持つ流体をインタラクティブにシミュレーションしレンダリングする方法を提案している。

この論文におけるMullerらの貢献は、

• Desbrunらの大変形体向けの手法を流体に拡張
• インタラクティビティの実現

の２点。

Desbrunの手法を流体に拡張するにあたり、Mullerらは流体が受ける力をナビエ=ストークス方程式から直接導出するようにし、さらに、表面張力をモデル化した項を加えることを提案している。

また、インタラクティブに流体を描画できるように、高速に計算できるカーネルを設計すると共に、キャッシュヒット率が向上するようデータ構造を工夫している。

レンダリングにはPoint SplattingとMarching Cubesの２通りの手法を用いており、2003年当時のPC(Pentium4+GeForce4)で2000個程度の粒子を計算した例では、Point Splattingで20fps、Marching Cubesで5fpsでの描画を実現している。

ハンニバル戦記

塩野七生の「ローマ人の物語」を読んでいる。

今日は、２巻の「ハンニバル戦記」を読み始めた。ハンニバル戦記は、そのタイトルの通りハンニバルのアルプス越えで有名なポエニ戦役について書かれている。

ポエニ戦役は紀元前264年〜同133年までローマとカルタゴの間で行われた戦争で、今日読んだのはそのうちポエニ戦役のきっかけとなったメッシーナでの戦い。

紀元前265年、シチリア半島の東端に位置するメッシーナは、南に位置するシラクサから攻め込まれていた。

メッシーナは、海峡を挟んで対岸にあるローマへ救援を求める。頼られたローマは迷ったあげく、クラウディウス率いる軍を派遣した。

ローマはシラクサ、カルタゴを撃破。シラクサとは同盟関係を結ぶことでシチリア島に影響力を伸ばすことに成功する。

しかし、これによってカルタゴの勢力圏と直接接することになり、その後のポエニ戦役へと繋がっていく…

次はいよいよポエニ戦役の始まり！

由比ヶ浜

夕方に鎌倉の由比ケ浜までふらっと行ってきた。

鎌倉には何度かいったことがあるけど、由比ヶ浜にいったのは初めて。

行ってみると、広いし、水もきれいだし、設備もあるし、駅から近い。すごく良い！

海というと子供のころの印象が強くて、車で何時間もかけて頑張っていかないといけないというイメージがあったけど、由比ケ浜は全然そんなことがない。

すごくいいところを見つけた♪

粒子法のシーンを２倍のサイズにしてみたが…

３次元の粒子法シミュレーションを、さらにシーンを２倍のサイズにして計算してみた。



あ…カメラの位置を調節していなかった…
再計算します。

追記：再計算しました。こちら


■関連する記事
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた

３次元の粒子法シミュレーション

粒子法のプログラム第１回（概要）や【粒子法】粒子を流体としてレンダリングでは粒子を２次元で計算していたので、今回はそれを３次元に拡張して計算してみた。平面だったものを、そのままz方向にも粒子を配置している。


【粒子表現】


【流体表現】



■流体シミュレーションに関するエントリ
粒子法のプログラム第１回（概要）
粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
【粒子法】粒子を流体としてレンダリング
３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
SPHによる巻き波のシミュレーション２
Haskell、OCamlでSPH法
このあとやりたいこと
カメラ位置を変えて流体をレンダリング
固液連成シミュレーション

■剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション

粒子法(SPH)のプログラム

粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。

粒子法(SPH)のプログラム一覧

動画

シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧

--