R言語をインストールしたので、練習に離散力学系の挙動をアニメーションにしてみました。以下の2つのMP4ファイルは、同じ現象を2つの方法で表現したものです。
- 3次元グラフによる表現 ddds-example-3d.mp4
メッシュを細かくしたヤツ ddds-example-3d_2.mp4 - ヒートマップによる表現 ddds-example-heatmap.mp4
内容：
- 行列計算と離散-離散・力学系
- グラフと遷移行列の生成
- 質量分布の描画
- 複数枚の絵をアニメーションにする
- 足し算・掛け算を変えてみたら
行列計算と離散-離散・力学系
計算としてやっていることは行列とベクトルの掛け算だけです。行列Tとベクトルvを与えて、v, Tv, TTv, TTTv, ..., Tnv という(n + 1)個のベクトルからなるシーケンスを作成します。それぞれのベクトルを絵にすると、(n + 1)枚の静止画が出来るので、それらをパラパラ漫画方式でアニメーションにします。
行列・ベクトルの掛け算の解釈（セマンティクス）として、離散力学的なモデルを想定します。そのほうが面白いし、自分が何を計算・描画しているかのイメージが持てますからね。
通常、離散力学系というと、時間は離散化しますが空間は連続のままです。ここでは、空間も離散化したモデルなので、離散-離散・力学系（discrete-discrete dynamical system）と呼ぶべきかもしれません*1。離散時間 0, 1, 2, ..., n に沿って離散空間内をナニカが動く現象を考えます。その動きをアニメーションで描き出すわけです。
ここで離散空間とは、有限有向グラフGだとします。GはN個の頂点を持つとします。G上を動くナニカを便宜上「物体」と呼ぶことにして、物体の状態は頂点集合上の質量分布で表されるとします。i番目の頂点に乗っている質量をv[i]とします。すると、v[1], v[2], ..., v[N] という長さNのベクトルによりある時点の物体が記述されます。
頂点j（j番目の頂点）に乗っている質量v[j]が次の時点で頂点iに移動する比率をT[i, j]とします。掛け算した量 T[i, j]v[j] は、「頂点j → 頂点i」と移動する質量です。次の時点に頂点iに乗る質量は、すべてのjからの寄与を足し合わせればいいので、ΣjT[i, j]v[j] です。この形を見ればわかるように、力学系の1ステップ時間発展は、行列Tの掛け算により与えられます。
グラフと遷移行列の生成
以下、記述にRの記法を使います。Rの書き方やプログラミングには独特の流儀（まだ僕も慣れてない）があるので違和感があるかも知れませんが、あまり気にしないでください。
まずは運動の舞台となるグラフを作ります。
> x <- runif(10)
> y <- runif(10)
これは、xとyにそれぞれ10個の一様乱数からなるベクトルを代入していて、結果としてランダムに配置された10個の点が得られます。x-y平面に10個の点をプロット plot(x, y) すると、例えば次のようになります。
辺も乱数で作ると、どうも具合が良くないので、決定性のアルゴリズムで作っています。出来るグラフは単純有向グラフ（多重辺を持たない）で、有向辺が下から上に向くようにしています。ただし、この前提はご都合主義的なもので何ら本質ではありません。
物体の移動を支配する行列Tを遷移行列と呼びます。遷移行列もグラフの形状から計算で求めています。物体の質量が保存するように、行列Tの各列の総和（縦方向に足した値） ΣiT[i,j] は1に正規化します。以下は、グラフの隣接行列（グラフのトポロジカルな形状を表す行列）とそれから作った遷移行列の例です。
> d10$adj
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 0 0 0 0 1 0 0 0
[2,] 0 1 0 0 0 0 0 1 1 0
[3,] 1 0 1 0 0 1 0 0 0 1
[4,] 0 0 0 1 0 0 0 0 0 0
[5,] 0 0 0 1 1 0 0 0 0 0
[6,] 1 0 0 0 0 1 1 0 0 0
[7,] 0 1 0 0 0 0 1 0 0 0
[8,] 0 0 0 0 0 0 0 1 1 0
[9,] 0 0 0 0 0 0 0 0 1 0
[10,] 0 0 0 0 0 0 0 0 0 1
> d10$trans
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.3333333 0.3333333 0 0.0 0 0.0 0.3333333 0.0 0.0000000 0.0
[2,] 0.0000000 0.3333333 0 0.0 0 0.0 0.0000000 0.5 0.3333333 0.0
[3,] 0.3333333 0.0000000 1 0.0 0 0.5 0.0000000 0.0 0.0000000 0.5
[4,] 0.0000000 0.0000000 0 0.5 0 0.0 0.0000000 0.0 0.0000000 0.0
[5,] 0.0000000 0.0000000 0 0.5 1 0.0 0.0000000 0.0 0.0000000 0.0
[6,] 0.3333333 0.0000000 0 0.0 0 0.5 0.3333333 0.0 0.0000000 0.0
[7,] 0.0000000 0.3333333 0 0.0 0 0.0 0.3333333 0.0 0.0000000 0.0
[8,] 0.0000000 0.0000000 0 0.0 0 0.0 0.0000000 0.5 0.3333333 0.0
[9,] 0.0000000 0.0000000 0 0.0 0 0.0 0.0000000 0.0 0.3333333 0.0
[10,] 0.0000000 0.0000000 0 0.0 0 0.0 0.0000000 0.0 0.0000000 0.5
>
質量分布の描画
質量分布v[i]（i = 1, 2, ..., N）を絵にしましょう。Rには、2変数の関数 z = f(x, y) を3次元にプロットして描いてくれるツールがあります。persp()関数です。見栄えがいいのでこれを使ってみます。
3次元プロットの仕組みを簡単に紹介します； 話を簡単にするために、x, yの範囲は、0≦ x, y ≦ 1 に固定します。区間[0, 1]をK等分して、0/K, 1/K, 2/K, ..., K/K = 1 と、両端を入れて(K + 1)個の点をサンプリングの基準にします。つまり、関数fの値は、f(0/K, 0/K), f(1/K, 0/K), ..., f(K/K, 0/K), f(0/K, 1/K), f(1/K, 1/K), ..., f(K/K, K/K) の合計(K + 1)*(K + 1)個の点で採取します。
α、βが 0, 1, 2, ..., K を動くインデックスだとすると、
- z[α, β] := f(α/K, β/K)
Rのインデックスは1始まりなので、s = 1, 2, ..., K + 1、 t = 1, 2, ..., K + 1 と動くインデックスs, tを導入して定義しなおしておきます。
- z[s, t] := f((s - 1)/K, (t - 1)/K) （1 ≦ s, tは整数 ≦ K + 1）
Rの関数persp()は、この行列zと、最初に決めたサンプリング点の情報を元に z = f(x, y) の関数グラフを描きます。x-y平面に方眼紙のような網目が描かれて、網目の交差点（格子点）におけるfの値が“z方向の高さ”に反映されます。次の図は、K=10としてx-y平面の網目だけ描いたものです。
さて、グラフGの頂点はx-y平面にばらまかれた質点と解釈できます。関数f(x, y)を、点(x, y)における質量だとして図示しましょう。
… んっ、あれ？ 困ったことになります。
グラフGの頂点がほとんど網目に引っかからないのですよ。頂点は大きさを持たない質点なので、よほどラッキーでないと網目の交差点に一致しません。サンプリングされなければ絵にもなりません。そこで、質点に大きさを持たせることにします。
質点をぼかす（あるいは、にじませる）処理に、次の関数を書きました。
blur <- function(vari, a, b, x, y) exp(-((x - a)^2 + (y - b)^2) / vari) blurredPoint <- function(vari, a, b) function(x, y) blur(vari, a, b, x, y)
質点（＝グラフGの頂点）の座標を (a, b) として、blur(vari, a, b, x, y) の (x, y) を動かすと、「ぼかした質量分布」となります。vari（分散；variance）は「ぼかす程度」を表すパラメータで、値を大きくすると「にじみ」が広がります。blurredPoint(vari, a, b) は、ぼかしパラメータと質点座標を固定した質量分布関数を返します。プログラミング的にはカリー化（引数の部分束縛）したblur()関数です。
blur()関数は、次の、簡略化したガウス分布を元にしています*2。
Vがぼかしパラメータ（通常、2σ2と書かれる量）で、rに、点(a, b)と点(x, y)の距離を代入すると、blur(V, x, y)になります。ぼかしは絵を描くためだけに行う処理で、離散力学系の計算には関係しません。したがって、ぼかした後の積分値を気にする必要はありません。ぼかし（にじみ）は質点の回りに出現するオーラで質量を持たないのです。
次の図は、2個の質点をぼかして描いたものです。
質量分布の可視化の方法としてヒートマップを使う場合も原理は同じです。質点をぼかしてからヒートマップを描きます。
実際に使ったグラフは20個の頂点と44本の辺*3を持つものです。一番下側に位置する頂点に質量1を乗せた初期状態から状態遷移させて、24枚の静止画を作って束ねています。
複数枚の絵をアニメーションにする
一枚の絵が描ければ、二枚も三枚も同じことです。初期状態となる質量分布ベクトルv[i]（i = 1, 2, ..., N、Nは質点の個数＝グラフの頂点数）を決めて、まずそれを描画します。遷移行列を掛け算して次の状態ベクトルv'[i]を求めて描画します。あとは繰り返し。
Rのグラフィックス機能そのものには、アニメーションを作る機能がありません。animationパッケージがアニメーション作成を支援してくれます。続けて描画した静止画を記録してくれて、それらの静止画像をまとめてアニメーションにしてくれます。
animationパッケージで生成できるデフォルトのアニメーションは、何枚かのスナップショットPNG画像とHTML／JavaScriptだけで実現するアニメーションです。これはブラウザだけで再生できる点はよいのですが、画像のプリロードでメモリを食うし、Firefoxではうまく再生できませんでした（Chromeでは大丈夫だった）。FFmpegをインストールすれば、それを使ってMP4やAVIファイルが作れます。結局、FFmpegを使う方法でMP4動画ファイルを作りました。
静止画や動画を作る際に調整できるファクターがたくさんあります。グラフの頂点の個数、辺の本数。頂点や辺、遷移行列を生成するアルゴリズム、ぼかしのパラメータ、色やレイアウトに関する種々のパラメータなど。素材となるデータは乱数で発生させているので、気に入ったデータが出来るまで何度も生成の繰り返し。これらの作業がけっこう大変！ 多くのパラメータはデフォルト値で済ませ、適当なところで「これでいいや」としました。
足し算・掛け算を変えてみたら
思いのほか手間がかかったので、2つのアニメーションを作っただけですが、当初はもっとたくさんのアニメーションを並べるつもりでした。
同じグラフ上に同じ初期状態を準備しても、遷移行列を変えれば現象（系の時間発展）は当然に変わります。さらに、遷移行列さえ同じであっても現象を変えることができます。行列計算のベースはスカラーの足し算と掛け算ですが、その足し算・掛け算を取り替えてしまうのです。行列係数の変更が量的変化だとすると、基本演算そのものを変えるのでよりドラスティックに質的な変化を引き起こします。あるいは、別な世界の現象をシミュレートすることになります。
通常の足し算・掛け算とは違う演算を持つスカラー系としては、一連のエキゾチック代数があります。これらのどれかを選んで、足し算・掛け算を入れ替えるとエキゾチックな世界の現象を扱えます。
マスロフ脱量子化では、通常の足し算を、連続パラメータ（プランク定数）に沿って変化させて、その極限としてmax-times代数（足し算が「最大値を取る演算」になったスカラー代数）に至ります。プランク定数hの値として、h = 1 以外に h = 0.5, h = 0 だけでもシミュレーションして、それらの挙動を比較できたら面白そう。
残念ながら、Rの関数／パッケージは通常のスカラー代数しかサポートしてないし、計算アルゴリズムがジェネリックに書かれているわけでもありません。エキゾチックな計算は自前で準備する必要がありそうですが、計算と可視化の環境は整備されているので、出来なくはないでしょう（やるとは言ってない）。
プログラミング言語としてのRを（ちょっとだけ）使った感想は、また別の機会に述べます。
[追記]いつも、サンプルのソースコードを記事内に貼り付けるようにしているのですが、今回は余りにもヤッツケでアニメーションを作ったので、少し整理してから貼ります。そもそも、Rの名前の付け方とかインデントの習慣とかも分からない。設計が古い言語なので、名前の管理がイマイチな気がするけど、僕の勘違いで何かうまいやり方があるのかも知れません。まー、そのうちいずれ。[/追記]
[追記 date="2015-04-28"]「Rの名前管理とか」にソースコードがあります。[/追記]
