マンデルブロ集合をOpenCLで描く(Mathematica)


CUDAを使ってマンデルブロ集合の描画を3桁速くする方法を以前紹介したのですが、同じことをOpenCLでやってみます。NVIDIAのGPUを搭載していないノートPCを使うことが多くなった自分用のメモでもあります。書き換え方は、OpenCLLink プログラミングで紹介されています。

書き換えたコードは以下の通りです(RSSリーダーでは表示されないかもしれません)。

最初に作成したコードをCore i7-3612QM 2.1GHz上で実行するのと比べて、ここで作成したコードをIntel HD Graphics 4000上で実行すると、3桁くらい速くなります。

Manipulateを使ってインタラクティブに描くコードは以下の通りです(RSSリーダーでは表示されないかもしれません)。

似たような話がMathematicaのマニュアルにも載っていますが、ここで書いたコードなら、描画の中心もインタラクティブに変わります。

realtime

CUDAを使ってマンデルブロ集合の描画を3桁速くする方法


追記:Mathematica 10で導入されたMandelbrotSetPlot[]は、GPUを使いませんがかなり速いです。この記事でやっていることは、以下のコードで再現できます。

追記:OpenCLの場合

描画が3桁速くなると、マンデルブロ集合をインタラクティブに楽しめるようになります。

はじめに

以前、Mathematicaでマンデルブロ集合を描く方法を紹介しました(「追悼、マンデルブロ博士」)。

Clear@f;
f[c_] := f[c, 0., 0]
f[_, _, 100] = -1;
f[c_, z_, i_] := With[{x = z^2 + c}, If[2 < Abs@x, i, f[c, x, i + 1]]]
DensityPlot[f[x + y I], {x, -2, 1}, {y, -1.5, 1.5}, PlotPoints -> 100]

描くのは確かに簡単だったのですが、「ここをもうちょっと拡大したい」という要望に応えるためのインタラクティブな仕掛けは、遅すぎて使えませんでした。「追悼記事なのに・・・」という思いがあったので、ちょっと改善することにしました。

目的

計算を高速化して、インタラクティブなマンデルブロ集合描画システムを作ります。

手法

(1)計算自体(上の例のf)高速化と(2)並列化で高速化します。(1)のために、再帰を反復に書き換えます。(2)のために、CUDAを使います。マンデルブロ集合の計算は、いわゆる「embarrassingly parallel」というやつなので、CUDAの練習に向いていると思います。

この記事のようにインタラクティブなものではありませんが、Mathematicaのマニュアルにもマンデルブロ集合の例が掲載されています。

条件

時間を計る条件は次のようにします。「有理数だと遅い」と思われるかもしれませんが、浮動小数点数を使うと、TableとParallelTableの結果が違ってしまう危険があることへの対応です。ちなみに、これはバグだと思うのですが、Wolfram社の回答は「バグでは無い」でした。

xMin = -2;
xMax = 1;
yMin = -3/2;
yMax = 3/2;
dx = (xMax - xMin)/300;
dy = (yMax - yMin)/300;

計算機環境(2013年3月更新)

  • CPU:Core i7 950(HTオフ)
  • 主記憶:12GB
  • GPU:NVIDIA GeForce GTX 580
  • OS:Windows 7 64 bit
  • Mathematica 9.0.1

結果

比較の基準

最初のバージョンをCUDAを使うバージョンで、計算時間を比較してみましょう。

最初のバージョンの実行時間を次のようなコードで計ると11.6秒でした。

CUDA

CUDAを使うバージョンの実行時間を次のようなコードで計ると0.003秒でした。4000倍くらい速くなることがわかります。

まあ、GPU無しの比較対象は最初に思いつくものですから、「まずGPU無しでどこまで速くなるのか追求してから」と思う方は、下の「おまけ」を読んでください。GPU無しでも最初のものより80倍くらいは速くできるので、それと比べれば、「GPUの効果は50倍くらい」というのが妥当なところでしょう。

