2011年12月31日土曜日

ケーパビリティを使ったリソース保護機構「Capsicum」


このエントリでは、ケーパビリティを使ったリソース保護機構である「Capsicum」を紹介します。Capsicumは、まもなく正式リリースされるFreeBSD 9に搭載される予定です。

Capsicumとは?
タネンバウムの「モダンオペレーティングシステム」では、オペレーティングシステムのリソース保護機構について、以下の3つを紹介しています。
  • 保護ドメイン
  • アクセス制御リスト
  • ケーパビリティ
これらのうち、ケーパビリティとは、改竄不可能なトークンを使ったアクセス制御ですが、これまでごく一部の研究用OSでしか実装されてきませんでした。

Capsicumとは、そのようなケーパビリティを使ったリソース保護機構を一般的なUNIXに導入しようというプロジェクトです。

Capsicumによるメリット
オペレーティングシステムのリソース保護機構については、DAC(Discretionary Access Control)やMAC(Mandatory Access Control)が広く使われていますが、これらは、1つのアプリケーションが様々な形式のリソースを扱う場合のリソース保護に向いたようには設計されていません。ここでいう様々な形式のリソースとは、たとえば、ウェブブラウザの場合、画像、動画、Flash、JavaScriptなどといったデータです。信用を受けたリソース元ではない可能性もあります。

そのような様々な形式のリソースを扱う局面おいてデータを保護するために、ウェブブラウザのChromiumでは、プロセスを複数に分けることでコンパートメント化を実現しています。ところが、この方法をとる場合、プラットフォームによってプロセス周りの仕様が変わるため、プログラマがそれぞれの環境に合わせてそれなりの量のコードを書いて実装しなければならないという問題があります。

この問題への解決策として、一般的なUNIXにケーパビリティを使ったリソース保護機構を補完的に導入しようというものが、Capsicumです。たとえば、Chromiumへのテスト実装では、わずか100行のコード追加で、コンパートメント化を実現できたといいます。

Capsicumは、FreeBSD 9に導入される予定です。12月9日にFreeBSD-9.0 RC3が公開されたので、来年の早いうちには正式リリースされるでしょう。

--

2011年12月28日水曜日

コンピュータシステムのアクセス制御における、2つの文脈


システムのアクセス制御といったとき、大きく2つの文脈に分けられます。オペレーティングシステムの分野における伝統的なアクセス制御と、企業内の情報システムにおけるアクセス制御の2つです。

後者は前者からの派生でありもちろん両者には関係性があるのですが、前者は悪意を持った攻撃者への対策という観点が主である一方、後者は複雑な組織内での適切な権限割当てという観点が大きくなっています。

オペレーティングシステムにおけるアクセス制御
タネンバウムの「モダンオペレーティングシステム」では、オペレーティングシステムのファイル保護機構は以下の3つが紹介されています。
  • 保護ドメイン
  • アクセス制御リスト
  • ケーパビリティ
保護ドメインは任意アクセス制御(Discretionary Access Control; DAC)、アクセス制御リストは強制アクセス制御(Mandatory Access Control; MAC)ともいわれます。

DACとMAC
従来のオペレーティングシステムにおけるアクセス制御として、DACとMACがあります。

DACはもともと、ユーザをお互いから守るように設計されていました。ファイルのようなオブジェクトの所有者が、そのオブジェクトへのパーミッションを指定できるという仕組みです。パーミッションは、そのオブジェクトがアクセスされるときに、OSによってチェックされます。

ところが、DACには脆弱性があります。root権限のような特権がアクセス制御を迂回できてしまうため、システムをクラックされて特権を取得されると重要なファイルを書き換えられてしまう可能性があるのです。これは、リソースの所有者とそのリソースへのアクセス権限設定者が分離されていないことによる問題です。

MACは、そういったDACの弱点を補い、システムのセキュリティポリシーを強化するよう設計されました。システム管理者がポリシーを決め、そのポリシーがOSのカーネルのあちこちに埋め込まれたフックによってチェックされます。

ケーパビリティ
DACやMACは、リソースに対してセキュリティレベルを設定するという考え方に基づく仕組みです。一方、プロセスに対して最小限のアクセス権限を与えることでシステムの安全性を高めるという考え方があります。それがケーパビリティです。

たとえば、UNIXにおけるプロセスには、一般ユーザ権限で動くものと特権で動くものがあります。もし、特権で動いているプロセスに脆弱性があった場合、不正な操作によりプロセスが持っている権限を取得されてしまう可能性があるのです。

この問題を解決する方法がケーパビリティです。特権をさらに細分化したケーパビリティと呼ばれる単位で取り扱えるようにし、プロセスに最小限のケーパビリティのみを与えて必要な処理のみを行わせるというものです。

この仕組みにより、もしプロセスに脆弱性がありそれが悪用されたとしても、そのプロセスが必要とする最小限のケーパビリティしか奪われないため、被害の範囲を狭めることが可能になります。プロセスが自らジェイルに入るようなものです。

ケーパビリティはこれまで、一部のセキュアOSでのみ使われていましが、最近の動きとして、Capsicumというケーパビリティ機構がFreeBSD 9に搭載されます。あらゆるシステムがインターネットに繋がれるようになる中、予期せぬ脆弱性からシステムを守るため、ケーパビリティのような保護機構の必要性が高まっているといえます。

企業内の情報システムにおけるアクセス制御
1990年ごろまで、軍事システムでMACが利用される一方、企業や行政組織ではMACではなくDACが適しているとされていました。しかし、DACには上述の脆弱性が存在します。企業や行政組織でITシステムが広く使われるようになるにつれ、それらの組織でも情報セキュリティがより重要となりました。

そんな中、DACの脆弱性を補いつつも、企業や行政組織でも運用が可能なアクセス制御の仕組みとして、ロールベースアクセス制御(Role-based Access Control; RBAC)が考えだされました。

DACやMACでは、アクセス制御の主体は個々のユーザやグループ、対象はファイルのような具体的なリソースでした。一方、RBACでは、パーミッションはユーザやグループではなく「役割(ロール)」について付与され、ユーザやグループはそれらの「役割」に紐づけられます。また、対象も具体的なリソースではなくそれらのリソースに対して一連の操作を行う「オペレーション」となります。

このように、「役割」や「オペレーション」といった、一段抽象度を上げた概念に対してアクセス制御を行うことで、RBACでは、企業および行政組織のポリシーや構造に合わせた柔軟な、しかし集中的なコントロールを可能としています。企業や行政組織での実運用が重視されているのがポイントです。


--

Compiling the Particle-based Simulation Description Language


On the former entry, I described simple N-Body simulation in the Particle-based Simulation Description Language as a embedded DSL and executed it with interpreter. In this entry, I introduce compiling and demonstrate it.

Implementation
The Particle-based Simulation Description Language is compiled to Haskell using Template Haskell, a GHC extension to Haskell that adds compile-time metaprogramming facilities.

With this transformation, a code corresponding to the N-Body simulation code written in plain Haskell on the former entry is produced, then it is compiled with GHC Haskell compiler.

In other words, the host language is the Particle-based Simulation Description Language and the target language is Haskell, though there is no clear line between the host and target languages because it is a embedded DSL.

Codes
Here the codes are, with the interpreter introduced on the former entry together, arranged as a module.

