An intuitive explanation for the "Russian Roulette" sampling technique in ray tracing

The "Russian Roulette" sampling technique (which is creepily named, but is known officially as so) is a sampling technique used in ray tracing to reduce the number of samples for evaluating some value, while maintaining the convergence. It was very mystifying at first sight, but became very intuitive when I came up of a graphical explanation.

Shortly put, the Russian Roulette sampling technique can be interpreted as making "copies" of high-order terms, to use in place of the dropped off terms that were not evaluated in other samples. This can be described graphically using the following figure.

f:id:woodrush:20150725201736p:plain

Let us consider a problem where we want to evaluate 4 orders of some value, with 5 samples. In ray tracing, this could be evaluating the radiance using terms of up to 4 reflections with 5 samples, for instance. In the figure above, the colors correspond to the orders of the value, and the number n in the circle indicates that it is the n-th sample.

The leftmost figure shows a direct method of evaluating this value, that is to just simply evaluate all of the orders in all of the samples. However, in ray tracing it can be time-consuming to do such evaluations.

This is where the Russian Roulette sampling technique comes in. In the Russian Roulette technique, we aim to evaluate only a few high-order terms. Basically, before evaluating each order of terms, we roll some dice - and proceed to evaluate high order terms at only a given probability, or else halt the evaluation and move on to the next samples. Each higher order term must survive this judgement before it gets evaluated, or otherwise the "game" (higher-term evaluation) ends and moves on to the next round (sample) - therefore, the name "Russian Roulette." In this figure, we've assumed a 50% probability of proceeding to evaluate the next higher order. The figure on the center shows the actually calculated terms in each sample, for the Russian Roulette. In other words, the circles in the figure are the circles that "survived" the Roulette. Notice that the number of evaluations is reduced to around half of that in the left figure, and many high orders are dropped off.

Now, since we have excluded many orders, we don't have enough circles (terms) to add up to calculate the average. (i.e., simply taking the average of all of the terms (circles) provides an inaccurate estimate of the true value, compared to the estimator shown on the left figure.) Specifically, each n-th order term shows up at a ratio of (1/2)^(n-1) compared to the number of them in the left figure.

How do we calculate the average with such few circles (terms)? The answer, which is the Russian Roulette sampling technique, is to use copies of other circles (terms) as placeholders. Mathematically, the technique is to divide the value of each evaluated term by the probability of its appearance when taking the grand average. This obscure technique can be understood by observing the right side of the figure. Dividing the term by its survival probability p is equivalent to considering as if 1/p copies of that term has existed in the evaluation. This is the whole idea of the sampling technique, and is depicted in the figure above as copies of the circle. (Copies are indicated by red index numbers.) By making "copies" in such a way, the number of samples of each order approximately match that of the direct technique shown in the left side of the first figure. Although the number of circles differ from its true value (they are short in the 3rd order and excessive in the 2nd and 4th order), if we take a sufficient number of samples, the law of large numbers helps us to diminish those discrepancies to an ignorable level. Therefore, the "magic" of dividing higher order terms in the Russian Roulette technique, can be interpreted as making the "right number of copies" of that term, so that it properly compensates for the higher order terms that did not survive to be evaluated (i.e. dropped in the sampling process).

MATLABで、for文・parfor文中で関数内のパラメータを様々に変えて計算をする方法(クロージャの利用)

たとえば以下のような状況があると思います。

  • ある微分方程式の解軌道を、for文中で異なるパラメータごとに描きたい
  • ある方程式の解を、for文中で異なるパラメータごとに求めたい

このようなとき、ode45やfsolveを使いますが、for文中で関数の中身のパラメータの値を直接いじることはできません。

あまりよくない方法(parfor内では使えない方法)

そのため、たとえばglobalを使ってグローバル変数を定義し、関数内から変えたりする方法を取ったりします。例として、関数odefunc内のpという変数を様々に変えてode45をする場合、次のようになります:

% odefunc.m
function ret = odefunc(t,x)
    global p;
    ret = [-p*x(2); x(1)];
end
% main.m
global p;
for p = 1:10
    [t, y] = ode45(@odefunc, linspace(0,10,100),[1; 1]); 
    figure;
    plot(t,y);
end

この方法は比較的簡単ですが、parfor内では使えないという重大な欠点があります。parforはfor文と同じ文法で簡単に並列計算・高速化を行うための方法(MATLABのparforの使い方メモ(高速化のために) - woodrush’s diary 参照)なのですが、globalな変数を用いていると並列計算が行えないため、この方法が使えません。

良い方法 - 「関数を返す関数」を使う方法(parfor内で使える方法)