うまくいかないときは、CUDAのための環境がちゃんとあることを確認してください。Visual StudioとCUDAがインストールされていればいいはずです(OpenCL版の方が簡単かもしれません)。私は、Visual Studio 2008 Professionalのデフォルトではx64のコンパイラがインストールされないために少しはまりました。本稿執筆時点でMathematicaが使うデフォルトのCUDA SDKは3.1で、これがサポートするのはVisual Studio 9 (2008)までだということにも注意が必要です(カスタマイズすればいいのでしょうが)。追記:2011年8月にCUDA SDK 4に対応したので、Visual Studio 2010がそのまま使えるようになりました。

リアルタイム描画

描画には、ArrayPlotを使うことにします。ListDensityPlotのほうが高度ですが、少し遅いのと、座標の扱いに問題があったため使えませんでした。Imageは速そうですが、やはり座標の扱いに問題がありました。

ユーザインターフェースはMathematicaのManipulateで構築します。

realtime

画像上をクリックすると、その場所を中心にして再描画します。Scaleで描画領域の広さを、Log[maxSteps]で反復の上限を変更できます。

ぬるぬる動くというわけにはいきませんが、「ここをもうちょっと拡大したい」という要望には十分応えられると思います。もっと速くしたい場合は、ArrayPlotをどうにかするということになるのでしょうが、もう少しコードを書かなければならないでしょう。

ところで、さすがにこれってCDFにしてもだめなんですよねえ。

おまけ(細かいことなので、これ以降は読まなくていいです)

追悼記事がダメだった理由

追悼記事がダメだった最大の要因は、DensityPlotにあったようです。DensityPlotは、下に示すように場所ごとにメッシュの細かさを自動調整してくれる便利な関数なのですが、この内部が並列化されていないため、あまり速くならないのです。

densityplot

CUDAなしでどこまで行くか

コンパイル+並列化

Mathematicaでは、機械精度の数を使うことに限定すれば、ユーザ定義関数をコンパイルできます。最初のバージョンをコンパイルし、並列化して実行すると30倍くらい速くなります。

fc = Compile[{{c, _Complex}, {z, _Complex}, {i, _Integer}},
   If[i == 100, -1, 
    With[{x = z^2 + c}, If[2 < Abs@x, i, fc[c, x, i + 1]]]], 
   CompilationTarget -> "C"];
DistributeDefinitions[fc];

time = AbsoluteTiming[
   result = ParallelTable[ fc[x + y I, 0, 0],
      {y, yMin, yMax, dy}, {x, xMin, xMax, dx}];][[1]];

time0/time

28.44619
反復

再帰を反復に書き直すともっと速く、80倍くらい速くなります。

hc = Compile[{{z, _Complex}},
   Module[{c = 0. + 0. I, i = 0},
    While[i++ < 100,
     c = c^2 + z;
     If[2 < Abs@c, Return@i]];
    Return@-1]];
DistributeDefinitions[hc];

time = AbsoluteTiming[
   result = ParallelTable[ hc[N@x + N@y I],
      {y, yMin, yMax, dy}, {x, xMin, xMax, dx}];][[1]];

time0/time

77.13116

ちなみに、この例ではコンパイル時に「CompilationTarget -> “C”」を付けてもあまり効果は無いようです。

CUDA無しでのリアルタイム描画

この方法を使ってリアルタイムに描画してみます。スケールや反復の上限もリアルタイムに変更できるように、少し修正します。

hc2 = Compile[{{xIndex, _Integer}, {yIndex, _Integer},
    {maxSteps, _Integer},
    {width, _Integer}, { height, _Integer},
    {xMin, _Real}, {xMax, _Real},
    {yMin, _Real}, {yMax, _Real}},
   With[{z = xMin + (xMax - xMin) xIndex/width
       + (yMin + (yMax - yMin) yIndex/height) I},
    Module[{c = 0. + 0. I, i = 0},
     While[i++ < maxSteps,
      c = c^2 + z;
      If[4 < Abs@c, Return@i]];
     Return@-1]],
   CompilationTarget -> "C"];
DistributeDefinitions[hc2];