The structure of the module is as below:
- Examples/
- SimulationDSL/
  - SimulationDSL.hs
  - SimulationDSL/
    - Compiler/
    - Data/
    - Interpreter/
    - Language/
    - Test/

Execution
To execute, compile CompilerRun.hs in Examples directory and run the output file as below:
$ cd Examples
$ ghc -i../SimulationDSL --make -O2 CompilerRun.hs
$ ./CompilerRun

Performance
Comparing performance on three codes:
  • with compiler - CompilerRun
  • in plain Haskell - NBody
  • with interpreter - InterpreterRun
The performance:
$ time ./CompilerRun > /dev/null

real 0m0.949s
user 0m0.927s
sys 0m0.011s

$ time ./NBody > /dev/null

real 0m0.972s
user 0m0.950s
sys 0m0.011s

$ time ./InterpreterRun > /dev/null
real 0m2.932s
user 0m2.906s
sys 0m0.022s

As chart:


As intended, the same performance to the code in plain Haskell is given with compiler, which is 2x-3x faster than with interpreter.

Conclusion
In this entry, I introduced compiling the Particle-based Simulation Description Language to Haskell and showed improvement in performance of simple N-Body simulation.

While this time I chose Haskell as the target language for this embedded DSL, it is also possible to compile it to C adapting OpenMP or CUDA with parallel computation.


--

2011年12月13日火曜日

粒子ベースシミュレーション記述言語のコンパイル


前回のエントリでは、粒子ベースシミュレーション記述言語を使ってN-Body(多体問題)シミュレーションを記述し、それをインタプリタとして実行しました。

今回は、それをコンパイルして実行するようにしたので紹介します。

実装方針
実装の方法としては、Template Haskellを使ってDSLをHaskellのコードに変換しています。この展開によって、前回のエントリでHaskellで直接書いた場合と同様のコードを生成し、さらにghcでコンパイルしています。

つまり、ホスト言語が粒子ベースシミュレーション記述言語、ターゲット言語がHaskellということになります。もっとも、組み込みDSLなので明確な線引きがあるわけではありませんが。

コード
コードはこちらに置いてあります。前回のインタプリタとあわせて、モジュールとして整理しました。

モジュールの構成は以下のようになっています。
- Examples/          実行例をまとめています
- SimulationDSL/
  - SimulationDSL.hs   モジュール外へのエクスポート
  - SimulationDSL/
    - Compiler/        コンパイラに関するコード
    - Data/            基本的なデータ構造
    - Interpreter/     インタプリタに関するコード
    - Language/        シミュレーション記述言語の構文と解釈
    - Test/            テストを追加していく予定

実行方法
以下のように、ExamplesフォルダのCompilerRun.hsをコンパイルし、ファイルを実行します。
$ cd Examples
$ ghc -i../SimulationDSL --make -O2 CompilerRun.hs
$ ./CompilerRun

実行効率の確認
実行速度を、以下について比べてみます
  • コンパイラを使った場合 - CompilerRun
  • Haskellで直接書いた場合 - NBody
  • インタプリタを使った場合 - InterpreterRun

実行速度は以下のようになりました。
$ time ./CompilerRun > /dev/null

real 0m0.949s
user 0m0.927s
sys 0m0.011s

$ time ./NBody > /dev/null

real 0m0.972s
user 0m0.950s
sys 0m0.011s

$ time ./InterpreterRun > /dev/null
real 0m2.932s
user 0m2.906s
sys 0m0.022s

グラフにすると以下のようになります。



狙い通り、コンパイラを使った場合は、インタプリタのx2〜x3の速度、つまり、Haskellで直接書いた場合と同等の速さで動いています。

まとめ
このエントリでは、粒子ベースシミュレーション記述言語をHaskellにコンパイルしました。その結果、実行速度が改善されたことを示しました。

今回は組み込みDSLとしてターゲット言語をHaskellにしましたが、Cにコンパイルすることや、さらにOpenMPやCUDAに対応させて並列計算を可能にすることも検討しています。

また、HackageDBにAccelerateという配列をCUDAで計算するライブラリがあるようなので、これを使う方法もありそうです。

逆の方向として、粒子ベースシミュレーション記述言語の表現力を強化することも考えられます。
  • 境界条件の指定
  • 効率的な近傍探索のサポート
  • 例外処理
  • デバッガ
  • プロファイラ
  • 階層構造の導入(剛体シミュレーション)
  • ドメインの導入(固液連成)

ラグランジュ的表現とオイラー的表現
この言語はラグランジュ的な表現を記述することができます。一方、オイラー的な表現を記述することはできません。

有限差分法や有限要素法、あるいは格子ガス法や格子ボルツマン法といった、格子およびメッシュをベースにした表現を使ったシミュレーションを記述するには、また別の言語モデルを設計することになるでしょう。


関連エントリ
このエントリに関連するページです。

粒子法シミュレーション記述言語によるn-Body(多体問題)シミュレーション
粒子法シミュレーションを記述するための言語を設計中
【まとめ】埋め込みDSL関係の論文(1)
【まとめ】埋め込みDSL関係の論文(2)


まとめページ

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

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




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

動画の一覧






--

2011年12月2日金曜日

N-Body simulation using the Particle Simulation Description Language


On the former entry (in Japanese), a free-fall simulation on point mass is shown as a starting point towards designing a language to describe particle simulation, or SPH.

On this entry, I introduce a improved Language and demonstrate N-Body Simulation on point mass using it.

Describing N-Body Simulation
N-Body simulation is represented in mathematical formula as below:



where xi, vi, ai and fi are position, velocity, acceleration and force on point mass i, respectively.