これは、「関数を返す関数」を使うとかなり綺麗に解決することができます。まずはソースを見てみましょう:

% odefuncFactory.m
function retfunc = odefuncFactory(p)
    function ret = odefunc(t,x)
        ret = [-p*x(2); x(1)];
    end
    retfunc = @odefunc;
end
% main.m
parfor p = 1:10
    odefunc = odefuncFactory(p);
    [t, y] = ode45(odefunc, linspace(0,10,100), [1; 1]); 
    figure;
    plot(t,y);
end

ここでodefuncFactoryという怪しい関数が出てきました。これはその名の通り、「パラメータpを与えると、そのpを使ったodefuncを返す関数」つまり「関数を返す関数」になっています。main.mの3行目にて、odefuncFactoryに所望のpを入れ、そのpを使ったodefuncを受け取っています。「odefunc = 」と書くことで、返ってきた関数をodefuncという変数に格納しています。(このように、MATLABでは関数も変数のように扱うことができるのです。)受け取ったodefuncをode45に渡し、シミュレーションを行っています。

この方法を使うと、異なるpごとに異なる関数を毎回作成することになるので、parfor内でも使えるようになるのです。このことを強調するために、2つ目のmain.mではfor文をparfor文にしました。

おまけ - 更に簡潔に

odefuncの程度のような短い関数ならば、無名関数を使うと更に短く書くことができます:

% main.m
odefuncFactory = (@ (p) (@ (t,x) [-p*x(2); x(1)]));
for p = 1:10
    odefunc = odefuncFactory(p);
    [t, y] = ode45(odefunc, linspace(0,10,100), [1; 1]); 
    figure;
    plot(t,y);
end

% 更に短く
parfor p = 1:10
    odefunc = (@ (t,x) [-p*x(2); x(1)]);
    [t, y] = ode45(odefunc, linspace(0,10,100), [1; 1]); 
    figure;
    plot(t,y);
end

main.mだけになってしまいました。一つ目の例のように、無名関数も「関数を返す関数」の作成に使えます。更に、odefuncくらいの関数なら、odefuncFactoryを作るまでもなく、二つ目の例のように直接odefuncを無名関数で構成することもできます。
ただし、最初の例では無名関数odefuncFactoryをparfor内で定義しているという違いに注意して下さい。無名関数を使うと何かとparforの制約に引っかかりやすくなるので、通常は個別のmファイルを作ることをおすすめします。

もう少し詳しく

この一連方法では「関数を返す関数」よりも、クロージャという概念が本質的に活用されています。2つ目のmain.mを見てみましょう。main.m中には「odefuncFactory関数の中にodefunc関数が包まれている」という階層構造があることに注目して下さい。次に、odefuncFactory.mのソースをよく見ると、「function ret = odefunc(t,x)」の中には変数pの定義が含まれていない(p= ... と書いていない)ことに気づきます。普通ならばどこにもpの定義が書いていなければエラーが起こります。それでも動くのは、odefuncの外のodefuncFactoryでpが定義されているからです。odefuncは、その外のodefuncFactoryに存在する変数pの値が「覗ける」ので、このプログラムは動くわけです。つまり、odefuncFactoryは、odefuncと共に、odefuncの外側の変数たちの知識、つまり外側の環境へのアクセスも与えています。このように、関数とその外側の環境のペアをクロージャと言います。odefuncFactoryは「関数を返す関数」と言いましたが、更に正確に言うと「クロージャを返す関数」だったわけです。

このようにMATLABの無名関数はクロージャを返してくれるという性質から、色々とマニアックなことができます。たとえば、条件分岐や複数行実行などが行えたりします。詳しくはMATLABの無名関数内で条件分岐・複数行実行を行う方法 - woodrush’s diaryを御覧ください。

MATLABのparforの使い方メモ(高速化のために)

MATLABのparfor構文を使うと、parfor内の処理を並列で行ってくれるため簡単に高速化ができます。たとえば多くのパラメータに対するシミュレーション結果を得たいとき、parforを使うと各パラメータに関する計算を簡単に並列化することができます。

しかし、parfor内の処理には、並列化のためのやっかいな制限があります。中でも「結果を配列に書き込むときは、異なる並列計算の書き込み先インデックスが競合しないことが保証されなければいけない」という制限がくせもので、これ自体は当たり前な事なのですが、これをMATLABに納得させるようにプログラムを書くのが少し面倒だったりします。
以下に、parforがちゃんと動くコードを示します。内容は、ある2パラメータtheta1とtheta2の直積の範囲で計算を走査するプログラムとなっています。