width = 300; height = 300;
Manipulate[
 With[{
   xMin = center[[1]] - 10^scale, xMax = center[[1]] + 10^scale,
   yMin = center[[2]] - 10^scale, yMax = center[[2]] + 10^scale,
   maxSteps = IntegerPart[10^logMaxSteps]},
  result = ParallelTable[
    hc2[xIndex, yIndex, maxSteps, width, height, xMin, xMax, yMin, yMax],
    {yIndex, 1, height}, {xIndex, 1, width}];
  ArrayPlot[result, DataRange -> {{xMin, xMax}, {yMin, yMax}},
   DataReversed -> True, ColorFunction -> "Rainbow"]],
 {{center, {-0.5, 0}}, Locator},
 {{scale, 0, "Scale"}, -13, 1},
 {{logMaxSteps, 2, "Log[maxSteps]"}, 1, 3}]

実行結果は上に掲載したものと同じです。速いCPUがある場合はこれで十分なのかもしれませんが、私の環境では、「一応動くがストレスフル」という感じでした。

数独の平凡な解法(C言語)


数独の解き方をいくつか紹介したことがあります。

困ったことに、「どれも変態でわかりにくい、ふざけんな」という批判を受けたので、みんな大好きC言語で、平凡に書いてみました(実行結果)。

数独の平凡な解法(Mathematica)

Excelを使っていると採点ミスで謝罪することになる?


事例が見つかるたびに騒ぎになっているようです。今回は、「8.2 – 7.2」の結果が1より小さくなってしまうというものです。最初に非難されたのは例によってExcelのようですが、他の環境でもこういう問題は起こります(C言語の例Rubyの例)。

#include <stdio.h>

int main() {
  double x = 8.2 - 7.2;
  if (x >= 1) printf("x >= 1");
  else printf("x < 1");
  return 0;
}

「コンピュータは数値を2進数としてしか表せなくて・・・」という説明をよく見ますが、この説明は間違いです。これでは、「コンピュータは正しい計算ができない」という誤解を生んでしまいます。

10進数の1, 2, 3, …を2進数の1, 10, 11, …で表すという具合に、数を2進数で直接表現しようとするとこういう問題が起きますが、別の表現方法を用いれば、このような問題は起きません。実際、Mathematicaは正しく計算します。

With[{x = 8.2 - 7.2}, If[x >= 1, "x >= 1", "x < 1"]]

x >= 1

コンパイルしても大丈夫です。(UMMでは動かないようです。)

f = Compile[{a, b}, With[{x = a - b}, Print[If[x >= 1, "x >= 1", "x < 1"]]]];
f[8.2, 7.2]

x >= 1

C言語やRubyではどうすればいいのかというと、許容可能な誤差を決めておいて、その範囲内はよしとするのです。

さらに、この許容可能な誤差は、扱う数によって変えなければならないことに注意してください。以下の例では、許容可能な誤差を1より大きくしなければなりません(実行結果)。

#include <stdio.h>

int main() {
  double a = 9007199254740992;
  double b = a + 1;
  if (a == b) printf("a==b");
  return 0;
}

4891006269こういうことは、『Microsoft Visual C++入門』の3.2.11項「浮動小数点についての注意」にちゃんと書いてあったりします。ですから、入門を終えたプログラマはこういうことは知っています。

Excelのような、こうことは知らないであろう人が使うソフトウェアでは、こういう工夫をしなくても正しく計算できるようにしておいてほしいものです。フェルマーの最終定理の「反例」のような冗談を言っているうちはいいのですが、「x >= 0.6 なら合格」というようなことって、よくやられているのではないでしょうか(ボーダーの学生が間違って不合格になる危険があります)。間違いが公になるときには、「採点ミス」になるのでしょうか。

世界一難しい数独


解けたら天才! フィンランドの数学者が作った「世界一難しい数独」というのがあるらしい。

こちらの問題、実はフィンランドの数学者Arto Inkala氏が作成した問題で、曰く「世界でもっとも難しい数独」とのこと

1847534511この問題の作者であるArto Inkalaさんは、これまで「史上最難」と言われいた問題「AI Escargot」の作者でもある。その彼が言うのだから、間違いないと言いたいところだが、そもそも、数独の難しさに客観的な指標なんてないはず。合理的なものに限っても、たまたま知っていた定石で一発、なんてこともあり得る。だから、「ある解き方においては最も難しい」という限定的なことしか言えないのだ。彼自身、著書“AI Escargot”の中で次のように述べている。

