OAuth認証でTwitterを利用するWebアプリケーション(PHPの場合)


OAuth認証が必要なAPIでは、Twitter APIとFacebook APIが有名ですが、ちょっと仕様が違うので、使い方を簡単に紹介しようとすると、記事を分けなければなりません。さらに、JavaとPHPの両方を学べるというお得な本を書いたこともあって、その読者をサポートするために、JavaとPHPの両方で記事を書かなければなりません。都合、以下の4パターンが必要です。

今回は、「PHPでTwitter API」です。「JavaでFacebook API」はまた別の機会に。

アプリの登録

http://dev.twitter.com/appsでアプリケーションを登録し、Concumer keyとConsumer secretを取得してください。その際、クライアントアプリケーションではなくブラウザウェブアプリケーションとして登録するように注意してください。

OAuthのためのライブラリ

PHPには、標準の「PECL/oauth」や準標準と言える「Pear HTTP_OAuth」、Twitterに特化した「Pear Services_Twitter」がありますが、例によって、PECLはWindowsで使いにくいし、HTTP_OAuthはベータ版、Services_Twitterはアルファ版なので、さっさと見限って「abraham / twitteroauth」を使いましょう。

OAuth.phpとtwitteroauth.phpをダウンロードして、ウェブサーバで配信できる場所に置きます。

cURL

XAMPP

PHPでcURLを使うために、php.iniに以下のように記述して、Apacheを再起動します。

extension=php_curl.dll

Ubuntu

次のコマンドを実行してcURLモジュールをインストールしてからApacheを再起動します。

sudo apt-get install php5-curl

4627847327お約束ですが、このあたりでもうついて行けないという場合は、拙著『Webアプリケーション構築入門』などを参照してください。以下のPHPプログラミングで必要な知識も本書にまとめてあります。

OAuth認証

以下の3画面で認証することにします(oauth-start.phpとoauth-end.phpは一つのファイルにまとめることもできますが、ここではわかりやすいように二つに分けています)。

  1. OAuth認証のためのURIを生成し、それにアクセスするためのリンクを表示するoauth-start.php
  2. アプリケーションを許可するTwitterページ
  3. Twitterからのコールバックを受信し、access tokenとtoken secretを取得するoauth-end.php

実装

oauth-start.php

OAuth認証の入り口となるoauth-start.jspは以下のようになります。「Consumer key」や「Consumer secret」、コールバックURIの部分は、状況に合わせて書き換えてください。

<?php
require_once('twitteroauth.php');

session_start();

//セッションを節約したい場合は別ファイルに書いておく。
$_SESSION['consumer_key'] = Consumer key;
$_SESSION['consumer_secret'] = Consumer secret;
$callbackUri = 'http://localhost/twitter/oauth-end.php';

$connection = new TwitterOAuth($_SESSION['consumer_key'], $_SESSION['consumer_secret']);
$request_token = $connection->getRequestToken($callbackUri);
$_SESSION['oauth_token'] = $token = $request_token['oauth_token'];
$_SESSION['oauth_token_secret'] = $request_token['oauth_token_secret'];

$authUrl = $connection->getAuthorizeURL($token);
?>
<!doctype html>
<html>
  <head>
    <title>OAuth Start</title>
  </head>
  <body>
    <p><a href="<?php echo $authUrl; ?>">Twitter OAuth認証開始</a></p>
    <p><a href="logout.php">logout</a></p>
  </body>
</html>

oauth-end.php

コールバック後の処理を実装するoauth-endは以下のようになります。

<?php
require_once('twitteroauth.php');

session_start();

$connection = new TwitterOAuth($_SESSION['consumer_key'], $_SESSION['consumer_secret'],
                $_SESSION['oauth_token'], $_SESSION['oauth_token_secret']);
$access_token = $connection->getAccessToken($_GET['oauth_verifier']);
unset($_SESSION['oauth_token']);
unset($_SESSION['oauth_token_secret']);