% 走査したいパラメータ
theta1 = 1:100;
theta2 = 1:100;

% parfor内から書き込みたい配列のサイズを、parforのインデックスの最大値として明示する
resultMatrix = zeros(length(theta1),length(theta2));
theta1indmax = size(resultMatrix,1);
theta2indmax = size(resultMatrix,2);
% theta1ind, theta2indを陽にインデックスに指定する・その1
parfor theta1ind = 1:theta1indmax
    for theta2ind = 1:theta2indmax
        % theta1ind, theta2indを陽にインデックスに指定する・その2
        resultMatrix(theta1ind,theta2ind) = theta1(theta1ind)^2 + theta2(theta2ind)^2;
    end
end

contourf(resultMatrix);

% parfor中で配列に書き出さない時は、上記の制限を気にせず書くことができる
parfor theta1ind = 1:length(theta1)
    for theta2ind = 1:length(theta2)
        [theta1(theta1ind) theta2(theta2ind)]
    end
end

画面を割ってしまいタッチ操作不能になったAndroidからadbでデータを救出・移行した記録

先日、Xperia Z3をアスファルトに落としてしまい、画面が操作不能になってしまいました。僕は画面ロックにパターンロックを使っていたのですが、これの解除に一苦労しました。結果的にadbとその仲間たち、特にパターンロックの解除にuiautomatorが大活躍しました。
以下に、パターンロックを開く→データを救出・移行→初期化までの流れをメモしておきます。

この方法の条件・注意点としては、

  • USBデバッグが必要。前提条件として「開発者向けオプション」を解禁しており、かつUSBデバッグモードをONにしている必要があります。
  • JDKAndroid SDKEclipse、およびantなど多くのツールを使います。
  • root化は必要ありません。

また、これは私の携帯の復旧記録なので、実機で行う場合は自己責任でお願いします。

1. uiautomatorを使ってパターンロックを解除する

adbではキー入力や一方向のスワイプ操作はできるのですが、連続した折れ線スワイプはできないため、adbだけではパターンロックが解除できません。ここを解決するのが最も難しい箇所でした。これにはuiautomatorを使って行いました。uiautomatorはJUnitを用いたAndroid向けのUIテスト用の機能で、複雑な操作を行うことができます。
ちなみに未確認ですが、スワイプパターンの入力だけならAndroid用のマウスなるものでもできるそうです。

uiautomatorviewerでスワイプ座標を取得

まずは、スワイプパターンの座標を取得するためuiautomatorviewerを使いました。
これには以下のサイトを参考にしました:
スマホ向け無料システムテスト自動化ツール(3):Androidテストで便利なuiautomatorviewer、UiScrollableの使い方、テキスト入力API制限事項の回避方法 (1/4) - @IT

uiautomatorクラスの作製

続いて、uiautomatorを使って実際にスワイプパターンを入力していきます。
これには、以下のサイトを参考にしました:
[Android] はじめてのuiautomator - adakoda

このサイト中の、次の4つの手順を行いました:

  • テストプロジェクトを作成する
  • テストクラスを作成する
  • ant を使用してビルドする
  • ビルドしたファイルを使用してテストを実行する

実際に作製したテストクラスは以下のものでした。(座標の値は変えてあります。)

// XperiaRedemption.java
package com.woodrush.xperiaredemption;
import android.graphics.Point;
import com.android.uiautomator.testrunner.UiAutomatorTestCase;

public class XperiaRedemptionTest extends UiAutomatorTestCase {
    public void testUnlockScreen() throws Exception {
        Point[] points = new Point[4];
        points[0] = new Point(100,200);
        points[1] = new Point(300,400);
        points[2] = new Point(500,600);
        points[3] = new Point(700,800);
        getUiDevice().swipe(points, 10);
    }
}

テストクラスの作製には以下のサイトを参考にしました:
Android - UiAutomatorを使用しての多段swipeの実現 - Qiita

次にantを使ってこれをビルドするため、ターミナルから以下を実行しました。

android create uitest-project -n XperiaRedemption -t 15 p .
ant build
adb push ./bin/XperiaRedemption.jar /data/local/tmp
adb shell uiautomator runtest XperiaRedemption.jar -c android.woodrush.xperiaredemption.XperiaRedemptionTest

これらの三行目までを行えば、あとは任意のディレクトリから4行目を実行するだけでロックパターンが入力できるようになります。

2. adbを使ってタッチ・キー操作

パターン入力後の操作は全てadbで行いました。これには以下のサイトを参考にしました:
Remote Input shell scripts for your Android Device, or my screen is cracked | The Grymoire

