照度分布計算その5 - デバグを簡単に

どんどん行こう、そしてさっさと済まそう、回折計算。
具体的な実装を始めよう。まず、まっさきにデバグエイドから。

数値計算の場合、いかに簡単にデバグできるようにしておくか、というのは重要である。どんなソフトを書くときでも重要だけど、数値計算は特に

1. 第3者的な視点が得られにくい
2. 正しい結果とそうでない結果の区別がつきにくい
3. 中間結果が単なる数字の羅列である場合が多い

と言う特徴がある。

汎用のソフトなら憎たらしくもありがたいバグレポートを期待できるけど、数値計算の場合不特定多数ではなく限られた人、特に良くあるのは書いた本人しか使わないと言う場合で、冷静にソースを見るという視点が得られないことが多い。他人が見れば「なにこれ」というような明らかなバグでも気がつかないと言うことがある。

なんで計算したいか、と言えば誰も結果を知らないから、なわけで、従ってコンパイルができてちゃんと実行できたとしても、その結果が正しいかどうか判断するのは難しい。

また、たいていの場合、2次元や3次元のメッシュの上で値を求める、ということが数値計算の場合多く、結果は数字の羅列で正しいのやら間違ってるのやらなんにもわからん、64メッシュ程度の3次元データでも26万点以上の数字の羅列になる。これを人間がいちいちチェックするのは絶対不可能。

自分しか読まないコードであっても、読みやすさを心がける

ではどうするか、というと、まず客観的にソースを眺める、ということ。一晩寝てから見る、と言うのもいいけど一週間経ってから見る、というのではすっかり忘れていて自分で見てもピヨピヨ、ということになってしまうので一晩くらいがちょうどいい。それよりもさらにいいのは、やっぱり読みやすいソースを心がけること。変数の名前は具体的に付けて、対応する閉じカッコが20行を超えて下にあるときはまとめられるところを関数にまとめて、というごく普通に言われている書き方を心がける。

僕の場合、昔うまく行ったときと同じスタイルで書く、ということと、関数やメソッドの自動変数（auto変数）の数がなるべく少なくなるような書き方をするのが良い。

auto変数のかずを減らす

ある計算をする場合にパラメータが10個あるとする。10個のパラメータを含んだひとつの関数（あるいはメソッド）では大きすぎる。普通パラメータ間には構造があって、10個のパラメータがフラットに（同時に）必要になることはない。

その構造に従って関数も構造化することで、関数の中でさらに別の関数を呼ぶことにすればそれぞれの関数で必要なパラメータの数を減らすことができる。まず一番外の関数でパラメータ4つ、それがまた別の関数を呼んでそれが3つのパラメータを持ち、さらにそれがまた別の関数を呼んでその関数が残ったパラメータを持つ、などというようにする。

まあ、こんなふうにいつもうまく行くとは限らないけど、ひとつの関数の内部で使うパラメータが2、3個でループ変数がせいぜい2、3個というのが望ましい。そんな風に書けないときはたいてい考えている問題が複雑なのではなくて、考えが整理されていないと言うことの方が多い。関数の中の変数が減らせるということは、数学が整理されている、ということだと思う。これはコーディングのスタイルとか言う以前のアルゴリズムがシンプルであるということで、関数の中の変数の数というのはその現れだと思う。

数学を整理してドキュメントする

今回、すごくはしょって書いてしまったけど、コーディングする前にまず先に数学を整理して、変数の名前等を決めてしまう、ということも有益。コーディングのときはその名前をそのまま使って書く。数学はちゃんとTeXなんかで書いてまとめておけば一週間経った後でもピヨピヨ状態の時間を減らすことができる。

トリビアルな結果も計算できるようにしておく

だれも知らない結果の正しさをどうやって判断するか、というと結果が分かっている場合の計算もできるようにしておくというのは重要。ごく特殊な条件では解析的に解が求まる、ということがたまにある。これは非常に大切で、その解析的な結果と同じ条件での数値計算の結果を較べられるようにしておく。よくあるのが、そういうトリビアルな場合の計算はするんだけど、解析的な解のほうを自分で導けない、あるいは数値計算と比較できるような形で書いておくことができない、と言う場合。教科書を引張り出して図を較べて「似てる」なんて言うことがある。これは計算以前の数学をもっと勉強しましょう、ということにつきる。数値計算万能主義者に多い。

デバグ用にグラフ表示の手段を持つ

膨大な数字の羅列を人間が解釈できるようにするには、もう視覚に訴えるしかない。とりあえず中間結果を含めてグラフ表示できるようにしておく。例えば僕は大きな数字の配列を持つクラスにはこんなプロトコルを強制している。

@protocol SMathematicaStyleDescription
- (NSString *)mathematicaStyleDescription;
@end

これは中身をmathematicaで読めるようなリストの形で文字列を書き出すと言うもの。デバグが必要になったら

    NSString    *debstr = [obj mathematicaStyleDescription];
    [debstr writeToFile:@"../../out.txt"
             atomically:YES
               encoding:NSASCIIStringEncoding
                  error:&error];

などと書いてMathematicaで

In[1]:= result=<<"result.txt";

みたいにして読めば、あとはMathematicaの上で表示できる。 例えば2次元の複素数の配列を持っているクラスの場合

- (NSString *)mathematicaStyleDescription
{
    NSMutableString *ret = [NSMutableString stringWithFormat:@"{"];
    int     h, v;
    complex *cp = cPlane;

    for (v = 0 ; v < sSize ; v ++) {
        [ret appendString:@"{"];
        for (h = 0 ; h < sSize ; h ++) {
            [ret appendFormat:@"%10.8f %+10.8fI, ",  // (1)
                        cp->re, cp->im];
            cp ++;
        }
        [ret appendString:@"Sequence[]},\n"];  // (2)
    }
    [ret appendString:@"Sequence[]}\n"];  // (3)
    return ret;
}

などと実装しておく。ここでcPlaneというのはこのクラスのインスタンス変数で2次元のパックされた複素数の配列へのポインタ。

ここでtipsをひとつ。Mathematicaで

{a,b,c,Sequence[d]}

は、Mathematicaに読み込んで式として解釈されたら

{a,b,c,d}

と言う意味になる。Sequence[]はListの中ではflatに繋ぎ合わされる。MathematicaのListに書き出そうとしたとき、要素同士を分けるコンマ","は要素の数よりちょうどひとつ少ない必要があるけど、最後にSequence[]（引数が空）を入れれば要素の数を知らなくてもコンマを要素と同じ数だけ書いてリストを閉じることができる。Mathematica式を文字列で書くとき特有のコツ。

上のコードでは(1)のループの中で要素を書き出して、その後ろにコンマ","を書いているが、最後の要素の後ろにもコンマがつく。(2)と(3)にSequence[]を付け加えることで最後のコンマを無かったことにできる。ループの中で最後の要素かどうかを判断する必要が無くなるので簡単になる。

僕はMathematicaを使える環境を持っていて慣れ親しんでいるのでこれが便利だけど、人によってはMATLABやEXCELが便利な人もいるだろう。その人は数字をタブ区切りやコンマ区切りで書き出せばいい。コンマ区切りのときはSequence[]なんで便利な物は無いのでコンマの数を合わせないといけないけど。

とにかく数値計算のソフトを書くときは、どんなに小さくてもこういうデバグのための道具は持っていた方がいい。今回も使う。

具体的な実装に入れなかったな。でも今回はの内容は重要。