//セッションに記録する。データベースなどを使ってもよい。
$_SESSION['access_token'] = $access_token;
?>
<!doctype html>
<html>
  <head>
    <meta charset="UTF-8">
    <title>Twitter OAuth認証完了</title>
  </head>
  <body>
    <pre><?php print_r($access_token); ?></pre>
    <p>access_tokenをセッションに保存した</p>
    <p><a href="post.php">送信テスト</a></p>
    <p><a href="logout.php">logout</a></p>
  </body>
</html>

つぶやく

認証が終わったら、以下のようなコードでつぶやけます。

<?php

header('Content-type:text/plain; charset=utf-8');
require_once('twitteroauth.php');

session_start();

$access_token = $_SESSION['access_token'];
$connection = new TwitterOAuth($_SESSION['consumer_key'], $_SESSION['consumer_secret'],
                $access_token['oauth_token'], $access_token['oauth_token_secret']);

$message = "テスト at " . time();
$parameters = array('status' => $message);
$status = $connection->post('statuses/update', $parameters);
print_r($status);

ブラウザを再起動すればはじめに戻りますが、次のようなlogout.phpを作っておいてもいいでしょう。

<?php
session_start();
session_destroy();
?>
<!doctype html>
<html>
  <head>
    <title>Logout</title>
  </head>
  <body>
    <p><a href="oauth-start.php">start</a></p>
  </body>
</html>

Wolfram CDF PlayerをMathematicaとして使う方法


http://www.unfindable.net/umm3/

計算可能ドキュメント形式(Computable Document Format, CDF)を閲覧するためのソフトウェアWolfram CDF PlayerとMathematicaの関係は、Portable Document Format (PDF)を閲覧するためのソフトウェアAdobe ReaderとAcrobatの関係に似ています。Wolfram CDF Playerで閲覧可能なCDF文書を作るにはMathematicaが、Adobe Readerで閲覧可能なPDF文書を作るにはAcrobatが必要です。

しかし、CDFとPDFには大きな違いがあります。PDF文書は内容が固定された静的な文書であるのに対して、CDF文書は内容を変化させられる動的な文書です。下はCDF文書の簡単な例です。Wolfram CDF Playerがインストールされている環境なら、aの値を変えながらSin[a x]をプロットしてみることができます。CDF文書の内容は計算によって変化するのです。

Sin[a x]のaの値を変えられるなら、もっと大胆に「Sin[a x]」という式全体を変えられるのではないかと考えるのは自然でしょう。Mathematicaの式を処理できるCDF文書、それはMathematicaとして使えるCDF文書です。使い勝手は多少悪くても、「Mathematicaを使いたいけれど高すぎて買えない」という人にとっては有用でしょう。みんな大好きWolframAlphaも、Mathematicaのすべての機能を使えるわけではありませんし。

残念ながら、直接的な方法はうまくいきません。CDF文書に入力できるのは数値だけであり、「Sin[a x]」のような文字列は入力できないからです。しかし、コンピュータ上で表現されるすべてのものは、メモリ上では数で表現されています。「Sin[a x]」のような式ももちろん数で表現されています(メモリのことがよくわからない人は、ゲーデル数を思い出してもいいでしょう)。ですから、Mathematicaの式を一度数値に変換してからCDF文書に入力し、CDF文書内でそれを元に戻すというような工夫をすれば、CDF文書で式を扱えます。このアイディアを実現したのが、以前紹介したUniversal Mathematica Manipulator (UMM)です。

UMMには、Mathematicaの式を変換してできる数値が長大なため入力に手間がかかるという問題がありました。Wolfram CDF Playerには、文書上でクリップボードからの貼り付けができないという制約があるため、長大になる数値をCDF文書上で入力するための面倒な仕掛け(VBScriptやAppleScript)が必要でした。(Mathematicaに付属するCDF Playerなら貼り付けられます。このように仕様がばらばらなことが後で混乱を生まないことを祈ります。)

ここでは別の方法を紹介します。

まず、Wolfram CDF Playerをインストールします。これがなくては始まりません。

次に、PHPを使えるウェブサーバをlocalhostに用意します。WindowsならXAMPPを導入するのが簡単でしょう。