usage: input ...
 input text <string>
 input keyevent <key code number or name>
 input [touchscreen|touchpad] tap <x> <y>
 input [touchscreen|touchpad] swipe <x1> <y1> <x2> <y2>

これらのコマンドを使うには、まずターミナルからadb shellを実行してAndroid側のシェルを起動し、そこでコマンドを入力していきます。
特に、戻るキーが押せなくなっていたのですが、これは

input keyevent 4

を実行して押しました。
キーコードはKeyEvent | Android Developersを参照しました。

3. adb backupを使ってデータを救出

バックアップは、アプリデータのバックアップと、写真などのバックアップを行いました。
写真については、Xperia付属のソフトで行いました。
後から同ソフトでアプリのバックアップも行えることに気付き、このソフトも使ったのですが、最初にadb backupも使ったのでこちらについても記しておきます。

adb backupの引数は以下のようになっています:

  adb backup [-f <file>] [-apk|-noapk] [-shared|-noshared] [-all] 
              [-system|-nosystem] [<packages...>]
            - write an archive of the device's data to <file>.
              If no -f option is supplied then the data is written
              to "backup.ab" in the current directory.
              (-apk|-noapk enable/disable backup of the .apks 
                 themselves in the archive; the default is noapk.)
              (-shared|-noshared enable/disable backup of the device's
                 shared storage / SD card contents; the default is 
                 noshared.)
              (-all means to back up all installed applications)
              (-system|-nosystem toggles whether -all automatically 
                 includes system applications; the default is to 
                 include system apps)
              (<packages...> is the list of applications to be backed 
                 up. If the -all or -shared flags are passed, then the
                 package list is optional. Applications explicitly 
                 given on the command line will be included even if 
                 -nosystem would ordinarily cause them to be omitted.)

今回は新しい携帯の本体がすでにあるため-nosystemとしました。

adb backup -f xperiaredemption.ab -apk -shared -all -nosystem

また、念のため新しい携帯の完全バックアップをとるべく、-systemフラグを入れた状態で以下も実行しました:

adb backup -f newandroid.ab -apk -shared -all -system


バックアップ後はデータの移行を行います。
これにはadb restoreを使いました:

adb restore xperiaredemption.ab

このコマンドでは、バックアップを行ったアプリのみが、バックアップのデータをもとに上書きされる、という具合に移行が行われました。

また、携帯の初期化はadbで画面操作を行いながら携帯の機能を使って行いました。

おしまい

今回はUSBデバッグモードをONにしておいて助かりました。しかしやはり携帯を壊してしまったのは前の携帯にも申し訳ないです。僕はカバーはあまり好きではないのですが今回の一件でカバーを買ってあげるのが一番だと思いました。

MATLABで自動的にfigureを重ならないように並べる

MATLABで、figureが互いに重ならないように並べるgridfigureというコマンドを作りました。
gridfigure - File Exchange - MATLAB Central

MATLABのfigureには'defaultfigurecreatefcn'という隠されたプロパティがあり、これを変えることでfigureウインドウ作成時に任意の関数を実行することができます。
関数はそのfigureが初めて作成された時に一度だけ呼ばれ、以降figure(1)などでフォーカスを当てても呼ばれません。
他のfigureウィンドウを作ったり、figureを一度閉じた後にもう一度開いた時も呼ばれます。

f:id:woodrush:20141212210244p:plain

この機能を使い、ここにfigureをグリッド状に並べる関数を登録しておくことで、何もコマンドを入力しなくてもfigureが勝手にグリッド状に並ぶようになります。
複数のfigureを出して見比べたい時にかなり重宝すると思います。

更に、手動でfigureを並べ替えられるgridfigureというコマンドを作り、

  • マルチディスプレイ対応(サブディスプレイに表示する順番を選べる→メインディスプレイで作業しながら、サブディスプレイにfigureを描画可能)
  • ウィンドウサイズを指定可能
  • 並べ方を指定可能(左から右、上から下)

という機能を付けました。
更に、一度gridfigureでウィンドウサイズなどを指定しておくと、以降出てくるfigureのサイズがすべて指定のサイズになります。

過去にも手動でfigureを並べるコマンドはあったのですが、defaultfigurecreatefcnを使って自動でfigureを並べるパッケージは無かったので、その点はかなり便利になっていると思います。

Generalized "Make 10 with 4 numbers game" solver

「4つの数字で10を作るゲーム」を一般化した、「与えられたn個の数字の四則演算(追記:および括弧)でNを作るゲーム」の解を一つ求めるコードをPythonで書きました。
ポイントは、Generator機能を再帰して使っていることです。これにより、計算式の生成がかなり簡潔に書け、更にn=100でもあまりメモリを使わずに動くようになっています。

