Zenn
🐠

反応する反応拡散系

2022/08/25に公開
#
CUDA
#
インタラクティブ
#
反応拡散系
#
セルオートマトン
#
tinkerboard
tech

(2017年10月〜2018年1月ころ) reaction-diffusion (反応拡散系) オートマトンのGray-Scottモデルを、いろんなマシン （PCとか単体専用機） でうねうね描いて、インタラクティブおもちゃを作ってみたはなし。

ソース:

https://github.com/tarohs/rd

※ 2022年8月現在、ソースみながら思い出して書いているところもあるので、やや不正確かもしれませぬ。

実行結果 （ARM Macbook Pro, SDL2, 200 x 120):

https://youtu.be/VdHEHO51TY4

※ 写真や画面キャプチャは、おいおい追加します。

反応する反応拡散系

ことのはじまりが何だったかはよく覚えていないのだが、反応拡散系というセルオートマトンをデモするプログラムを作りたいと思った。
横浜国大の数学の今野先生が『数理芸術』という催しを個人的に開催されている。数学に関連する制作物……一発芸的なものから「へぇ〜」と感心するような工芸まで、主に研究室の学生さんが作って展示している。横浜元町のおしゃれなカフェで数日間展示して、最後はパーティーで盛り上がる、という、なかなかナイスなイベントなのだが、知り合いを通じて2016年ころから、わしの3Dプリントや電子工作などを披露させてもらっている。
んで、その年のネタは、ちょーどハマっていた反応拡散系をデモることにした。

マシンはラップトップPCとか、デスクトップPC + 外付けディスプレイではなく、単体機であることを目指した。
セルオートマトンの単体機 (専用機) ってのがすでに意味不明だけど、まあそういうおもちゃというかオブジェだと思っていただければよろしい。
あと専用機なので、画面にタッチするとインタラクティブに状態を変えたり、傾けると反応するようにしたい。音なんかも出るといいな。

つまり、リアクション (反応) する反応拡散系、である。
ということで、いろいろ作ってみた記録である。

反応拡散系 (Gray-Scottモデル) についてちょっと

反応拡散系というのは、有名なライフゲームみたいなセルオートマトンである。

ただしライフゲームの状態が生き死にの2値であるのに対し、作ったGray-Scottモデルの反応拡散系というのは実数値を取る2次元のプレーン2枚u, vである。またライフゲームが周囲8近傍の生の個数を数えて次の状態を決定するのに対し、Gray-Scottモデルでは周囲4近傍のu(x\pm 1, y\pm 1), v(x\pm 1, y\pm 1)の実数値から次のu(x, y), v(x, y)の実数値を決定する。

これが計算機で動く人工生命っぽいということで、かのアラン・チューリングさんもハマって、いろいろ貢献したりしている。実際にわしが取り上げたGray-Scottモデルってのは、2種類の薬物 (おい) が隣の細胞から染み出してきて、ある細胞内で反応してまた隣の細胞に散っていく、という意味があるらしく、与える定数によってはシマウマの模様 (発生のときに白黒のシマがランダムにできる) とか熱帯魚の模様とかサンゴとか、それっぽい模様ができる。

オリジナルのGray-Scottモデル

https://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/

は、連続では次のように定義されるが

\frac{\partial_u}{\partial_t} = r_u \nabla ^2 u - u v ^2 + f(1 - u)\\ \frac{\partial_v}{\partial_t} = r_v \nabla ^2 v - u v ^2 + (f + k)v

(f, k はてきとーな定数)、離散すなわち計算機上でのマス目を変化させるセルオートマトン （ライフゲームみたいなの） ではよーするに

u(x, y) <= u(x, y) + dt * (du * (
                           u(x, y - 1) + u(x - 1, y) +
			   u(x + 1, y) + u(x, y + 1) - 4 * u(x, y)
			         ) - u(x, y) * v(x, y) * v(x, y)
				   + f * (1 - u(x, y))
		          ) 
v(x, y) <= v(x, y) + dt * (dv * (
                           v(x, y - 1) + v(x - 1, y) +
			   v(x + 1, y) + v(x, y + 1) - 4 * v(x, y)
			         ) - u(x, y) * v(x, y) * v(x, y)
				   + (f + k) * v(x, y)
		          )

