粒子ベース剛体シミュレーション（プレビュー）

SPH法による粒子法流体シミュレーションに続き、粒子ベースの剛体シミュレーションを書いています。

少しずつ進めているのですが、いまの時点できているものを少しだけ載せてみます。

【回転するキューブ】


【並進衝突するキューブ】


粒子ベース剛体シミュレーションについても、流体シミュレーションと同様に説明エントリを書きたいと思っています！
-

「粒子法」で検索すると１ページめにこのブログが！

いつのまにか、Googleで「粒子法」を検索するとこのブログが１ページめに出てくるようになってました！

Googleで「粒子法」を検索

いろんな人に見てもらっているからですね。ありがたいことだー
-

シーグラフアジアに行きたい！

12月に横浜でシーグラフアジアが開かれるので、19日の土曜日に行ってこようと思ってます。

シーグラフアジアは、コンピュータグラフィクスとインタラクティブ技術を中心とした国際カンファレンス＆エキシビジョンで、CGの国際学会シーグラフのアジア版です。

エキシビジョンは1,000円でみれるので、そっちメインで楽しんできます！

■このブログのCG系エントリ
【粒子法】粒子を流体としてレンダリング
粒子法のプログラム第１回（概要）

そのほか、エレクトロニックシアターというのもあって、世界中からエントリされた優れたCG作品をみることができる２時間のショーです。こちらは2,500円と少し高いのでどうしようかなぁ？？

あと、FF13のメイキングセッションもあって見てみたいんだけど、16日の平日なので見れない。残念。。
-

bitly経由のクリックがGoogle Analyticsにあがらないのはなぜ？

Twitterにbitlyで短縮したリンクを貼って、そのリンク経由のアクセスをGoogle Analyticsで見てみようとしました。

ところが、bitlyには６クリックと出ているのに、AnalyticsにはTwitterやbitlyからの流入が１件もあがりません。

なぜだろう？

まずbitly経由のアクセスの場合のリファラがtwitterになるのかbitlyになるのかを調べてみると、リファラはhttp://bit.lyでした。なので、もしAnalyticsにあがるとすれば、bit.lyからの流入のはずです。しかし、bit.lyからのアクセスは１件もあがっていません。

もしかすると、TwitterからのアクセスはTwitterをクロールするbotからのもので、それをAnalyticsが除外しているのだろうか？

あるいは、TwitterクライアントについているURL Expanderによって短縮前のリンクだけが取得されており、それを同様にAnalyticsが除外しているのだろうか？

だれか、ご存知の方、教えてください！
-

お砂糖エステが流行るよ！

先日髪をカットしたときに、新しく始めたんですよと「お砂糖エステ」を薦められました。

エステ用の特別な砂糖をつかってマッサージしてもらったんですが、ほんとに肌がスベスベに！

で、つい昨日渋谷を歩いていると試供品を渡されて、それが同じくお砂糖エステでした。もらった試供品を使うと、まえと同じように即効でスベスベ！

シュクレ（SUCRE）という北海道産の新商品らしくて、渋谷のランキンランキンで売ってます。

香りもいろいろあるみたいで、シトラスとかバニラとか、北海道ならではのラベンダーも。

ぜったい流行るよ、これ！
--

Love Paranoia：柴咲コウ

柴咲コウのニューアルバム「Love Paranoia」

ジャケットの柴咲コウが可愛すぎる！
完全にジャケ買いしてしまった。

初回盤には、来春行われるライブツアーの優先抽選権がついていたので、さっそく申し込みました。あたるかなー？

iTunesからアクセスするにはこちら
-

コンピュータ将棋のいま

今日、東大駒場祭に行って、コンピュータ将棋の展示を見てきました。

展示されていたのはGPS将棋というソフトウェアで、東京大学大学院総合文化研究科の教員・学生が開催しているゲームプログラミングセミナー（Game Programming Seminar = GPS）のメンバーが中心になって開発しているものです。

プロ棋士の棋符を学習させたり手の評価を改善することで、ここ数年で将棋プログラムはかなり強くなっており、プロ４段相当といわしめるレベルまできているそうです。

将棋の他はどうなんですかと聞いたところ、チェスではカスパロフがディープブルーに負けてから熱気が下がってしまったとのこと。一方で囲碁はいま熱いそうで、モンテカルロ法を使ったアプローチによってここ１〜２年で急激に強くなっているそうです。

また、どうぶつ将棋についての展示も併設してあり、そちらは全局面での勝敗がすでに解析し終わっているそうなのですが、それを逆に活かして「絶対人間に勝たない」プログラムを組んで、将棋の普及に協力していました。

AdsenseとAmazonアフィリエイトを置いてみました

実験的にAdsenseとAmazonアフィリエイトを置いてみました。

目的は、どのくらいクリックされるのかと、その傾向をみるためです。

(2009.11.20追記　Adsenseの審査がおりなかった…)
-

このあとやりたいこと

このエントリでは、ここまでやってきたことを一度振り返ってみて、それからこの後やっていきたいことを書き綴ってみます。

ここまでやってきたこと

ここまで、粒子法、中でもSPH法を使って、流体のシミュレーションとその可視化をやってきました。

まず２次元でダムブレイクのシミュレーションをして、それをマーチングキューブを使って陰関数曲面を抽出しPovrayを使ってレダリングしました。
→粒子法のプログラム第１回（概要）
→【粒子法】粒子を流体としてレンダリング

そのあと、３次元に拡張して、同様にダムブレイクのシミュレーションをし、可視化しました。
→３次元の粒子法シミュレーション

３次元のダムブレイクについては、いくつかパラメータを変えてみたりしてどのように動きが変わるかを調べてみました。
→粒子法のシーンを２倍のサイズにしてみた

また、論文を参照して、プログラム中の式の意味を深掘ってみたりもしています。
→p = k ( ρ - ρ0 )でのρ0の意味

ダムブレイクのあとは、もっと他のシーンもレンダリングしてみたいと思い、巻き波のシミュレーションをしてみました。巻き波によっておきるしぶきは、粒子法の得意とするところです。論文を参考にしたりして、ムービーを作りました。
→SPHによる巻き波のシミュレーション第１回

SPH法のコードを、関数型言語のHaskellやOCamlでも実装してみました。
→Haskell、OCamlでSPH法

このあとやっていきたいこと 

このあとやりたいと思っていることはいろいろあります。

いまやろうとしているのは、粒子ベースの剛体シミュレーションです。たとえばたくさんのトーラスが落下するようなシーンを、年内くらいに実装できるよう進めています。それができたら、固液連成シミュレーションでいろいろなシーンを作りたいと思っています。

その他、
・弾性体シミュレーション
・気液連成シミュレーション
・SPH法の高速化
・PointSplattingなどの可視化方法
・Lattice Boltzmann Methodによる流体シミュレーション
・炎や煙のシミュレーション
・美しいシーン
・人工生命
のようなことを、やりたいこととして思っています。

何よりも、手法をいろいろ持っているだけでなく、その手法を使って綺麗なシーンをいっぱい描いていきたい。たとば、https://graphics.stanford.edu/wikis/cs348b-09/のような感じです。
-

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)のプログラム一覧

動画

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

動画の一覧

-