it is often impossible to clearly identify the clearly most difficult puzzle in the world.

それにも拘わらず「世界でもっとも難しい」なんて言うことの大人の事情は放っておいて、せっかくだから手元にあるアルゴリズム(単純なしらみつぶし)でどうなるかを見てみよう。解を見つけるまでにどれだけの候補をチェックしなければならないかを調べる。Mathematica版なら全体をTraceで囲んで、結果の長さを見ればよい。

problem={
  {0, 0, 5, 3, 0, 0, 0, 0, 0},
  {8, 0, 0, 0, 0, 0, 0, 2, 0},
  {0, 7, 0, 0, 1, 0, 5, 0, 0},
  {4, 0, 0, 0, 0, 5, 3, 0, 0},
  {0, 1, 0, 0, 7, 0, 0, 0, 6},
  {0, 0, 3, 2, 0, 0, 0, 8, 0},
  {0, 6, 0, 5, 0, 0, 0, 0, 9},
  {0, 0, 4, 0, 0, 0, 0, 3, 0},
  {0, 0, 0, 0, 0, 9, 7, 0, 0}};

Reap@Length[Trace[search[{problem}, Join[#2, #1] &, False]]]

{80070, {{1 4 5 3 2 7 6 9 8}}}
          8 3 9 6 5 4 1 2 7
          6 7 2 9 1 8 5 4 3
          4 9 6 1 8 5 3 7 2
          2 1 8 4 7 3 9 5 6
          7 5 3 2 9 6 4 8 1
          3 6 7 5 4 2 8 1 9
          9 8 4 7 6 1 2 3 5
          5 2 1 8 3 9 7 6 4

となるので、この問題をMathematicaで解くときの、Traceの長さは80070。

一方、AI Escargotはと言うと、

problem={
  {1, 0, 0, 0, 0, 7, 0, 9, 0},
  {0, 3, 0, 0, 2, 0, 0, 0, 8},
  {0, 0, 9, 6, 0, 0, 5, 0, 0},
  {0, 0, 5, 3, 0, 0, 9, 0, 0},
  {0, 1, 0, 0, 8, 0, 0, 0, 2},
  {6, 0, 0, 0, 0, 4, 0, 0, 0},
  {3, 0, 0, 0, 0, 0, 0, 1, 0},
  {0, 4, 0, 0, 0, 0, 0, 0, 7},
  {0, 0, 7, 0, 0, 0, 3, 0, 0}};

Reap@Length[Trace[search[{problem}, Join[#2, #1] &, False]]]

{71766, {{1 6 2 8 5 7 4 9 3}}}
          5 3 4 1 2 9 6 7 8
          7 8 9 6 4 3 5 2 1
          4 7 5 3 1 2 9 8 6
          9 1 3 5 8 6 7 4 2
          6 2 8 7 9 4 1 3 5
          3 5 6 4 7 8 2 1 9
          2 4 1 9 3 5 8 6 7
          8 9 7 2 6 1 3 5 4

となる。解くときのTraceの長さ71766は、先の値80070より小さい。ということは、ここで採用しているアルゴリズムを基準にすれば、今回発表された問題は、確かにAI Escargotより難しい。

しかし、このアルゴリズムを基準にするなら、もっと難しい問題はある(Mathematicaでは時間がかかるからC++版で試す。オプション「-pg」を付けてコンパイルして実行し、gprofでプロファイルを見ればよい)。MathematicaでもC++でも、わざわざコードを書き換えたりしないのがポイント。

要するに、どの隅から解いていくかという話なのだが、ここで使っている方法の場合、どちらの問題も、右下から左下に向かって解くときに最も難しくなる。再帰的な関数searchが、今回発表された問題180度回転したものの場合は110874回、AI Escargot180度回転したものの場合は1502267回呼ばれる(どちらもCodepadでは解けない)。

つまり、ここで採用した基準では、世界一難しいのは今回発表された問題ではなく、依然としてAI Escargot(を180度回転したもの)だということになる。

追記:どちらも自力で解けた友人によると、AI Escargotの方が難しかったそうです。もう一つの方は30分ほどだったとか。