(dt, du, dvは更新する時間間隔や解像度に応じた定数 (1でよい)、f, kはてきとーな定数) を、画面上のセル(x, y)全域に渡って計算するだけである。目がちらちらすると思うが、前半はよーするに上下左右の和からその位置の数値を引いて、あと何倍かして現在の値に足しているだけである。
プレーンは2組用意して、状態を更新するまで前のプレーンを変更しないのがミソなのは、ライフゲームと同じだ。

画像へのマッピング

得られたu(x, y), v(x, y)に応じて、カラーで画素(x, y)にマッピングすれば、なかなか綺麗なグラフィックスが得られる。
マッピングはいろいろ試してド派手かつちょっとキモいw 感じになるように、下記のように決定した。

そのときどき、あるいはパラメータによって、u(x, y)とかv(x, y)が取りうる実数値は異なるので、まずu, v全体を見渡して、最大値と最小値を得る。
個々のセルのu(x, y), v(x, y)の値が最小値〜最大値の間でどのくらいかを0〜255にマッピングする。uは赤、vは青、緑は赤と青の差の絶対値、としている。

定数の決定

定数f, kをてきとーに変えてやると、シマウマ・熱帯魚・珊瑚などいろいろな模様ができる。定数によって安定したり不安定だったりする……つまり、ある状態に収束してしまって動かなくなる (模様が安定したり、逆に不安定で世界全体が消滅したり) とか、逆にいつまでも不気味にぷるぷるしていたり、ダイナミックに模様が変化することもある。

f, kについては下記あたりで詳しい。

http://mrob.com/pub/comp/xmorphia/uskate-world.html

トップの海岸マップっぽい画像は、fとkをどー変化させるとどーいう模様になるか、みたいなマップになっているので、実際にマシンに流し込む定数は、これを参考にさせてもらった。

初期状態とか

クリーンな状態では、uは全部のマスが1、vは0にしている。
その上にランダムに『アクティブセル』みたいなのをばらまきたい。
ライフゲームなんかは全面クリーン状態の0 (死) の上にてきとーに1 (生) をばらまけば、まあ良い感じの密度に収束して、あと死滅したり安定したりフリップフロップやグライダーとか動いたりしちゃうのだが、Gray-Scottモデルでは初期にばらまくものが点ではダメなようである。u = 0.5, v = 0.25とかの値を一定の大きさの正方形の塊にしとかないと、すぐに消滅してしまう。

また上記の定数f, kによっても、一個アクティブな正方形があればそこから段々模様が増えていく定数だったり、かなりの数のアクティブをばらまかないと死滅してしまう定数とか、アクティブな正方形がトリガとなってダイナミックにループがぐわんぐわん動くなんてのもあるので、初期状態でどのくらいの個数・密度の正方形をばらまいておくか、は定数f, kによっても異なる。

後述するように、インタラクティブなマシンに仕上げたときには、途中で定数f, kのセットを変えたり、指でなぞった軌跡の通りアクティブ値を書き込んだり、マシンの傾きによって場所によって定数f, kを変えてみたり (傾いている画面の、上の方が不安定で下のほうが安定にすると、模様が傾けた下の方に流れていくようになる) した。

作ってみたいろいろなバージョンたち

これらはすべて、ひとつのツリー内で#ifdefで切り分けたりしているので、ソースはやや煩雑になっている。

※ 後述の2021年『デムパ時計』に搭載したreaction & diffusionルーチンと、2022年になってからESP32・Teensy + 超小型のOLEDや液晶で超小型単体機を作って『マイコンの速度に挑戦』したりもした。これらは別ソースにしてあるので、いずれまた別の記事で。

プロトタイプ: PC上でのエミュレーション

さいしょは、Python3 + OpenCV3で実装するのが簡単だろ、と思って作った。まあこれはすぐできる。
ついでに、マウスやキーボードで動作中に模様を描いたり (一定の大きさの四角形を落とす) 定数のパラメータセットを変えたりできるようにしてみた。
Python3は……こういうiterativeなセルオートマトンとか、ほんっっっっとに遅いですね(笑)。