下のような内容のexpression.phpをhttp://localhost/umm3/expression.phpというURLで呼び出せるようにしておきます。ディレクトリumm3は、ウェブサーバから書き込めるようにしてください。

<?php

$file = 'expression.txt';

if (isset($_GET['body'])) {
  file_put_contents($file, $_GET['body']);
  echo $_GET['callback'].'()';
} else {
  if (is_file($file)) echo file_get_contents($file);
  echo '(* OK *)';
}

PHPのmagic_quotes_gpcがOffであることを確認します。XAMPPの場合はデフォルトでOffになっています。Macの場合はphp.iniを編集してApacheを再起動する必要があるかもしれません。

4627847327お約束ですが、上の手順がよくわからない人は、拙著『Webアプリケーション構築入門』などを参照してください。

以上の準備ができたらhttp://www.unfindable.net/umm3/にアクセスします。ピタゴラス3体問題のような比較的大きなプログラムも実行できることを、Wolfram CDF Player 8.0.4(Win, Mac)で確認しています。

円周率の近似値を手元に置くこと


B00264YQ80

円周率の近似値が印刷されたマグカップを手元に置くのが一つの手です。

487310002X

円周率の近似値が印刷された書籍、牧野貴樹『円周率1000000桁表』を手元に置くのも一つの手です。

必要な時にすぐ計算できるように、スーパーπをダウンロードしておくのも一つの手です。スーパーπはかなり前からある有名なプログラムで、ベンチマークなどでもよく使われています。本稿執筆時点で、419万桁の計算時間を集めたサイトスーパーπランキングのトップは、Intel i7 2600K@5.454GHz (Sandy Bridge)の35秒とのこと。私のマシン(Core i7 930)では121秒でした。

Mathematicaを使うのも一つの手です。Mathematicaが手元にない人は、UMMで代用してもいいでしょう。スーパーπのベンチマークと同じだけ計算したい場合は、「pi = ToString@N[Pi, 4194307];」とします。有効な結果はStringDrop[pi, -2]です(最後の2桁を捨てる)。私のマシン(Core i7 930)では16秒と、スーパーπよりかなり速く計算できました(もちろん計算結果は同じです)。

結果をそのままファイルに書き出しても4MB程度なので、「記録しておくより、必要な時に計算したほういい」という状況ではありませんが。

Audio Function Generator


4873114837Wolfram CDF Playerがインストールされていれば動きます(UMMでも動きます)。RSS Readerでは何も見えないかもしれません。

