Fortranにやってもらう

Rをある程度習熟すれば、考えたことをなんでも実行できる。塩基配列の取扱いはBiostringsパッケージが強力だし、parallelパッケージの並列処理で多くの作業を高速化できる。

だけれども、やりたい計算があまりにもマニアックでapply系が使えないとき、とても困る。for()を多用せざるをえないような計算だ。

一応、for()も高速化はできる。先に適切なオブジェクトをつくっておき、その中のある座標だけを入れ替えるようにfor()ループを組み立てておけば、それなりに速く計算が終わる。宣言する感じ。そのつどオブジェクトを作り直すように組み立ててしまうと極端に遅い。それに、巨大なデータを用いて計算すると、速いといっても大したことない。

そこでFortranの出番となる。

RとFortranをつなぐのは.Fortran()という関数。Fortranのコードはプログラムではなくサブルーチンとして用意する。RからFortranに変数を渡せるし、FortranからRに変数を返してもらえる。ただし、サブルーチンの記述通りの順番で変数を渡す必要があるし、渡した変数しか返してもらえない。

サブルーチンの最後まで行った時点で保持していた変数のうち、Rから受け取ったものと同じもののみが、Rに帰ってくる。リストとして束ねられている。Fortranサイドで得た計算結果を受け取るためには、あらかじめRから初期値を入れた変数を与えて、それをFortranサイドで上書きしてしまえば良い。

RからFortranに変数が渡されると、Fortranが計算する。Fortranから変数が戻ってくるとき、Rのオブジェクトの中身だけを入れ替える形で戻ってくるらしい。だから、Fortranサイドでの宣言の仕方によってはおかしなことになる。

たとえば character(3) x みたいな感じで３個しか文字を格納できないようにしていると、Rから「x = "ABCDE"」というふうに５文字を受け取ったとしても Fortranで使うのは最初の３文字だけとなる（"ABC"）。この３文字を仮に別の文字に入れ替えてみるとする。たとえば、「x(3:3) = 'a'」とすると"ABa"になるのだが、このxをRに返してやると、上書きで返すのでRが受け取る変数は最初の３文字だけが改変された５文字、つまり"ABaDE"ということになる。

数字の並びを部分的に置換したいとき

数字の並びを部分的に置換することを考える。

たとえば、あるinteger vectorの中に存在する「3, 4, 5」を「5, 4, 3」に置き換えたいとき。
ここでは、文字列用の関数であるpaste()とgsub()を使ったトリックを紹介する。

適当なinteger vectorを作る。
int <- 0:9
int.vec <- c(int, int, int)

> int.vec
 [1] 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9


【for()ループでやる場合】

int.vec.new <- int.vec
for(i in 1:length(int.vec)){
if(int.vec[i] == 3 & int.vec[i+1] == 4 & int.vec[i+2] == 5){
int.vec.new[i:(i+2)] <- c(5, 4, 3)
i <- i + 2
}

> int.vec.new
 [1] 0 1 2 5 4 3 6 7 8 9 0 1 2 5 4 3 6 7 8 9 0 1 2 5 4 3 6 7 8 9


int.vecの中の「3, 4, 5」を「5, 4, 3」に置き換える仕事をfor()ループでやろうとすると、
このように一応はできるのだけれども、for()ループは遅い。
とくにvectorが長いほど遅さを実感する。

そこで、paste()とgsub()を使ったトリックの登場。


【paste()とgsub()でやる場合】

chr.vec <- paste(int.vec, sep = "", collapse = "")
int.vec.new <- as.integer(strsplit(gsub("345", "543", chr.vec), split = "")[[1]])

> int.vec.new
 [1] 0 1 2 5 4 3 6 7 8 9 0 1 2 5 4 3 6 7 8 9 0 1 2 5 4 3 6 7 8 9


まず、int.vecの値を文字列として繋げてしまう。
paste() はintegerを与えられたとしても、それをcharacterとして扱う。
collapseをうまく使えばすべての文字が繋がって、値がひとつのcharacter vectorになる。

chr.vec <- paste(int.vec, sep = "", collapse = "")

> chr.vec
[1] "012345678901234567890123456789"

こうしてひとつなぎの文字列にしてしまえば、gsub()が使える。
sub()は１つ目のヒットだけを置換するが、gsub()はすべてのヒットを置換する。

gsub("345", "543", chr.vec)

> gsub("345", "543", chr.vec)
[1] "012543678901254367890125436789"

あとは、strsplit()で１文字ずつ切り出し、integerに戻す。
strsplit() はリストで返すので、その１つ目を採用する。

> strsplit(gsub("345", "543", chr.vec), split = "")
[[1]]
 [1] "0" "1" "2" "5" "4" "3" "6" "7" "8" "9" "0" "1" "2" "5" "4" "3" "6" "7" "8"
[20] "9" "0" "1" "2" "5" "4" "3" "6" "7" "8" "9"


回りくどい仕事のように見えるが、対象が大きくなるほど早さを実感する。