Gray-Scottモデルは、定数にもよるが、高速で計算すると模様がじわ〜っと拡がるか、またはダイナミックなうねうねアニメーションになったりもする。だから、なるたけ高速で計算したい。
単純な計算である、と書いたが、これがなかなかの計算量である。実数であるうえ、セル全体を順番に走査して真面目に計算しないといけない (当たり前だ)。その上、ライフゲームみたいに1ステップごとにちゃきちゃき更新されるわけではなく、おおむね100ステップくらい計算して「うねうね」の「う」の動きがみられる程度である (描画に計算より時間が掛かる場合は、数〜数十ステップの計算ごとに画像1フレームを作っても、それほど動きが飛んだりしないほどである）。
Python3でプロトタイプを書いた当初は、プログラムの誤りというよりは、前述の初期値として点をばらまいてしまったり (すぐ消滅する)、遅すぎて「なんか更新されてない」と思ったり、で正しく動いているようにみえなかった。

ということで、すぐC (後でいろんなライブラリを使う都合からC++っぽく書いてはあるけど) に書き直すことになる。

ただし、目標は単体機の (能力が低い) マイコンでの計算で、ディスプレイも解像度が低いのがくっついているから、少なめのドット数でどう見えるか、というシミュレーションをしたい。セル数は32\times 32 (拡大率10倍で画面表示する)、64\times 64 (拡大率8)、128\times 128 (拡大率4)、252\times 252 (拡大率3)、200\times 120 (拡大率4)、256\times 178 (拡大率4)といったプリセットから選べるようにしておいた。任意のx, yも与えられるが、後述するSIMDのおかげで4の倍数だとかの制約がある。
つか、解像度が画面があるからといって、例えば4096\times 2160ドットとかまじめに計算したら爆速コンピュータでも相当低いFPS値になっちまうだろうし、模様が細かすぎて観ていたら頭が痛くなるに違いない。

実は最初に単体機の表示装置として考えていたのが、32\times32ドットのマトリクスLEDである。
HUB75という6ビットくらいのパラレルインターフェースを備えていて、RGBのLED自体はON/OFFだが、超高速のマイコンやFPGAでPWM制御すれば中間色も出せる (これを連接したのが、球場の表示パネルだったりする)。速度の関係 (つまり膨大な計算をしながら超高速PWMで表示もこなす)、および32\times 32セルというのが少なすぎる (ほとんどのパラメータでは、周辺およびリングワールドの反対側の影響 (後述) を受けて、十分面白い変化にならない)、といった理由で途中で挫折したのだが、はからずもこれは2021年に実現している (『デムパclock』の製作記事参照)。

RaspPiへの実装

単体機といえば!! やっぱりRaspberry Piシリーズであろう。
これで作るのは、Arduinoマイコンなどと違って実はまったく難しくない。OSとしてLinuxが載ったただのPCだし(笑)、ボードにハット形式で一体化する専用のディスプレイも発売されているし、それがタッチパネルになってたりする。つまり、OpenCVなりOpenGLなりをapt getしといて、ただのC++とかPythonの『アプリケーションソフト』を書けば (電源オン→自動ログイン→自動実行時→ウィンドウを最大化すれば)、それで『専用機』は完成してしまう。キーボード入力やマウスでの描画(?)に反応するプログラムにすれば、ちゃんとリアクションする反応拡散系になってくれる。

で、つくりはじめた当時はRaspPi 3Bだったのだが、すぐ後にASUSの爆速互換機 (? コネクタの位置などは同じ。電源がUSBでなく通常のACアダプタだった) Tinker Boardを入手して、それで最終形態の単体機を組むことにした。SoCがRockchip RK3288というCoretex A17系のARMで、GPUでMali T764が載っている (ちなみにこのGPUもハックしてやろうと思ったが、資料が十分に手に入らなかった)。
タッチパネルディスプレイと一体化してテーブル様の台を付けるために、ちょろっとアクリル工作をする。

CUDAとマルチコア

このころから高速化に目覚めるようになる。
つまり……いくらラップトップとはいえ、4コアだか6コアのGHzクラスのCPU (当時Intel) を積んでいるMacbook Proでプロトタイプを開発していたので、RaspPiや互換機のTinker Boardが遅い遅い。

で、どうにかならないか……と思っていたら、手元に転がっていたのが、NVIDIAのSBC、Jetson TK1である。
これまた、CPUなんかはRaspPiとかとほとんど同じARMマイコンなのだが、CUDA GPUを積んでいるのが違う。
すぐ後にJetson TX1も使えるようになるのだが、こっちはちょっと大きめである。
TK1もTX1も、外付けディスプレイやキーボード・マウスを必要としたので、まあこれは参考程度というか、GPUプログラミングへのチャレンジとして移植してみた。

まあでも、CUDA搭載とはいってもCUDAコア数が256とかのマイコン、というのは、セルオートマトン単体機 (だからなんじゃそりゃw) としては立ち位置が微妙である。つまり、ラップトップのMacbook ProとかでCPU演算してる方がまだ高速で、単純に爆速を目指すなら2048 CUDAコアとかのNVIDIA GTXやRTXシリーズを挿したデスクトップ用のGPUカードのほうがよい。
そもそもJetsonシリーズは単体でdeep learningなどの高速計算をせざるを得ない機器のために開発されている。自動運転車や障害物を避けるラジコンマルチコプターにMacbook ProとかタワーPCを搭載して走らせたり飛ばしたりするわけにいかないので (いやTeslaとか積んでるかもw)、Jetsonの存在意義がある。
別に単体機で動く必然性も需要もない、セルオートマトンを表示するだけのおもちゃのためにあるわけではない。

さておき。
どーしてCUDAだと (つまり単純計算の超並列化) 速くなるかというと、画面が分割できるから・分割したそれぞれを並列に計算できるから、である。
当たり前である。
当たり前だけど、どう分割するかなどをうまく考えてやらないと、それほど速くならない。
例えば隣り合ったセルをすべて違うCUDAコアに割り当てちゃったりすると、上下左右のセルの値が入った共有メモリをそれぞれが参照するわけだから、ぶつかってそれほど速くならなかったりする。

まあここはちょっと悩んでいろいろ試しましたけど、最終的には単純に、連続した一定の領域に等分してますな。あと、後述する斜め演算も併用できるので、それもやってる。

  • CUDA nvccでコンパイルするときは、CMakeLists.{TX1,TK1}.txtを使ってください。

IntelのMacbook Proでも、4だか6だかのマルチコアなので、CUDAほどばきばきにプレーンを割るわけではないが、いちおうstd::threadで並列化して速くなっている (#define THREAD)。

もっと速く・もっと計算量を減らせ

GPUを使うと使わないとに関わらず、CPUでの計算が速いか遅いかに関わらず、計算の工夫で多少速くなることはある。

まず気づいたのが、意外にコンパイラ (RaspPiのGCCとかmacOSのLLVM) がアホであること。
例えば-Oxすれば、オリジナルの定義式通り (ソースを読みやすく) 書いておいても、ループで変化しない値なんかはとっといてくれるのかと思ったら、いちいちループ中で計算してたりするので (例えば dt * du / dx / dxなんてのは変化しない値だが、定義式通り書くといちいち計算しやがる)、自分で分配法則適用して括弧展開したほうが速い(笑)。

あと『上下左右の和から中央のマス目の4倍を引く』なんてのも、工夫すれば速くなる (ソースiter.ccで#define OPTEQ)。つまり、あるセルの上と左の和を計算したプレーンを作っておけば、『上下左右のマス目の和』は『そのマス目の上・左の和』プラス『右斜下のセルの上・左の和』で計算できる。メモリの使用量は2倍になるが、セルひとつに付き加算が2回減る……だけではなく、これがGPUやスレッドなんかだと共有メモリの競合が減るためか?、結構 (※ 当社比) 速くなる。

あとはSIMD、つまりパイプライン命令だ。
float 32bitでの単純な加減乗除くらいなら、4つの値 (4つの連続するセルの値) を並べて、128bitのALUでえいやっ!! と一発で計算する、という仕組みがある。
Intel (Macbook Pro --- 当時) ではSSE命令、ARM (RaspPiや互換機・Jetson) ではNEON命令というのが使えたので、定義式 (を手動で展開・最適化したもの(笑)) を睨んで、そのへんが使えそうな行をアセンブラで書き直したりしている (SSEは#define SSE、NEONは#define NEON)。

まあ、GPUの計算手順やブロックサイズ、どこをSIMD化するか、なんてのはさじ加減なので、この辺は実際に動かして試行錯誤で一番速いところをみつけたりしている。

さらに描画に関しては、初期のOpenCVから、OpenGLを使ってみたり (Jetson、#define USEGL)、SDL2を使ったり (macOS・RaspPi、#define USESDL) して高速化している。確か記憶では、OpenGLはそれほど速くならなかったような……あとSDL2は、後述の音を出すライブラリの関係もあったと思う。

世界の終わりとリングワールド

ここはちょっと、セルオートマトンそのものの挙動に関わることなのだが。
世界、つまりプレーンの端っこの処理をどうするか。

ここが反射だ吸収だ、ってのは数学的にも興味ある問題らしいのだが、とりあえず計算機、ってか限られた画面上での実装は、端っこをどうするかの問題を解決しないと、そもそも配列からはみ出してしまう。
いちおう、最終的には動作中に吸収端とドーナツ世界を切り替えられるようにした。

吸収端は、端の一個外側に全0の値のセルがあると仮定する。これは上下左右に1セルだけ大きめのプレーンを取って0を突っ込んでおき、元のセルの部分だけを計算・更新すればよいので簡単だ。

ドーナツ世界では、上辺と下辺がつながっており、また左辺と右辺がつながっている。つまり、最上セルの上は最下セル (または逆)、最右セルの右は最左セル、ということにする。これまた、ライフゲームなどでは定番の手法である。
ただし計算を速くするため、『ねじれドーナツ世界』にすることにした。これは、2次元のセルを1次元のメモリ上で単純に考える。上下はxサイズだけ前または後ろをアクセスする。そうすると、最右セルの右は一行下の最左セル (またはその逆) となる。最上行の上に1行分だけ最下行のセルの値を、最下行の下に1行分だけ最上行のセルの値をコピーしといてから計算をはじめる (前述の斜め加算 (上と左) をする場合には、さらに考えなければならないが)。

反応拡散系の音

音を出したい。
なんでもいいから、画面の模様の変化によって変化する音が欲しい。

まあ時間があって、計算パワーがあれば、画面の色、っていうか(u, v)の値を音程や音色にマッピングしてMIDI楽器でも鳴らしたかった。模様パターンを、用意しておいたリズムパターンやフレーズにマッピングしてもよいなあ。
なんて妄想していたが、最終的には単純に、青の画素値を左上から右下まで、ステレオ外付けスピーカ (単体機だがスピーカだけは別売りだ) から出している。単純すぎ (笑)。

どういうことかというと。
セル数が200\times 120で、サンプリング周波数22050Hzなので、これを単純にスキャンすると、約1秒周期で「ぎっ ぎっ」みたいな音がでる。例えば画面に縦縞が1本あったとすると120Hzくらいの音になり複数の模様があるとその数倍の周波数になるわけで、模様がちょーど耳に聴こえる音になる。
それが模様が変化すると「ぎゃあ ぎゃあ」になってみたり「じゃーじゃー」になってみたり、なかなかそれっぽい。隣り合った画素をステレオLRに振り分けているので (というより、なんも考えなくてもそういうオーディオバッファになってるので)、なんとなく似たようななんとなくずれた音が左右から出て、立体感もばっちりだ(笑)。

※ トップ動画の音 (Macbook Proのオーディオデバイス） は、1chモノラルですね。すみません。

ちなみにオーディオはSDL2で簡単に扱うことができる (#define SDLAUDIO)。表示する青の値をオーディオバッファにも突っ込んでおいて、再生バッファが空に近くなるとかかるコールバックでオーディオバッファを転送するだけだ。手抜き(笑)。

反応する単体機

さて、単体機『反応する反応拡散系』のテーマは、周囲の環境、つまり手で画面に模様を描いたり、傾けたり揺らしたりに応じて、画面のうねうねが反応してほしい、ということである。

で実は、これは最初から織り込み済みなのだが、デバッグ用のPCでもキーボードとマウスで

  • マウスで画面に描画すると、そのセル周囲の正方形に(u, v) = (0.5, 0.25)を置く (初期にばらまく『アクティブ値』と同じである)
    • 初期状態では全セルクリア後、パラメータに応じてランダムにこの正方形が置かれる
  • キーボードで
    • 『0』『1』『9』パラメータセットの切り替え
    • 『w』で周縁部を吸収にするかリングワールドにするか切り替え
    • 『y』でアクティブ値正方形をランダムに置く
    • 『x』で画面クリア・『z』で衝撃を与える （画面上のすべてのセルの(u, v)を一定量いじる）……トップ動画の1分14秒あたり
    • 『場を傾ける』 （左右によって・上下によって画面上のfの値の分布が変わる --- abc...j(中央)...qrsで最左傾き→水平→最右傾き、ABC...J(中央)...QRSで最上傾き→水平→最下傾きというように）……トップ動画の45秒〜1分09秒あたり

などのインタラクションを与える機能をつけてある （トップ動画参照）。

マウス入力はタッチパネルでも同様にできるので、あとは衝撃や傾きによって、対応するキー入力をUSBポートから模倣してくれるデバイスを外付け （実際には単体機『内部』に外付けしている） を作ればよい。

これには、以前使ったことがある （パーツ箱に余っていた）、プログラムにも使うオスUSBポートがUSBホストにもなり、USB HID (keyboard)クラスを吐けるmicro beetleという基板 (ATmega32u4) を使った。
本体のLinuxマシンとmicro beetle上で動くArduinoプログラムは、それぞれ独立したコンピュータのプログラムとして開発すればよいし、デバッグ入力と同じ方法でこのふたつのコンピュータ間通信が行なわれているので、最も開発が楽なやり方である。

単体機を傾けるとうねうねの模様が下に下がっていったり、どかんと机に置くとセルが死滅したり、というのは、なかなか面白いものである。

実はマイク入力とかカメラ入力もつけてみたかったのだが、まあ同様の外付けデバイスを追加すれば、そんなに難しいことではなかろう。

展示とその後

まあ、『数理芸術』 (2017. 12. 10--) では、その道のひとたち （数学のプロ） がいろいろ興味を示してくれた。
何を思ったか、展示期間中にふつーの （Conwayの） ライフゲームを単体機に追加したり、Macbookのカメラから入力した顔画像とかを初期値にしてそっからうねうね……とかのルーチンを追加している（笑）。

ここまで1ヶ月位でいろんなアーキテクチャで、血を吐きながら （逆流性食道炎が悪化してほんとに胃から出血した……w。大したことなかったんだけど） 最適化とかGPUプログラミングとかアセンブラ書いたりとかしたのは、面白かった （なおJetson TX1バージョンや、Macbook Proで動くバージョンも参考に展示した）。
その後も2ヶ月くらい?、微妙なバグ取りやオーディオのレイテンシ （つまり模様の変化が生じてから発音の変化が聴こえるまでの遅れ……わかるわけないんだけど（笑）） を改善したりしている。

だがなんといっても、2022年現在すげーなと思うのは、後年このGray-Scottモデルのプログラムのコア部分が、別のハードウェアでまた動いていることだ。

ひとつめは、速度の問題で （があると思って） 挫折した、HUB75信号で駆動するフルカラーLEDマトリクスパネルへの表示が、2021年に製作した、『ボケをかますLEDクロック』こと『デムパ時計』の一アニメーションとして成功してしまったことだ。ドットマトリックスLEDパネルへの表示ルーチンを公開してくれたひとがいて、あとESP32が2コアのマルチコアである (計算と表示に1コアずつ振り分けている) ことが大きい。

もうひとつは、SPIやI2C接続の、いわゆるマイコン用の超小型OLEDやLCDパネルに、ESP32とかTeensyの速度のベンチマークとして実装し、そこそこうねうね動かしたことである。これまた面白い （ディスプレイとの通信など） 最適化をかましてしまった。

これらはまた別口で記事にしようと思う。

Discussion