(11/27追記:コードを簡潔に書き直しました。)

provideNumsFromUserFlag, numOfNums, randmaxをいじると色々モードが変えられます。「1から100までの(ランダムな)100個の数字で10を作るゲーム」の解の出力は

[1, 2, 3, 4, 5, 9, 10, 11, 12, 12, 13, 13, 14, 17, 19, 20, 20, 20, 20, 21, 21, 21, 21, 22, 23, 23, 23, 24, 25, 26, 27, 28, 29, 30, 32, 35, 37, 37, 41, 41, 42, 43, 46, 47, 49, 50, 50, 50, 50, 50, 50, 51, 54, 54, 55, 55, 58, 59, 60, 62, 63, 64, 66, 67, 67, 67, 69, 69, 71, 72, 72, 73, 73, 75, 75, 76, 78, 80, 81, 82, 82, 82, 84, 86, 86, 86, 88, 89, 92, 94, 95, 96, 98, 98, 98, 98, 99, 99, 100, 100]
Calculating...
10 = (((((1-2)-((((3-4)-5)-9)-((10-11)-(((12-12)-(13-(13-14)))-(17-19)))))-(((20-20)-((((20-20)-(((21-(21-(21-21)))-(22-((23-23)-(23-24))))-((25-((26-27)-(28-29)))-30)))-((((32-(35-37))-(37-(41-41)))-(42-(43-46)))-(((((47-(49-50))-((50-50)-50))-(50-50))-(51-54))-54)))-((((55-55)-58)-(((59-((60-62)-63))-(64-(66-67)))-(67-67)))-(69-69))))-((71-72)-(72-73))))-73)-(((((75-75)-(((76-78)-80)-81))-82)-(((((82-82)-(84-86))-((86-86)-(((88-89)-92)-(94-95))))-(((96-98)-98)-(((98-98)/99)-99)))-100))-100))

のようになります。
かなり時間がかかると思っていたのですが、なんと数秒で結果が出たので驚きました。

関数allEquationsFromが式を(ある程度)重複無く生成する心臓部で、PythonのGeneratorを使っています。Generatorの再帰については、
再帰とジェネレータ
yieldのちょっとした理解
を参考にしました。
また、ここで使っている二分木の再帰的な構成方法については、
Enumeration of Binary Trees
をかなり参考にしました。

実行時間について

何度か試したところ、n=100でも5秒程度で解が求まる場合がほとんどでした。n=50では1秒未満で求まりました。
nが大きいと、式中の演算をほぼ全てマイナスに選び、ランダムに括弧付けを行うだけのアルゴリズムでもすぐに10が作れました。
数を多く使うと10は作りやすいそうです。

しかし、ほとんどの数が1であるような場合などは、マイナスだらけになると10にならないため、探索時に+を優先的に使って式を作らないとなかなか10になる式にたどり着かないので、数の組み合わせによっては工夫が必要です。(1が10個ある場合など。)

少し気になったので詳しく調べてみると(なぜやってしまったのか…)、数の数n=50, 目的の数N=50の場合について、演算子をマイナスだけに限定し、式の値を69937サンプル取った時の分布が下図になります。
f:id:woodrush:20141125061930p:plain
平均がほぼゼロ(-9.5949)の、正規分布に近いような分布になっています。引き算だけで式を作っているので、1から50までの50個の数を選んだとき、-25から25までの25個の数の和を得ているような感じになり、分布が正規分布に近づいているのだと思います。

ちなみに最初は木を作りやすそう&直接式を評価できる&Lisp練習しようという理由でLispで書いていたのですが、loop、let、progn、returnなどを多用していてあまりそれっぽくなりませんでした。
Lisp版
n=4ならば十分な時間で全数探索できるプログラムになっています。

UnityのTransformの子要素機能は大量に要素を詰めると重い

Sceneロード実行時に128^2個(16384個)のGameObjectをInstantiateし、それを全てUnityのGameObjectに直接の子Transformとして指定するようにしたところ、Sceneのロードがかなり重くなりました。

更に、各子要素のあるメソッドを呼んでいたのですが、これはInstantiateするのと同じくらい時間がかかっていました。

回避策として子要素に追加する処理を除いたところ、動作が従来の3倍も速くなりました。
Hierarchyは見難くなるものの、Unityでは大量のGameObjectをある要素の子にするのはやめた方が良いようです。

// 遅いコード
GameObject p = ...;
for(int i = 0; i < 128*128; i++){
    ((GameObject) Instantiate (...)).transform.parent = p.transform;
}