With the introduced language, as embedded DSL in Haskell, this simulation can be described as below:
nBody =
  do define         "x" (integral (var "v" <*> dt))
     define         "v" (integral (var "a" <*> dt))
     defineWithType "a" (var "f" </> m) ExpTypeVector
     define         "f"
                    (let r = norm (var' "x" <-> var "x")
                         k = m <*> m <*> g </> r </> r
                         n = (var' "x" <-> var "x") </> r
                     in sigma (k <*> n))
     initialConditionV "x" x0
  where dt = constantS 0.01
        m  = constantS 1
        g  = constantS 9.8
        x0 = V.fromList [ (0,0,0), (2,0,0), (0,2,0) ]

Explanation
define declaration defines an equation on physical amount on point mass i.

defineWithType declaration also defines an equation on physical amount, instructing type with it. When we use define declaration, the type of the equation we define is inferred with other equations. On the other hand, if the type can not be inferred, we must instruct it using defineWithType declaration. The types we can use are Scalar and Vector for now.

initialConditionV declaration is to give the initial condition of the physical amount defined with a "integral" expression. The last letter "V" means that the physical amount has Vector type. As the argument of it, an array of Vector is given as the initial condition of the position "x". Besides, the initial condition of the velocity "v" can be omitted, and an array of zero vector is given implicitly.

integral expression means integration, calculated explicitly inside the module.

var expression and var' expression are variables indicated by the symbol given as its argument. var' expression must be in a sigma expression. The difference in var and var' is described in the next paragraph.

sigma expression means summation of the calculated result of its operand expression on point mass j excluding itself. In the operand, var "x" means xi and var' "x" means xj.

dt, m and g are constant values lifted with the constantS function. The last letter "S" shows that the constant has Scalar type.

Execution
The nBody defined in the code can be evaluated as below:

main = mapM_ (printRow . flip machineRegisterValueV "x")
     $ take 100 $ runSimMachine nBody
* for the whole code, see Simulation.hs in GitHub

About performance efficiency
For now, this language is implemented as interpreter, so the performance of codes written in it is not very fast. I've compared it with the similar N-Body simulation code written in plain Haskell (see NBody.hs in GitHub). The result is below. Simulation is in the language and NBody in plain Haskell.

$ time ./Simulation > /dev/null

real    0m2.900s
user    0m2.869s
sys     0m0.020s

$ time ./NBody > /dev/null 

real    0m0.846s
user    0m0.829s
sys     0m0.011s

The program in the language is x1/3 - x1/4 as fast as in plain Haskell.

To improve the performance, using Template Haskell, we can compile, in the broad sense, the language to plain Haskell and make it work more efficiently. Also, compiling it to generate C or CUDA code would be pretty nice :)

--

2011年12月1日木曜日

粒子法シミュレーション記述言語によるn-Body(多体問題)シミュレーション


以前のエントリでは、粒子法シミュレーションを記述するための言語として、まず1つの質点についての自由落下を記述できるようにしてシミュレーションする例を示しました。

今回は、それを進めて、複数質点についてのn-Bodyシミュレーションを記述できるようにしたので紹介します。

n-Bodyシミュレーションの記述
n-Bodyシミュレーションを、数式で以下の様に表現した場合を考えます。



ここで、xi, vi, ai, fiは、それぞれ、ある質点iの位置、速度、加速度、力です。mは質点の質量、gは重力定数です。質点の質量はすべての質点について同じとします。

提案する言語を使うと、このシミュレーションを、以下のようにHaskellの組み込みDSLとして記述できます。
nBody =
  do define         "x" (integral (var "v" <*> dt))
     define         "v" (integral (var "a" <*> dt))
     defineWithType "a" (var "f" </> m) ExpTypeVector
     define         "f"
                    (let r = norm (var' "x" <-> var "x")
                         k = m <*> m <*> g </> r </> r
                         n = (var' "x" <-> var "x") </> r
                     in sigma (k <*> n))
     initialConditionV "x" x0
  where dt = constantS 0.01
        m  = constantS 1
        g  = constantS 9.8
        x0 = V.fromList [ (0,0,0), (2,0,0), (0,2,0) ]

記述の説明
"define"によって、ある質点iの物理量について方程式を定義しています。

"defineWithType"は、"define"と同じように物理量を定義しますが、併せて型を指定する点が異なります。"define"を使った場合には、その式の型は他の式から推論されますが、他の式から型を推論できない場合には、"defineWithType"を使って型も併せて指定します。現在使える型は、スカラーとベクトルの2種類です。

"initialConditionV"は、integralを使って定義された物理量の初期条件を与えるためのものです。最後の"V"は、その物理量の型がベクトルであることを示しています。引数として、ベクトルの配列を渡しています。速度vについては初期条件が省略されており、その場合は0ベクトルが使われます。

"integral"は、積分を表します。内部では陽解法で計算しています。

"var", "var'"は、引数で与えたシンボルで表される変数です。var'はsigmaの中でのみ使うことができます。varとvar'の違いについては、以下のsigmaの中で説明します。

"sigma"は、自身を除いた質点jについてオペランドの式を計算した値の合計です。オペランド中の var "x" がxiを、var' "x" がxjを表します。iが一定で、jをi以外の各質点について変えた場合の値を合計したものを返します。

dt、m、gといった定数は、constantS関数でリフトさせて使います。最後の"S"は、その定数がスカラーであることを示しています。

実行
上記で定義したnBodyを、以下のようにrunSimMachine関数に渡すことで計算します。
main = mapM_ (printRow . flip machineRegisterValueV "x")
     $ take 100 $ runSimMachine nBody
※コード全体についてはこちらを見てください

実行速度について
現時点では、この言語はインタプリタとして実装されているため、実行速度はあまり速くありません。

同じシミュレーションをHaskellで直接書いた場合(ソースコード)と速度を比較すると、以下のようになりました。SimLangがこの言語、NBodyがHaskellで直接書いた場合です。

$ time ./SimLang > /dev/null

real    0m2.900s
user    0m2.869s
sys     0m0.020s

$ time ./NBody > /dev/null 

real    0m0.846s
user    0m0.829s
sys     0m0.011s
Haskellで直接書いた場合に比べて、1/3〜1/4の速度になっています。

実行速度については、Temlate Haskellを使って、Haskellで直接書いた場合と同等のコードに(広義の)コンパイルすることで効率化できます。

課題
現時点での実装には、以下のような課題があります。
  • 0割算などの例外処理
  • デバッグのしやすさ、記述にエラーがあったときのメッセージの改善

ロードマップ
今後のロードマップとして以下のようなことを考えています。
  • より強力な型推論
  • Template Haskellによる(広義の)コンパイル
  • テスト
  • 設計の解説
  • 組み込みDSLや言語抽象についての解説
  • SPHの記述(分岐や近傍探索の導入)
  • 剛体シミュレーションの記述(階層表現の導入)
  • 固液連成シミュレーションの記述(ドメインの導入)
  • ベンチマーク

関連エントリ
このエントリに関連するページです。

粒子法シミュレーションを記述するための言語を設計中
【まとめ】埋め込みDSL関係の論文(1)
【まとめ】埋め込みDSL関係の論文(2)


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

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



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

動画の一覧





--

2011年11月29日火曜日

このブログのGitHubレポジトリを作りました


このブログのGitHubのレポジトリを作りました。

blog-codes@GitHub

これまでのブログエントリで公開していたコードを、GitHubのレポジトリに移しました。

粒子法のプログラム第1回(概要)
sph.cpp - GitHub

エントリはこちら。粒子法のプログラムを解説しています。論文や本の数式をプログラムに落とすところやパラメータの設定が分かるよう、ソースコードも合わせて公開しています。

Haskell、OCamlでSPH法
sph.hs - GitHub
sph.ml - GitHub
sph.cpp - GitHub

エントリはこちら。粒子法のコードを、Haskell、OCaml、C++で書いて比較しています。

【OCaml】スタンフォードバニーの読み込みライブラリ
PlyReader.ml - GitHub
PlyReader.mli - GitHub

エントリはこちら。OCamlで、スタンフォードバニーの読み込みライブラリを書きました。

粒子法シミュレーションを記述するための言語を設計中
Simulation.hs - GitHub

エントリはこちら。粒子法シミュレーションを記述するための言語を考えており、それを紹介しています。

CPUキャッシュによるプログラム高速化の例
code1.c - GitHub
code2.c - GitHub

エントリはこちら。CPUのキャッシュメモリを使うで、プログラムの実行性能をあげることができます。キャッシュメモリを活用することでどの程度性能が向上するのかを紹介しています。


--

2011年11月24日木曜日

Acceleration of programs using CPU caches effectively


Effective use of CPU caches make the performance of programs improved. This entry shows how much improvement is given with effective use of CPU caches.

Motivation
Computer data storages have hierarchy which consists of registers, cache, main memory and so on, distinguished by response time. Registers and cache memory are smaller but faster than main memory.

CPUs can directly operate registers only. So, between registers and main memory, there is about x100 gap in access speed, and to fill the gap, CPUs have cache memory inside. With effectively using CPU cache, we can improve the performance of programs.

However, programmers can not make CPUs directly to use the cache memory in the way they want, so they must write programs considering how the CPUs read/write the data on the memory to indicate CPUs indirectly how to use the cache.

* SPARC64 VIIIfx used in the K computer of RIKEN in Japan has "sector cache" functionality and programmers can directly indicate it how to use the cache. With sector cache, the cache in a CPU is divided in two sectors and programmers can instruct which sector they use to store data.

This entry demonstrates an example of performance improvement using CPU cache effectively.

Example
I show the optimization in efficiency of a program which computes the multiplication of 1000x1000 matrices using cache blocking technique.

The code before optimized (code1.c):
#include <stdio.h>
#include <stdlib.h>
 
#define N 1000
 
int main ()
{
  double A[N][N], B[N][N], C[N][N];
  int i, j, k;
  
  // initialize the matrices
  for ( i = 0; i < N; i++ )
    for ( j = 0; j < N; j++ )
      {
        A[i][j] = 1.0;
        B[i][j] = 1.0;
        C[i][j] = 0.0;
      }
  
  // matrix multiplication
  for ( i = 0; i < N; i++ )
    for ( j = 0; j < N; j++ )
      for (k = 0; k < N; k++ )
        C[i][j] += A[i][k] * B[k][j];
  
  printf ( "%f\n", C[0][0] );
  
  return 0;
}

I rewrite this code as below with cache blocking (code2.c). Of course, the result of multiplication is same as to that of code1.c.
$ cat code2.c
#include <stdio.h>
#include <stdlib.h>
 
#define N 1000
 
int main ()
{
  double A[N][N], B[N][N], C[N][N];
  int i, j, k;
  
  // initialize the matrices
  for ( i = 0; i < N; i++ )
    for ( j = 0; j < N; j++ )
      {
        A[i][j] = 1.0;
        B[i][j] = 1.0;
        C[i][j] = 0.0;
      }
  
  // matrix multiplication
  int ibl = 100;
  int ib, jb, kb;
  for (ib=0; ib<N; ib+=ibl)
    for (jb=0; jb<N; jb+=ibl)
      for (kb=0; kb<N; kb+=ibl)
        for (i=ib; i<ib+ibl; i++)
          for (j=jb; j<jb+ibl; j++)
            for (k=kb; k<kb+ibl; k++)
              C[i][j] += A[i][k] * B[k][j];
  
  printf ( "%f\n", C[0][0] );
  
  return 0;
}

Demonstration
The times taken to execute the unoptimized code and the optimized code are as below respectively:
$ gcc -Wl,-stack_size,1000000000 -O2 -o code1 code1.c
$ time ./code1 > /dev/null

real    0m7.595s
user    0m7.494s
sys     0m0.035s

$ gcc -Wl,-stack_size,1000000000 -O2 -o code2 code2.c
$ time ./code2 > /dev/null

real    0m2.509s
user    0m2.462s
sys     0m0.030s
* The -Wl option is to designate ld to allocate large arrays on the stack

The result represented as a chart:



Cache blocking improve the performance by about x3.

Comparing the performance as the size of cache blocking
Next, I show how the performance changes as the size of block. I changed the size of the block to 100, 200, 250, 500, then timed to calculation.

As a chart:


With the block size of 500, the performance degrades massively, which represents that the expansion of the block size leads to overflow of data from the cache.

The bytes of the 500x500 sized block is 2MB, with double typed array, so, considering other programs running in background, it agrees with the CPU cache size of my computer :)

--

2011年11月23日水曜日

CPUキャッシュによるプログラム高速化の例


CPUのキャッシュメモリを使うで、プログラムの実行性能をあげることができます。このエントリでは、キャッシュメモリを活用することでどの程度性能が向上するのか見てみます。

動機
コンピュータの記憶領域は、レジスタ、キャッシュメモリ、メインメモリなどからなる階層構造になっています。レジスタやキャッシュメモリはメインメモリより容量が小さいものの、非常に高速です。

CPUが直接触れるのはレジスタのみです。レジスタとメインメモリでは、データアクセスの速度に約100倍の差があり、その差を埋める目的で、CPU内部にキャッシュメモリがあります。キャッシュメモリを活用できれば、プログラムの実装性能を大きく高めることができます。

ただし、基本的に、CPUがどのようにキャッシュメモリを使うのかをプログラマが直接指示することはできません。なので、キャッシュメモリを活用するには、CPUがどのようにメモリ上のデータを使うのかを考慮したプログラムを書くことで間接的にキャッシュメモリを使うことになります。

※理化学研究所の京で使われているSPARC64 VIIIfxにはセクタキャッシュという機構があり、プログラマがキャッシュを制御できます。セクタキャッシュを使うと、キャッシュメモリを2つの領域に分割し、データごとにキャッシュに使う領域をプログラマが指定できます。

このエントリでは、キャッシュメモリを活用することで、実際にどの程度性能が向上するのか、一例を見てみます。


サンプルコード
1000×1000行列の積を求めるプログラムについて、キャッシュブロック化によって実行性能を最適化します。

以下が、最適化前のCのコードです。(code1.c)
#include <stdio.h>
#include <stdlib.h>
 
#define N 1000
 
int main ()
{
  double A[N][N], B[N][N], C[N][N];
  int i, j, k;
  
  // initialize the matrices
  for ( i = 0; i < N; i++ )
    for ( j = 0; j < N; j++ )
      {
        A[i][j] = 1.0;
        B[i][j] = 1.0;
        C[i][j] = 0.0;
      }
  
  // matrix multiplication
  for ( i = 0; i < N; i++ )
    for ( j = 0; j < N; j++ )
      for (k = 0; k < N; k++ )
        C[i][j] += A[i][k] * B[k][j];
  
  printf ( "%f\n", C[0][0] );
  
  return 0;
}

これを、キャッシュブロック化によって以下のコードに書き換えます(code2.c)。もちろん、計算の結果そのものは変わりません。
$ cat code2.c
#include <stdio.h>
#include <stdlib.h>
 
#define N 1000
 
int main ()
{
  double A[N][N], B[N][N], C[N][N];
  int i, j, k;
  
  // initialize the matrices
  for ( i = 0; i < N; i++ )
    for ( j = 0; j < N; j++ )
      {
        A[i][j] = 1.0;
        B[i][j] = 1.0;
        C[i][j] = 0.0;
      }
  
  // matrix multiplication
  int ibl = 100;
  int ib, jb, kb;
  for (ib=0; ib<N; ib+=ibl)
    for (jb=0; jb<N; jb+=ibl)
      for (kb=0; kb<N; kb+=ibl)
        for (i=ib; i<ib+ibl; i++)
          for (j=jb; j<jb+ibl; j++)
            for (k=kb; k<kb+ibl; k++)
              C[i][j] += A[i][k] * B[k][j];
  
  printf ( "%f\n", C[0][0] );
  
  return 0;
}


実行
最適化前、最適化後それぞれについて、実行にかかる時間は以下のようになりました。
$ gcc -Wl,-stack_size,1000000000 -O2 -o code1 code1.c
$ time ./code1 > /dev/null

real    0m7.595s
user    0m7.494s
sys     0m0.035s

$ gcc -Wl,-stack_size,1000000000 -O2 -o code2 code2.c
$ time ./code2 > /dev/null

real    0m2.509s
user    0m2.462s
sys     0m0.030s
※大きな配列をスタック上に確保するために、ldにオプションを渡してスタックサイズを大きくしています

グラフにすると以下のようになります。



キャッシュブロック化によって、約3倍の性能向上が得られました。


ブロックサイズによる性能比較
次に、ブロック化のサイズによって性能がどのように変わるかを調べてみます。ブロック化のサイズを100、200、250、500に変えて、実行時間を計ると、以下のようになりました。
100 : 2.311s
200 : 2.390s
250 : 2.509s
500 : 6.588s
※ブロック化のサイズの単位は、配列の要素数です。バイト数に直すと、double型の配列なので8を掛けたものになります

グラフにすると、以下のようになります。


ブロック化のサイズが500になると、急激に性能が悪化しています。このことから、ブロック化のサイズが大きくなりキャッシュにデータが収まらなくなった様子が分かります。

500×500のブロックのデータ量が2MBなので、このプログラム以外が動いていることも勘案すると、使っている環境のCPUの持つキャッシュメモリサイズと一致した結果といえます。


--

2011年11月22日火曜日

Common Lispでの、依存関係を考慮したファイルロード


Common Lispにおいて、asdfを使って、依存関係を考慮してプロジェクトの構成ファイルをロードする方法をまとめました。

動機
よっぽど小さなプログラムでない限り、プログラムの規模が大きくなってくると、ファイルを複数に分けることになります。

困ったことに、Common Lispには、依存関係をみて必要なファイルをロードする機能がありません。なので、ファイルを複数に分けると、必要なファイルを1つ1つロードするはめになります。

asdfを使うことで、ファイルの依存関係をみて必要なファイルを1度にロードすることができます。

asdfとは?
asdfは、"Another System Definition Facility"の略で、Common Lispプログラムを1つの「モジュール」として管理し、読み込むことができるツールです。現在はASDF2が広く使われています。

makeや、JavaのAntに相当するもので、PerlのCPANやRubyのRubyGemsの1つ下のレイヤに位置づけられます。CPANやRubyGemsにあたるものとしては、asdf-installがあります。

以下のような「sample」プロジェクトがあるとします。
sample
 - caller.lisp  ;; callee.lispで定義されている関数を使う
 - callee.lisp  ;; caller.lispから呼ばれる
 - sample.asd   ;; asdfモジュール定義ファイル

各ファイルの内容は以下のようになっています。
$ cat caller.lisp
(in-package :common-lisp)

(defpackage sample.caller
  (:use :common-lisp
        :sample.callee))
(in-package :sample.caller)

(hello-world)

$ cat callee.lisp 
(in-package :common-lisp)

(defpackage sample.callee
  (:use :common-lisp)
  (:export :hello-world))
(in-package :sample.callee)

(defun hello-world ()
  (print "Hello, World!"))

$ cat sample.asd 
(defsystem "sample"
  :components ((:file "callee")
               (:file "caller"))
  :serial t)

sample.asdがasdfから参照されるように設定し、common lispを起動します。

以下の関数を評価することで、sampleプロジェクトをロードすることができます。sample.asdの:serialオプションによって、ファイルの依存関係を考慮した形でファイルがロードされます。
(asdf:load-system 'sample)


--

2011年11月7日月曜日

TPPと混合診療の関係を簡単にまとめ


TPPと混合診療の関係を簡単に調べてみたので、そのまとめ。

重要なポイントとして、TPPによるアメリカの狙いの最も大きなものの1つが、サービス(金融)の自由化によって日本の個人資産を喰うこと。これは間違いない。わざわざサービス(金融)を後から捩じ込むだけのボリュームがあるターゲットは、日本だけだから。

アメリカ観点(混合診療を解禁しろ派)
混合診療解禁
→皆保険枠と混合診療枠ができる
→混合診療枠が医療保険のターゲット
→医療保険のノウハウを持っているのはアメリカの保険会社
→アメリカの保険会社はウマウマ

ref. http://haetarou.web.fc2.com/TPP/TPP_Kaisetu.html

※サービス(金融)はもともとTPPには入っておらず、アメリカが後から追加したもの。間違いなく、日本の個人資産を喰おうと狙っている

医師会観点(混合診療を解禁するな派)
混合診療解禁
→「政府は、財政難を理由に、保険の給付範囲を見直そうとしている」
→「混合診療を認めることによって、現在健康保険で見ている療養までも『保険外』とする」
→高額な診療が保険外になる
→高額な診療を受ける患者が減る
→医者の取り分が減る。困る

※医者は、国の医療費からたくさん診療費をもらいたい

ref. http://www.med.or.jp/nichikara/kongouqa/

厚生労働省観点(混合診療を解禁しろ派)
建前上、財政健全化したい
→医療費を減らしたい
→混合診療を解禁することによって、保険対象を減らすことで、医療費を減らす

※現状では、混合診療を認めていないため、その分手厚く保険対象にしてある

--

2011年10月25日火曜日

粒子法シミュレーションを記述するための言語を設計中


粒子法シミュレーションを記述するための言語を考えている。

最初のステップとして、質点の自由落下を表現できるようになったので、紹介する。

いま考えている言語を使うと、自由落下は以下のような記述で表現できる。Haskell上のEDSLとして実装されており、質点の座標の時間変化を出力する。

import Simulation

main = print $ take 100
             $ map (flip machineRegisterValue "x")
             $ runMachine (freeFall 0 0)

freeFall x0 v0 = do
  define "x" (constant x0 + integral (var "v" * constant dt))
  define "v" (constant v0 + integral (var "a" * constant dt))
  define "a" (constant g)

g  = -9.8
dt = 0.01

上でインポートしているSimulationモジュールはこちらからダウンロードできる。まだまだ作りかけのプロトタイプである。

言語の設計にあたって狙いとしているのは、「数式を記述すると、それがそのままシミュレーションとして動く」ことである。freeFallの定義において、数式に近い形で宣言的にシミュレーションを記述できていることに注目して欲しい。

位置x, 速度v, 加速度aの定義の順番は任意である。プログラマは、計算の順番ではなく、数式の意味に従って定義を記述できる。計算の順番は、モジュールが自動的に判断する。このような「表現と計算の分離」がこの言語の基本となるコンセプトである。

数式を記述するという点では、Mathematicaに代表されるようなシンボル計算と似ているが、そういったソフトウェアはナンバークランチングな数値計算のためというよりは、数式を正確に解くためのものである。その意味で、ここで考えている言語はMathematicaのようなものとは異なる。

現時点では、記述したシミュレーションはHaskellプログラムとしてインタプリタ的に実行されるが、将来的にコンパイラを作ることも可能である。コンパイラを作れば、たとえば、記述したシミュレーションをCUDAを使って実行するといったことができる。

いったん、上記では質点1つのシミュレーションとなっているが、直近の拡張として、物理量を集合として扱うことで、多粒子のシミュレーションを記述できるようにする予定である。多粒子のシミュレーションを記述できるようになれば、このブログでこれまで行ってきたようなシミュレーションを、この言語で改めて記述し直すことができる。

--

2011年10月19日水曜日

シミュレーション動画一覧


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


Falling Rigid Bunnies
たくさんのスタンフォードバニーを落下させる剛体シミュレーションを行い、結果をレンダリングしました。

動画を見る


Bunny-shaped fluid simulation
スタンフォードバニーの形をした流体を落下させて、結果をレンダリングしました。

動画を見る



粒子法による固液連成シミュレーション
粒子法で固体と液体の連成シミュレーションをおこない、結果をレンダリングしました。

動画を見る



固液連成シミュレーション
粒子法で固体と液体の連成シミュレーションをおこない、結果を粒子でレンダリングしました。

動画を見る



傾斜平面への剛体落下シミュレーション
傾斜平面への立方体の落下を、粒子ベースの剛体シミュレーションとして計算しました。

動画を見る



引き続き、粒子ベース剛体シミュレーション
粒子ベースの剛体シミュレーションとして、平面への立方体の落下を計算しました。

動画を見る



粒子ベース多体衝突シミュレーション
粒子ベースのアプローチで多体衝突の剛体シミュレーションをしました。

動画を見る



SPHによる巻き波のシミュレーション第3回
粒子法シミュレーション(SPH)によって、巻き波のシミュレーションをしてみました。砕波の様子が見られます。

動画を見る


粒子法のシーンを2倍のサイズにしてみた
3次元の粒子法シミュレーションを、スケールを2倍にして実施し、結果をレンダリングしました。

動画を見る


3次元の粒子法シミュレーション
粒子法シミュレーションを3次元で行い、結果をレンダリングしました。

動画を見る



【粒子法】粒子を流体としてレンダリング
2次元の粒子法シミュレーションの結果を、流体としてレンダリングしました。

動画を見る





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

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



論文まとめ



--

2011年10月8日土曜日

2011年9月24日土曜日

Falling Rigid Bunnies


Falling rigid bunnies, with particle-based simulation.

400 bunnies are falling, which consist of about 100 particles respectively.



粒子法のプログラム
粒子法のプログラム第1回(概要)
粒子法のプログラム第2回(プログラムの大枠)
粒子法のプログラム第3回(データ構造)
粒子法のプログラム第4回(密度と圧力の計算)
粒子法のプログラム第5回(力の計算)
粒子法のプログラム第6回(境界条件と粒子位置の更新)
粒子法のプログラム最終回(粒子の出力)
SPHによる巻き波のシミュレーション1
SPHによる巻き波のシミュレーション2
SPHによる巻き波のシミュレーション3
このあとやりたいこと

固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング
3次元の粒子法シミュレーション
粒子法のシーンを2倍のサイズにしてみたが…
粒子法のシーンを2倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション(プレビュー)
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

--

2011年8月4日木曜日

L-systemでのモデル生成



新しくPixarの中の人になったらしい人のエントリが面白かった。

Homage to Structure Synth

L-systemでモデルを生成して、RenderManでレンダリングしている。

--

2011年7月31日日曜日

構文チェックツールのFlymakeを、Haskell用に設定する方法


Emacsでは、Flymakeを使うことで、コードを書くと随時文法チェックが行われ、エラーのある行がバッファ上に表示されるようになる。

Haskell用にFlymakeを設定したので、その設定を書いておく。.emacsを編集するだけでOKだ。

Flymakeでは、基本的にコードと同じディレクトリにあるMakefileを使って文法チェックを行う形になっているが、いちいちMakefileを用意するのは煩わしいので、Makefileなしでも動作するようにした。特定のパラメータが必要な場合などは、Makefileを用意すれば自動的にそちらを使って文法チェックが行われるようにしてある。

.emacsの設定


以下のelispを.emacsに追加するだけでOK。

※Emacs21以前の場合は、flymake.elが標準で付属していないので、自分でインストールする必要がある。Emacs22以降であれば、標準で付属している。
;; Flymake Haskell

(require 'flymake)

(defun flymake-haskell-make-command (temp-file)
  (list "make"
        (flymake-haskell-make-parameters temp-file)))

(defun flymake-haskell-make-parameters (temp-file)
  (list "-s"
        "-C"
        "."
        (concat "CHK_SOURCES=" temp-file)
        "SYNTAX_CHECK_MODE=1"
        "check-syntax"))

(defun flymake-haskell-default-ghc-command (local-file)
  (list "ghc"
        (flymake-haskell-default-ghc-parameters
          (file-name-nondirectory local-file))))

(defun flymake-haskell-default-ghc-parameters (local-file)
  (list "-fno-code" local-file))

(defun makefile-exists-p (path)
  (file-exists-p (concat path "Makefile")))

(defun flymake-haskell-init()
  (let* ((temp-file  (flymake-init-create-temp-buffer-copy
                       'flymake-create-temp-inplace))
         (local-file (file-relative-name
                       temp-file
                       (file-name-directory buffer-file-name))))
    (if (makefile-exists-p (file-name-directory buffer-file-name))
        (flymake-haskell-make-command temp-file)
        (flymake-haskell-default-ghc-command local-file))))

(push '(".+\\hs$" flymake-haskell-init) flymake-allowed-file-name-masks)
(push '(".+\\lhs$" flymake-haskell-init) flymake-allowed-file-name-masks)
(push '("^\\(\.+\.hs\\|\.lhs\\):\\([0-9]+\\):\\([0-9]+\\):\\(\.+\\)" 1 2 3 4)
      flymake-err-line-patterns)

(add-hook 'haskell-mode-hook
    '(lambda ()
       (if (not (null buffer-file-name))
           (flymake-mode))))

説明


flymakeの有効化
add-hookによって、haskell-modeに入ったときに(flymake-mode)が実行される。

拡張子の対応付け、エラーメッセージの正規表現の定義
add-hookの上でpushしているのは、Haskellの.hsファイル、.lhsファイルとflymakeの初期化関数との紐付け、及び、ghcの出力から行番号やエラーメッセージを抽出するための正規表現である。

flymake haskellの初期化
flymake-modeに入ると、flymake-haskell-initが実行される。flymake-haskell-initでは、構文チェック用にコードをコピーし、そのコードを構文チェックするためのコマンドをflymake本体に返している。このとき、チェック対象のコードと同じディレクトリにMakefileがあるかどうかでコマンドを変えている。

チェック対象のコードと同じディレクトリにMakefileがある場合
この場合、構文チェックのコマンドはflymake-haskell-make-command関数によって生成される。flymakeがデフォルトで使用するmakeコマンドと同じものを生成している。

チェック対象のコードと同じディレクトリにMakefileがない場合
一方、この場合は、Makefileなしで構文チェックを行うよう、GHCを使ったコマンドを、flymake-haskell-default-ghc-command関数で生成している。生成するコマンドの形式は、
("ghc" ("-fno-code" 構文チェック用ファイルの名前))
となる。外側のリストのcarがghcコマンド、cadrがghcコマンドに渡すパラメータのリストである。これがflymakeに渡され、実行される。

補足


flymakeではエラーメッセージも表示してくれるのだが、ghcが複数行にわたるエラーメッセージを返す場合にはうまく処理できない。

flymakeは、それだけで完全にコンパイルをとおすためではなく、コードを書くときのタイポ回避程度に位置づけると良いだろう。

参考


http://www.emacswiki.org/emacs/FlymakeHaskell




--

2011年7月23日土曜日

【まとめ】埋め込みDSL関係の論文(2)


引き続き、埋め込みDSL(Domain Specific Language)関係の論文をまとめたので書いておく。

Template Meta-programming for Haskell (2002) - Tim Sheard, and Simon Peyton Jones
→Haskellでコンパイル時メタプログラミングを行うTemplate Haskellの論文

A Special-Purpose Language for Picture Drawing (1997) - Samuel N. Kamin, and David Hyatt
→MLを使って、2次元の絵を表現するための組み込みDSLを作る

Compiling Embedded Languages (2000) - Conal Elliott, Sigbjorn Finne, and Oege de Moor
→Kaminのプログラム生成のアイデアを元に、組み込みDSLの最適化コンパイラを示す。計算効率が重要な分野の応用例として、この方法を画像合成・操作のための言語Panに適用する。

Functional Implementations of Continuous Modeled Animation (1997) - Conal Elliott, and Paul Hudak
→連続的アニメーションのモデルを設計し、それを試作として実装してみせる。

Programming Graphics Processors Functionally (2004) - Conal Elliott
→3Dグラフィクスのための組み込みDSLであるVertigoを紹介する。VertigoはHaskellで書かれている。Vertigoは、GPUのコードを生成する最適化コンパイラである。

Functional Reactive Animation (1997) - Conal Elliott, and Paul Hudak
→インタラクティブでリッチなアニメーションを構築するためのフレームワークFRAN(Functional Reactive Animation)を紹介する。

Building Domain-Specific Embedded Languages (1996) - Paul Hudak
→良いソフトウェアを書くには抽象化が重要であり、最も"理想的"な抽象化は問題に合わせて設計されたプログラム言語を作ることであるとして、その障害である言語設計と実装の困難さを、組み込みDSLというアプローチで回避することを提案する。

Applicative programming with effects (2008) - Conor McBride, and Ross Paterson
→副作用を持つプログラムの抽象化の方法として、Monadよりも弱いため広く使えるApplicativeファンクタを紹介する。

Arrows, Robots, and Functional Reactive Programming - Paul Hudak, Antony Courteny, Henrik Nilsson, and John Peterson
→FRP(Functional Reactive Programming)の実例として、Haskellによる組み込みDSLであるYampaを紹介する。Yampaではモナドの一般化としてアローを使用する。Yampaによって、時間変化を扱うプログラムを宣言的なスタイルで実現できる。


■関連エントリ

【まとめ】埋め込みDSL関係の論文(1)

--

2011年7月8日金曜日

【まとめ】埋め込みDSL関係の論文(1)


埋め込みDSL(Domain Specific Language)関係の論文をまとめたので書いておく。

Plugging Haskell in (2004) - Andre Pang, Don Stewart, Sean Seefried, and Manuel M. T. Chakravarty
→Haskellをプラグインとして使う方法を紹介している。紹介されているプラグインは、HaskellでもHaskell FFIがサポートしている他の言語でも使用できる。GHCで提供されている仕組みを使って、プラグイン、実行時コンパイル、動的リンクを実現する方法を示している。また、サンプルも紹介されている。

Standard ML as a Meta-Programming Language (1996) - Samuel Kamin
→Standard MLをメタプログラミング言語として使うことを提案し、4つの実例を示している。Standard MLのコードからCのコードを生成。

Building program generators the easy way - Samuel Kamin
→関数型プログラミング言語で簡単にプログラムジェネレータを作ることができることを紹介している。この方法によって、ソフトウェア開発に関数型プログラム言語を便利に応用できる可能性を示している。

Effecient Image Manipulation via Run-time Compilation (1999) - Oege de Moor, Conal Elliott, and Sigbjorn Finne
→DSLとして定義したイメージ操作システムを、実行時に最適化・コンパイルすることで、その効率化する方法を紹介している。

Optimizing Embedded DSLs using Template Haskell (2004) - Sean Seefried, Manuel Chakravarty, and Gabriele Keller
→コンパイル時メタプログラミングによってEDSLの実行効率を改善する方法を示している。EDSLとしてのシンプルさを保ったまま実行効率を改善する。

GPU Kernels as Data-Parallel Array Computation in Haskell (2009) - Sean Lee, Manuel M. T. Chakravarty, and Vinod Grover
→HaskellにGPUカーネルを埋め込んでデータ並列の配列計算を行う方法を示している。

Accelerating Haskell Array Codes with Multicore GPUs (2011) - Manuel M. T. Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell, and Vinod Grover
→GPUコーンピューティングの抽象度を上げるために、配列計算の高レベルDSLを提案する。HaskellにこのDSLを埋め込み、NVIDIAのCUDAのためのコードを実行時に生成する。


■関連エントリ

【まとめ】埋め込みDSL関係の論文(2)

--

2011年7月6日水曜日

Language-oriented programming


A kind of programming that's quite common in Common Lisp and relatively uncommon in non-Lisp Languages.

Rather than provide an API built primarily out of functions, classes and macros, language-oriented programming provides language that you can embed in your Common Lisp programs.



--

2011年7月5日火曜日

LispマクロとTemplate Haskellの違い


LispマクロとTemplate Haskellには違いがある。

LispのマクロはS式を変換する。一方、Template Haskellが変換するものはHaskellのAST(Abstract Syntax Tree)である。

従って、Template Haskellで変換できる式は、まずASTとして有効でなければならない。ASTとして有効な式でなければTemplate Haskellに渡せない。Lispの場合はS式でさえあれば良い。シンボルをただ並べたリストでも良い。

制限が小さい分、Lispマクロの方が強力だといえるだろう。



--

2011年7月4日月曜日

【書評】実用Common Lisp

実用Common Lisp では、人工知能(AI)、コンピュータプログラミング技術、Common Lispを扱っている。原題はParadigms of Artificial Intelligence Programming : Case Studies in Common Lisp。

「第7章 STUDENT:代数の文章題を解く」
代数の文章題を解く。方法は、文章題を代数の数式に変換してそれを解くというもの。自然言語(英語)と数式の翻訳表をプログラムとして持つ。

「第8章 記号計算:簡単化プログラム」
Common Lispの数式処理システムであるMaximaのサブセットを作る。

「第9章 効率性の問題」
高速化の手法の1つとしてコンパイルをとりあげる。宣言的に記述された問題を、マクロを使ってコンパイル時にLispコードに変換する。

こういうのを「知能」だと見なしていた時代があったんだなぁ。

Lispは人工知能の分野で使われるというが、その理由が腹落ちした。この本によると、解き方が不明な問題に対するアプローチとしては、問題をそのまま宣言的に記述しそれをプログラムコードに変換するプログラムを書くのが適切であるという。そのような問題を扱う代表例がAIだった。そして、宣言的に記述された問題をプログラムコードに変換するのにLisp(のマクロ)が向いていたということだった。

2011年5月30日月曜日

【ドラフト】物理シミュレーションを対象としたDSELについての試み

Abstract


物理シミュレーションを対象としたDSELについての試みを紹介する。このDSELを使うことで、数式に極めて近い形で宣言的に定義を記述するだけで物理シミュレーションのプログラムを書くことができる。

組み込みドメイン固有言語(DSEL)とは


ドメイン固有言語(Domain-Specific Language, DSL)は、特定の問題領域に特化したプログラム言語を提供することでソフトウェア開発の生産性を向上させるものである。DSLは問題領域に適した抽象レベルでの記述を可能にするため、コード量の低減やそれに伴うメンテナンス性の向上といったメリットがある一方、新しく言語を作ること自体が高コストであるという問題がある。

その問題に対する解として、組み込みドメイン固有言語(Domain-Specific Embedded Language, DSEL)は、既存の汎用言語に組み込む形でドメイン固有言語を実装するものである。このアプローチによって、ホスト言語のライブラリやコンパイラ、デバッガ、プロファイラを再利用できるというメリットが得られる。

DSELについては、以下の文献がある。
・Building domain-specific embedded languages
・Huddakの、Haskore music notation
・conal elliot, Functional Image
・金融契約
・Haskell School of Expression

ここでは、物理シミュレーションを対象としたDSELについての試みを紹介する。ホスト言語にはDSELを作るのに向いているという理由でHaskellを選んでいる。

物理シミュレーションを対象としたDSEL


まず具体的な例を挙げよう。簡単な例として、質点の自由落下のシミュレーションを考えてみる。

自由落下を数式で表現すると以下のようになる。

... 数式 ...

提案するDSELを使うことで、質点の位置の時間過程を出力するプログラムを以下のように記述できる。(これは、そのまま実行可能なコードである)

import Process                -- DSELを定義したモジュール

main = print $ take 10 $ process x
                              -- 10ステップ分の位置を表示

x = x0 +: integral' (v * dt)  -- 質点の位置
v = v0 +: integral (a * dt)   -- 質点の速度
a = g                         -- 質点の加速度

x0 = 10.0                     -- 質点の初期位置
v0 = 0                        -- 質点の初期速度
dt = constant 0.01            -- 時間刻み
g  = constant (-9.8)          -- 重力加速度

これだけである。このコードをコンパイルして実行すると、10タイムステップ分の質点の位置が出力される。質点の位置、速度、加速度について、数式そのままの定義を宣言的に記述するだけであることに注目されたい。計算順序について考慮する必要はない。つまり、定義さえ記述すれば「実装」はもはや不要なのである。

もしこれをDSELを使わずにそのまま書くと、「aからvを計算し、vからxを計算する。それをタイムステップごとに繰り返す」という実装を自分で書くことになるだろう。ごく普通の、シミュレーションのプログラムの流れである。

DSELの設計


...DSELの設計について

DSELの実装


...DSELの実装について

DSELの適用可能条件


DSELの課題


・計算効率
 ・持ち上げられているため、コンパイラの最適化が効かない
  →メタプログラミングによるコンパイル時最適化 or コードジェネレータ
・スペースリーク
 ・condAを使うとスペースリークする

DSEL、現状での利用可能範囲


・配列まるごとをProcessの対象とする。そのなかの要素1つ1つは特定のタイムステップについて明示的に記述
 →上記の2つの問題は生じない
 →ただし、フルにDSELで書くことができず、DSELとpure Haskellのハイブリッドでの記述となる
  ・配列まるごと→DSELの世界での記述
  ・配列内の要素1つ1つについて→あるタイムステップについて明示的にpure Haskellの世界に下ろして記述

2011年4月26日火曜日

【まとめ】Haskellでの正格評価とWHNF


Haskellは遅延評価を特徴とする言語ですが、seq関数や!パターンを使うことで正格評価をさせることもできます。

ただし、最終的な値まで完全に評価(deepseq)されるわけではなく、WHNF(Week Head Normal Form)に簡約されるまでの評価です。

WHNFとは、簡単にいうと、以下のいずれかの形をもつ式です。
  • プリミティブである
  • 式の一番左にデータコンストラクタがある
  • 式の一番左にλがある
正格評価を使ったときにコードの挙動がどう変わるのか、具体的には以下のようになります。

ケース1(遅延評価)

$ cat case1.hs
loop 0 i = i
loop n i = loop (n-1) (i+1)

n = 100000000

main = let x = loop n 0
       in putStrLn "1"
$ ghc --make -O case1.hs
$ time ./case1
1

real    0m0.008s
user    0m0.001s
sys     0m0.003s
遅延評価の場合、let式の中のxは評価されません。評価されるのは文字の出力のみのため、実行時間はごくわずかです。

ケース2(xを正格評価)

$ cat case2.hs 
loop 0 i = i
loop n i = loop (n-1) (i+1)

n = 100000000

main = let !x = loop n 0
       in putStrLn "1"
$ ghc --make -O -XBangPatterns case2.hs
$ time ./case2
1

real    0m3.470s
user    0m1.873s
sys     0m0.032s
xに正格フラグ(!)を付けて正格評価にします。xは実際には使われませんが、評価はされるため、その分実行時間が長くなっています。

ケース3(代数データ型を使い、遅延評価)

$ cat case3.hs
data X = X Int

loop 0 i = i
loop n i = loop (n-1) (i+1)

n = 100000000

main = let x = X $ loop n 0
       in putStrLn "1"
$ ghc --make -O case3.hs
$ time ./case3
1

real    0m0.006s
user    0m0.001s
sys     0m0.003s
ケース1と同様に遅延評価ですが、loop関数の結果をデータコンストラクタXに適用してからxに束縛しています。ケース1と同様、xは評価されません。評価されるのは文字の出力のみのため、実行時間はごくわずかです。

ケース4(代数データ型を使い、正格評価)

$ cat case4.hs
data X = X Int

loop 0 i = i
loop n i = loop (n-1) (i+1)

n = 100000000

main = let !x = X $ loop n 0
       in putStrLn "1"
$ ghc --make -O -XBangPatterns case4.hs
$ time ./case4
1

real    0m0.006s
user    0m0.001s
sys     0m0.002s
ケース2と同様に、xに正格フラグ(!)を付けて正格評価にします。注目するのは、ケース2と異なり実行時間がほとんどかかっていない点です。

評価されるのはWHNFまでであるため、xにはデータコンストラクタについての関数適用がサンクとして束縛され、それ以上の評価はされません。そのため、loop関数は評価されていないのです。

もう1つ、ケースを追加してみましょう。

ケース5(代数データ型のフィールドを正格評価にする)

$ cat case5.hs 
data X = X !Int

loop 0 i = i
loop n i = loop (n-1) (i+1)

n = 100000000

main = let !x = X $ loop n 0
       in putStrLn "1"
$ ghc --make -O -XBangPatterns case5.hs
$ time ./case5
1

real    0m1.788s
user    0m0.997s
sys     0m0.020s
ケース4に加えて、型Xのフィールドに正格フラグを付けます。今度は、型Xのフィールドも評価されるため、実行時間が長くなっています。


以上のように、Haskellで正格評価を使う場合には、評価されるのがWHNFまでである点に注意が必要です。

ちなみに、WHNFまででなく完全に評価したい場合には、Control.Strategies.DeepSeqライブラリというものがあるので、そちらを検討してみると良いでしょう。

--

2011年3月3日木曜日

iPad2について検討


iPad2について。

重さをみてみると、600gある。iPad1とあんまり変わってない。

Kindleが、小さい方で240g、DXでも540gだから、iPadは重たいな。

2万円程度で、防水で、重くなくて、7インチくらいのウェブブラウジング端末なら買うのに。残念。

--