With[{functions = {Sin,
    TriangleWave[#/2/Pi] & -> "Triangle",
    SawtoothWave[#/2/Pi] & -> "Sawtooth",
    SquareWave[#/2/Pi] & -> "Square"}},
 Manipulate[Play[
   volume1 function1[2. Pi frequency1 t] +
    volume2 function2[2. Pi frequency2 t] +
    volume3 function3[2. Pi frequency3 t],
   {t, 0, 3}, SampleRate -> sampleRate],
  {function1, functions},
  {function2, functions},
  {function3, functions},
  {{frequency1, 440}, 20, 8000},
  {{frequency2, 550}, 20, 8000},
  {{frequency3, 660}, 20, 8000},
  {{volume1, 0.5}, 0, 1},
  {{volume2, 0}, 0, 1},
  {{volume3, 0}, 0, 1},
  {{sampleRate, 8000}, 1000, 44100}, ContinuousAction -> False]]

ボリュームは相対評価です。

バットマン方程式


バットマン方程式というのがあるそうです(Is this Batman equation for real?)。なんでも、この式を満たす(x, y)をプロットするとバットマンのマークになるのだとか。

batman equation

さっそくMathematicaのContourPlotで描こうとしたら、うまくいきませんでした。

仕方がないので、こんな感じにx方向に走査しながら方程式を満たす点を探して、一応描けました。

f[x_, y_] := ((x/7)^2 Sqrt[
      Abs[Abs[x] - 3]/(Abs[x] - 3)] + (y/3)^2 Sqrt[
      Abs[y + (3 Sqrt[33])/7]/(y + (3 Sqrt[33])/7)] -
    1) (Abs[x/2] - ((3 Sqrt[33] - 7)/112) x^2 - 3 +
    Sqrt[1 - (Abs[Abs[x] - 2] - 1)^2] -
    y) (9 Sqrt[
      Abs[(Abs[x] - 1) (Abs[x] - 3/4)]/((1 - Abs[x]) (Abs[x] -
           3/4))] - 8 Abs[x] -
    y) (3 Abs[x] + .75 Sqrt[
      Abs[(Abs[x] - 3/4) (Abs[x] - 1/2)]/((3/4 - Abs[x]) (Abs[x] -
           1/2))] -
    y) (9/4 Sqrt[Abs[(x - 1/2) (x + 1/2)]/((1/2 - x) (1/2 + x))] -
    y) ((6 Sqrt[10])/
     7 + (3/2 - Abs[x]/2) Sqrt[
      Abs[Abs[x] - 1]/(Abs[x] - 1)] - (6 Sqrt[10])/14 Sqrt[
      4 - (Abs[x] - 1)^2] - y)

points =
  Flatten[Table[{x, #} & /@ Cases[y /. NSolve[f[x, y] == 0, y], _Real],
    {x, -8, 8, 0.1}], 1];

ListPlot[points]

batman equation1

Mathematicaを持っている人は、pointsの計算の前にDistributeDefinition[f]を実行し、Tableの代わりにをParallelTableを使った方が速いでしょう。Mathematicaを持っていない人でもUMMで試せますが、ParallelTableは使えないようです。

こんなもんかなと思っていたら、Mathematicaできれいに描く方法がPlaying With Mathematicaで紹介されていました(The Batman curve)。パーツに分けてしまえば描けて当然という気もしますが、こっちのほうがきれいですね。

batman equation1

Universal Mathematica Manipulator—Poor Man’s Mathematica


Universal Mathematica Manipulator

Newer Version(新しいバージョン)

Wolfram CDF Player is a free of charge software to view documents in Computable Document Format (CDF). Users cannot edit codes using the software. For editing, Mathematica is required. Consequently, the relationship between CDF player and Mathematica is similar to the one between Adobe Reader and Acrobat.

Since CDF Player has almost all Mathematica functions to allow rich user experiences, users should take full advantage of those functions.

Although users cannot edit CDF documents with CDF player, they can still input data. If there were a document with a universal manipulator that receives program codes and outputs the evaluated results of them, then the document can perform as poor man’s Mathematica, or a complimentary alternative for Mathematica.

Users cannot use the ordinary Mathematica expressions as inputs, because, according to Wolfram research, CDF does not support non-numeric input fields [FAQ]. However, any expressions can be represented with numbers (e.g., Gödel number), and CDF supports numbers. Consequently, if there were a converter that converted Mathematica expressions to numbers and vice versa, CDF documents could process Mathematica expressions by using the converter.

Then, I developed such a converter. As a result, a CDF document that can process any Mathematica documents with the converter became reality. Therefore, I call the CDF document “Universal Mathematica Manipulator (UMM).”

Yet, one little problem remains. Users cannot paste on CDF documents unless he/she has Mathematica. As numbers that represent Mathematica expressions are very large, it takes a great deal of trouble to input them digit by digit. Luckily, both Mac OS X and Windows have scripting languages that can simulate key strokes: AppleScript for Mac OS X and VBScript for Windows.

The script in AppleScript is given below. It obtains a number from clipboard and emulates key strokes to input it digit by digit.

tell application "System Events"
     set str to get the clipboard
     keystroke "a" using command down
     keystroke str
end tell

This script is created on AppleScript Editor and saved in User Script Folder (~/Library/Scripts). Make sure that Script menu is shown in the menu bar in AppleScript Editor’s preferences. One downside to the script in AppleScript is to have to hit either “Enter” or “Tab” key to resume the script when it stops. Though, the way to eliminate the extra step is unknown.

Also, the script in VBScript is below. It is better than the one in AppleScript, as we do not need to press the enter key after the number is inputted. Nevertheless, this script has a limitation where it cannot handle an input with 5,000 or longer in length.

Set objIE = CreateObject("InternetExplorer.Application")
objIE.Navigate("about:blank")
text = objIE.document.parentwindow.clipboardData.GetData("text")
objIE.Quit

Set objShell = WScript.CreateObject("WScript.Shell")
WScript.Sleep(3000)
objShell.SendKeys("^a0{Enter}")
WScript.Sleep(1000)
objShell.SendKeys(text & "{Enter}")

In executing the script in VBScript, a warning message of Internet Explorer appears. You can disable the message at your own risk. (Tools > Internet Options > Security tab > Custom level… > Scripting > Allow Programmatic clipboard access > Enable)

This system is tested on Wolfram CDF Player 8.0.3.0 Mac OS X x86 and Microsoft Windows (32-bit). Relatively large codes like a solution for Burrau’s problem of three bodies can be executed.

solution for Burrau's problem of three bodies on UMM

Newer Version(新しいバージョン)

恩師の条件—魚偏の漢字


中学3年間、中勘助『銀の匙』1冊だけを教科書に行われたの国語の授業を紹介した『奇跡の教室』。この本に出てくる「魚偏の漢字は678字ある」という話を確かめられなかった、ということを以前書きました

4576050516予想はしていたのですが、「諸橋大漢和だと『恩師の条件』書いてある」と教えてもらったので、さっそくチェックしてみました。

恩師の条件自体は私にはあまり関係がないのですが、「私は高校生の勉強を“孤独な戦い”に終わらせたくなかった(p.170)」というのに、我が意を得たりなどと思いながら見ていくと、次のような記述がありました。

漢字を知らなければその魚が口にはいらぬわけでもないが、魚の漢字を見てどんな魚かわかるというのも日本人的教養というものであろう。

どんな字があるか、辞書を写せば簡単だがそれでは面白くない。趣味と教養を両立させたいものである。すし屋で出すマッチ、箸ぶくろ、茶呑みなどに、魚つくしを使っているのがある。ちょっと気をつけていると何種類も集まってくる。その読解をやってみるとなかなか面白い。大漢和辞典の六七八字の中にさえ含まれていないで、いくら考えても読めないものもあり、まるで難しいクイズを解くようなおもしろさが味わえる。「銀の匙研究ノート」

「魚偏」とは書いていませんが、諸橋大漢和であることは確認できました。重くてちょっとあれなのですが、引っ張り出して確認してみると、魚部の漢字が679字あるので、おそらくこれのことだったのだろうと思います(魚偏だけに限定すると625字)。補巻で追加された魚部が23字(うち魚偏は20字)などを考えるとまたよくわからなくなるので、「魚部は約700字」ということにするのがいいかと思います。

諸橋大漢和の情報ならUnicodeのデータベースに入っているじゃん」と思って、いつものようにUnihan.zipを漁(≠魚偏)ってみしたが、

grep kMorohashi Unihan_DictionaryIndices.txt | cut -f 1 | sort > morohashi.txt
grep "\\s195'\{0,1\}\." Unihan_RadicalStrokeCounts.txt | cut -f 1 | sort | uniq > 195.txt
join morohashi.txt 195.txt | wc -l

284

このように、284字しか出てきませんでした。CJK統合漢字拡張漢字に関しては、諸橋大漢和との対応はとられていないようですね。

魚​魛​魜​魝​魞​魟​魠​魡​魢​魣​魤​魥​魦​魧​魨​魩​魪​魫​魬​魭​魮​魯​魰​魱​魲​魳​魴​魵​魶​魷​魸​魺​魻​魼​魽​魾​魿​鮀​鮁​鮂​鮃​鮄​鮅​鮆​鮇​鮈​鮉​鮊​鮋​鮌​鮍​鮎​鮏​鮐​鮑​鮒​鮓​鮔​鮕​鮖​鮗​鮙​鮚​鮛​鮜​鮝​鮞​鮟​鮠​鮡​鮢​鮣​鮤​鮥​鮦​鮧​鮨​鮩​鮪​鮫​鮬​鮭​鮮​鮯​鮰​鮱​鮲​鮵​鮶​鮷​鮸​鮹​鮺​鮻​鮼​鮽​鮾​鮿​鯀​鯁​鯂​鯃​鯄​鯅​鯆​鯇​鯈​鯉​鯊​鯋​鯌​鯍​鯎​鯏​鯐​鯑​鯒​鯔​鯕​鯖​鯗​鯘​鯙​鯚​鯛​鯜​鯝​鯞​鯟​鯠​鯡​鯢​鯣​鯤​鯥​鯦​鯧​鯨​鯩​鯪​鯫​鯬​鯭​鯮​鯯​鯰​鯱​鯲​鯶​鯷​鯸​鯹​鯺​鯻​鯼​鯽​鯾​鯿​鰀​鰁​鰂​鰃​鰄​鰅​鰆​鰇​鰈​鰉​鰊​鰋​鰌​鰍​鰎​鰏​鰐​鰑​鰒​鰓​鰔​鰕​鰖​鰗​鰘​鰙​鰚​鰜​鰝​鰞​鰟​鰠​鰡​鰢​鰣​鰤​鰥​鰦​鰧​鰨​鰩​鰪​鰫​鰬​鰭​鰮​鰯​鰰​鰱​鰲​鰳​鰴​鰵​鰶​鰷​鰸​鰹​鰺​鰻​鰼​鰽​鰾​鰿​鱀​鱁​鱂​鱃​鱄​鱅​鱆​鱇​鱈​鱉​鱊​鱋​鱌​鱍​鱎​鱏​鱐​鱑​鱒​鱓​鱔​鱕​鱖​鱗​鱘​鱙​鱚​鱛​鱜​鱝​鱞​鱟​鱠​鱡​鱢​鱣​鱤​鱥​鱦​鱧​鱨​鱩​鱪​鱫​鱬​鱭​鱮​鱯​鱰​鱱​鱲​鱳​鱴​鱵​鱶​鱷​鱸​鱹​鱺​鱻​鷠​魯​鱗

4582128041ちなみに、一家に一冊『字通』には、魚部の漢字は48字あるようで、このくらいが平和でいいなあと思いました。さらにちなみに、.NET版『字通』はすてに古くなっていて、残念ながらWindows 7にはインストールできません。紙の『字通』はこの先もずーと読めるでしょうから、改めて「紙は強い」と思いました。とはいえ、諸橋大漢和は重くて場所取ってあれなので、早く電子化してほしいのですが、

CD-ROM化に関しては、コンピュータで扱える漢字の数の問題や、『大漢和辞典』そのものの巨大さがネックとなって、一朝一夕には実現できない状況です。(「大漢和辞典 よくある質問」のQ13

4469800007などと言うなら、とりあえず画像で、番号での検索のみ可能にしたバージョンを、DVDやiPad版で出せばいいんじゃないかと思います。ポスターならできるくせに。自炊している人とかいるんですかねえ。

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


描画が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;

計算機環境

  • CPU:Core i7 930(HTオフ)
  • 主記憶:6GB
  • GPU:NVIDIA GeForce GTX 465
  • OS:Windows 7 64 bit
  • Mathematica 8.0.1

結果

比較の基準

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

最初のバージョンの実行時間を計ります。

time0 = AbsoluteTiming[
   result = Table[ f[x + y I],
      {y, yMin, yMax, dy}, {x, xMin, xMax, dx}];][[1]]

23.5279876

CUDA

CUDAの準備をします(MathematicaからCUDAを使う時、書かなければいけないCのコードはこの部分だけなので簡単です)。

Needs["CUDALink`"]

kernelCode = "
  __device__ int check(Real_t x, Real_t y, int maxSteps) {
    Real_t re = 0, im = 0;
    int i = 0;
    while (i++ < maxSteps) {
      Real_t newRe = re * re - im * im + x;
      Real_t newIm = 2 * re * im + y;
      if (4 < newRe * newRe + newIm * newIm) return i;
      re = newRe;
      im = newIm;
    }
    return -1;
  }

  __global__ void mandelbrot(int* result, int maxSteps,
      int width, int height,
      Real_t xStart, Real_t xEnd, Real_t yStart, Real_t yEnd) {
    int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
    int yIndex = threadIdx.y + blockIdx.y * blockDim.y;
    if (xIndex < width && yIndex < height) {
      Real_t x = xStart + (xEnd - xStart) * xIndex / width;
      Real_t y = yStart + (yEnd - yStart) * yIndex / height;
      result[xIndex + yIndex * width] = check(x, y, maxSteps);
    }
  }
  ";

kernel = CUDAFunctionLoad[kernelCode, "mandelbrot",
   {{_Integer, _, "Output"}, _Integer, _Integer, _Integer,
    _Real, _Real, _Real, _Real}, {16, 16}];

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

実行してみると、4000倍くらい速くなることがわかります。

time = AbsoluteTiming[
    width = 300; height = 300;
    buffer = CUDAMemoryAllocate[Integer, {height, width}];
    kernel[buffer, 100, width, height, N@xMin, N@xMax, N@yMin, N@yMax];
    result = CUDAMemoryGet[buffer];][[1]];

CUDAMemoryUnload[buffer];

time0/time

3920.8

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

リアルタイム描画

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

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

width = 300; height = 300;
buffer = CUDAMemoryAllocate[Integer, {height, width}];

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

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がある場合はこれで十分なのかもしれませんが、私の環境では、「一応動くがストレスフル」という感じでした。

簡単な画像処理でライフゲームの初期データを作る方法


Mathematicaにはセルオートマトンのための汎用関数が用意されているので、簡単にセルオートマトンの一種であるライフゲームで遊べます。

しかし、初期状態を用意するのが面倒なので、ウェブの画像から初期状態を作ってみましょう。やり方は、QRコードの数値化と同じです。画像処理なんかしなくても、LifeWikiを探せば有名なものはだいたい見つかるでしょうが、こういうことは、別の機会に役立つかもしれません(Wolfram CDF Playerがインストールされていれば実際に動かすことができます。RSS Readerでは何も見えないかもしれません)。(UMMでも動きます。)

(* 画像を読み込む *)
gun = Import[
   "http://upload.wikimedia.org/wikipedia/commons/thumb/e/e0/Game_of_\
life_glider_gun.svg/610px-Game_of_life_glider_gun.svg.png"];

(* 白黒にし、数値化する。 *)
data = ImageData@Binarize@gun /. {0 -> 1, 1 -> 0};

(* 同じ値がいくつ続くかを調べる *)
unit = Take[Union[Flatten[Length /@ Split@# & /@ data]], 2];

(* サンプリングし直す *)
init = data[[Min@unit + 1 ;; Length@data ;; Total@unit,
    Min@unit + 1 ;; Length@data[[1]] ;; Total@unit]];

(* アニメーションにする *)
ListAnimate[ArrayPlot /@
  CellularAutomaton[{224, {2, {{2, 2, 2}, {2, 1, 2}, {2, 2, 2}}},
    {1, 1}}, {init, 0}, 150]]

この方法にも限度はあるので、人類が作り上げてきたものを堪能したいという人は、Gollyのような専用ソフトウェアがおすすめです。

Javaでサニタイズするときの注意—Apache Commons Lang 2.xと3.0の場合


Java 5をサポートする「Apache Commons Lang 3.0」がリリース、というニュースがありました。正確には、Java 1.4以下をサポートしない「Apache Commons Lang 3.0」がリリース、だと思いますが、それはまあいいとして、これに関連する少し深刻な話を忘れないうちに書いておきましょう。

「Javaプログラマは文字列の操作方法を確認した方がいいかもしれない」のつづきです。

HTML文書を書く時には、「𠮷野家(>_<)」のような文字列は、「𠮷野家(&gt;_&lt;)」のように書かなければなりません。「<」や「>」、「&」などは特殊な文字なので、文書中にそのまま書くことはできないのです。

効率を気にしないなら、次のように書いて置換するのが簡単です。

str = "\uD842\uDFB7\u91CE\u5BB6(>_<)";
str = str.replaceAll("&", "&amp;")
        .replaceAll("<", "&lt;").replaceAll(">", "&gt;");
System.out.println(str);
//𠮷野家(&gt;_&lt;)

しかし、このようなコードをいつも書くのは面倒ですし、順番を間違える危険もあります。実は「'」や「"」を置換したかったりもします。

4627847327そこで、HTML中に直接書けない文字を置換するためのライブラリを使います。Apache Commons Langはそのようなライブラリの一つで、拙著『Webアプリケーション構築入門』でも紹介しました。次のように使います。

まずは最近リリースされたCommons Lang 3.0の場合です。

str = "文字列(>_<)";
str = org.apache.commons3.lang.StringEscapeUtils.escapeXml(str);
System.out.println(str);
//文字列(&gt;_&lt;)

3.0は出たばかりですし、Java 5以上が必須なので、まだ2.xの方がほとんどでしょう。2.xでは次のようになります。

str = "文字列(>_<)";
str = org.apache.commons.lang.StringEscapeUtils.escapeXml(str);
System.out.println(str);
//&#25991;&#23383;&#21015;(&gt;_&lt;)

Commons Lang 2.xでは、ASCII以外の文字は「&#コードポイント(10進表記);」という形(数値文字参照)に置換されるので、人間にはわかりにくくなっていますが、ウェブブラウザでなら「文字列(>_<)」と表示されます。

2.xと3.0の間には、(1) パッケージ名(2.6はcommons、3.0はcommons3)(2) 数値文字参照の使用(2.6では使う、3.0では使わない)、に関して違いがあります。

ここまでは問題ありませんが、対象文字列がサロゲートペアを含むようになると、いやな感じになります。

まずは3.0の場合。

str = "\uD842\uDFB7\u91CE\u5BB6(>_<)";
str = org.apache.commons3.lang.StringEscapeUtils.escapeXml(str);
System.out.println(str);
//𠮷?野家(&gt;_&lt;

変なところに「?」が入って、最後の「)」が無くなっていますが、これは間違いです。数値文字参照を使わない3.0では、「𠮷野家(&gt;_&lt;)」にならなければなりません。

バグレポートを送ったので、この問題は3.0.1で解決されるはずです。

サロゲートペアを使う可能性が少しでもある人は、Commons Lang 3.0はスルーしましょう。

パッチを見てもらえばわかりますが、「1文字ずつ処理する」部分にバグがあったのです(参考:Javaプログラマは文字列の操作方法を確認した方がいいかもしれない)。

とりあえずの結論:Javaの文字列処理は難しすぎる。

次にCommons Lang 2.xの場合です。

str = "\uD842\uDFB7\u91CE\u5BB6(>_<)";
str = org.apache.commons.lang.StringEscapeUtils.escapeXml(str);
System.out.println(str);
//&#55362;&#57271;&#37326;&#23478;(&gt;_&lt;)

「𠮷」が「&#55362;&#57271;」という2つの数値文字参照になっていますが、これは間違いです。「&#134071;」という1つの数値文字参照にならなければなりません(参考:Using character escapes in markup and CSS)。

この問題を解決するパッチを送りましたが、「3.0で修正されている」としてスルーされてしまいました。間違いが2.xで修正されるとしたら、動作をそろえるためにStringEscapeUtils.unescapeXml(str)にも修正が必要だと思い、そのためのパッチも送ったのですが、こちらは現時点では間違った動作はしていないので、やはりスルーされてしまいました。

「3.0で修正されている」とは言っても、先に述べたように、2.xと3.0ではサポートするJavaのバージョンが違いますし、パッケージ名も変わっているので、2.xでも直しておいた方がいいと思ったのですが、まあ、仕方ありません。

Commons Lang 2.xを使っている人は、「サロゲートペアはサポートしない」と仕様書に明記しておきましょう。

最終結論:Javaの文字列処理は難しすぎる。