アルゴリズム（written by Taku Nakahara）
第８回は「アルゴリズム」をやりたいと思います。
wikipediaによるとアルゴリズムとは、
アルゴリズム (algorithm) とは、数学、コンピューティング、言語学、あるいは関連する分野において、問題を解くための手順を定式化した形で表現したものを言う。「算法」と訳されることもある。
「問題」はその「解」を持っているが、アルゴリズムは正しくその解を得るための具体的手順および根拠を与える。さらに多くの場合において効率性が重要となる。
コンピュータにアルゴリズムを指示するための（電子）文書をプログラムという。人間より速く大量に正しい結果を導くことができるのがコンピュータの強みであるが、そのためにはプログラムは正しく効率的なアルゴリズムに基づくことが必要である。
と定義されています。この定義にあるようにあらゆるプログラムは何らかのアルゴリズムを実装したものということができます。アルゴリズムの中でも特によく使われるものには、囲碁における｢定石｣にあたるような有名な指し手があります。たとえばめちゃくちゃに並んだ数字の配列を順番に並べる方法とか、乱数を用いて数値計算的を行う方法とか。。。
たとえば、我々のよく使うところで言うと、blastやclustalwなどのアラインメントを計算の一部もしくは全部で行うようなプログラムで使われているアルゴリズムに”dynamic programming (DP)”というのがあります。
Nature Biotechnology 22, 909 – 910 (2004)
doi : 10.1038/nbt0704-909
What is dynamic programming?
Sean R Eddy
The filled dynamic programming matrix for two DNA sequences, x = TTCATA and y = TGCTCGTA, for a scoring system of +5 for a match, -2 for a mismatch and -6 for each insertion or deletion.
The cells in the optimum path are shown in red. Arrowheads are ‘traceback pointers,’ indicating which of the three cases were optimal for reaching each cell. (Some cells can be reached by two or three different optimal paths of equal score: whenever two or more cases are equally optimal, dynamic programming implementations usually choose one case arbitrarily. In this example, though, the optimal path is unique.)
本当はもう一個、加算していないmatrixがあるととても説明しやすいのですが、まあ興味を持った人はこの文献をあたったり、配列解析の初歩が載っている参考書でも読んでみるといいでしょう。
今日は簡単なアルゴリズムをいくつか実装してみましょう。
モンテカルロ法で円周率πを求める
まずは乱数を用いて円周率を求めるアルゴリズムからやってみましょう。中にはたくさん覚えている人もいるかもしれませんが、
円周率 ≒ 3.14159 26535 89793 23846 26433 83279 50288 …
です。この値をコンピューターを用いてできるだけ正しく計算してみましょう。
乱数について軽く説明しておくと、これは「何の規則性もなくでたらめに発生する数｣と定義されます。適当に思いついた数をぱらぱら並べていけばそれが乱数のセットになる訳ですが、コンピューターもこれと似たようなことができます。ただし、コンピューターなのでいつかは前に言った数をもう一回言ってしまって同じことの繰り返しになります。できるだけ一様にばらけた乱数を発生させるアルゴリズムはいろいろあるようですが、完全な乱数をコンピューターで発生させることはできないので、コンピューターが作り出した乱数のことを「疑似乱数｣と呼ぶようです。
C言語では乱数を発生させるには、rand()という関数を使います。
#include int ransuu = rand();
という使い方をします。rand()をよびだすと0~32767までの整数を適当に返してくれます。
準備はこれくらいにして、下のような1/4の円を考えます。半径は１です。
このなかにx=(0~1), y=(0~1)の乱数を発生させて、それに対応する点をぽつぽつ打っていったとします。そうすると、全体の面積１で、円の中に入る確立はπ/4になるはずですね。全部の点の数をN、円の中に入った点の数をIとすると、
なので
ということになりますね。
まとめると、
この範囲に適当に多数の点(x,y)をぽつぽつ打つ。
円に入ったもの(x2+y2 < 1になった点）がいくつあるか数える。
４×（円に入った数）÷（打った点の数）を計算するとπがわかる。
という手順になります。実際にコーディングしてみましょう。
(pi.c)
#include <stdio.h>
#include <stdlib.h>
double getPoint();
int main(int argc, char *argv[])
{
int N = atoi(argv[1]);
int I = 0;
int i;
double x, y, pi;
double delta;
for (i = 0; i < N; i++) {
x = getPoint();
y = getPoint();
printf("(%f,%f)\t", x, y);
if ((x*x+y*y) < 1) {
printf("It's in!\n");
I++;
} else {
printf("Out!\n");
}
}
pi = (double)4*I/N;
delta = 3.141592653589793 - pi;
printf("%d/%d were in the circle area.\n", I, N);
printf("pi = %f\n", pi);
printf("d = %f\n", delta);
return 0;
}
double getPoint()
{
double rnd = (double)rand();
double point = rnd/RAND_MAX;
return point;
}
このプログラムを実行する時には引数としていくつ点を打つか(N)を指定します。この結果、πがいくつと計算されたか(pi = …)と、その値が既に分かっている円周率と比べてどれくらいずれているか(d = …)が表示されます。Nが大きくなるとpiがより3.1415…に近づいてdが小さくなるが分かるでしょう。
このように乱数をもちいて行う数値計算法のことを「モンテカルロ(MC)法」と呼びます。タンパク質のダイナミクスを計算するのにはよく分子動力学法(MD)が用いられますが、MCも分子計算ではよく用いられる手法です。分子計算の場合、ひとつ前の状態に依存して次の状態が決まるマルコフ鎖(MC)的な性質を持つのでMCMC(マルコフチェインモンテカルロ)法と言ったりします。この依存関係はメトロポリスの基準に従います。それが何かは自分で勉強してみてください。
バブルソート
ソート(sort)とは何らかの規則に従ってデータを並べ替えることです。たとえば{4, 7, 2, 9}という整数の配列があったとき、小さいものほど前にくるというルールでソート（要するに小さい順に並べる）と、{2, 4, 7, 9}になりますね、みたいのがソートです。これをプログラミング上で実現するアルゴリズムはたくさん提案されています。違いは計算量です。もっとも単純な直接選択法と呼ばれる方法は、n個の要素をもつ配列を並べる計算量はn2に比例します。それが改良された方法（たとえばクイックソートだとかヒープソートだとかの方法）だとnlog2nに比例するようになります。たとえば１００万個の要素(106)の要素を並び替えるとき、遅いアルゴリズム(O(n2))だと計算量が1012になりますが、速いアルゴリズム(O(nlog2n))だと2*107となり、遅いのと比べると５万倍もはやく並べ替えることができます。ただ、今日やってみるのは簡単な遅い方のアルゴリズムで、バブルソートと呼ばれます。下の図を見てください。
これは{80, 50, 56, 30, 51, 70}という整数の配列があった時にバブルソートが小さい順に並べ替える様子を模式的にあらわしたものです。パス１を見てください。この方法は下の方から上に向かって進み、自分の一個上の要素と比較して自分の方が小さいとき、自分の方が上に上がるという方法です。まず、70と51を比較して順番はOKなのでそのままです。次に51と30を比較してこれも小さい方が上なのでOKです。30と56を比較すると30の方が小さいのに下にいるのでここはひっくり返します。次に30と50を比較して、これも30の方が小さいのでひっくり返します。さらに30と80を比較し、これもひっくり返します。こういうことをやっていると、小さい値の要素がぷかぷか上に上がっていきます。だからバブルソートと言うのです。では実際にcodingしてみましょう。
(bubbleSort.c)
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char *argv[])
{
/*
* prepare integer array
*/
int i, j, tmp;
int N = argc-1;
int array[N];
for (i = 0; i < N; i++) {
array[i] = atoi(argv[i+1]);
}
/*
* bubble sort
*/
for (i = 0; i < N; i++) {
for (j = N-1; j > i; j--) {
if (array[j] < array[j-1]) {
/* swap */
tmp = array[j];
array[j] = array[j-1];
array[j-1] = tmp;
}
}
}
/*
* show results
*/
for (i = 0; i < N; i++) {
printf("%d ", array[i]);
}
printf("\n");
return 0;
}
結果は各自で確認してください。こういうアルゴリズムが自分で書けるようになると、例えばひとつの要素に複数の属性があるとき（タンパク質の質量と総電荷と糖鎖修飾数とか）があった時に、まず質量の順を優先して、質量が同じ時には総電荷で勝負を決めて、だけど糖鎖がある場合には、それはそれだけで別ランキングにするとか、そういう細かい設定でいろんなデータをソートできるんじゃないかと思います。
今日は三つの有名なアルゴリズムを紹介しましたが、実はこういった有名なアルゴリズムを自分自身で実装しなきゃならんということはないのです。世の中にはよく使うプログラムをライブラリやモジュールとして公開しているものもよくあるので、そういったものを利用すればよくある面倒なことはスキップできます。Bioinformatics系で言えば、bioperl, biojava, bioruby, biopythonなんかは有名です。ただ、こういった外部のライブラリに頼るにしても中でどんな計算が行われているのか全く分かっていないのはよくないです。今後もいろいろなプログラムやライブラリを利用して研究を行っていくことになりますが、ぜひ中でどんなアルゴリズムが使われているのか意識して、分からなければ調べるようにしましょう。