5hun.tumblr.com

William Stafford Nobleの"A Quick Guide to Organizing Computational Biology Projects"を和訳した。

計算生物学という生物学の分野の研究者に向けて、計算機実験におけるファイルの管理方法やスクリプトの書き方について書かれている。しかし、計算生物学に限らず、色々なパラメータでスクリプトを走らせてベンチマークを取ったり、データ分析を行ったりするような、計算機上で実験を行う人には参考になる部分があると思う。

boost::pythonを使用する際に参考にした記事や、サンプルコードを載せておく。

全般的に参考となるドキュメント

boost::pythonの利用に際して全般的に参考になる情報源である。公式のリファレンスの他、wiki.python.orgの記事が良かった。

boostのリファレンス

wiki.python.orgの記事

python の iterable なオブジェクトと c++ の vector などの暗黙的な変換

pythonのリスト・タプルなどの iterable なオブジェクトと c++ の std;:vector のような互換性がありそうなクラスのオブジェクトを暗黙に変換する方法が下記で紹介されている。

http://stackoverflow.com/questions/15842126/feeding-a-python-list-into-a-function-taking-in-a-vector-with-boost-python

C++の抽象(Abstract)クラスを使用する場合

wiki.python.orgのclassの記事のAbstract Classの欄

抽象クラスを継承する抽象クラスの場合はこんな感じ

C++の仮想関数をpython側に公開したい場合

boostの公式ドキュメント

Non-virtual interface idiomを使用する方法

関数のオーバーロード

wiki.python.orgの記事

boost::serializeを用いてpickle処理を実装する

Boost.Python Pickle Support

boost::python を使用して仮想関数を持つクラスを python 側に公開する場合、ラッパークラスを用意する必要がある。（詳しくはboostのドキュメントなどを参照。）

しかし、Non-virtual interface (NVI) と呼ばれるイディオムを使用すれば、ラッパークラスが不要になる。

また、NVIは boost::python とは関係なく使用されているイディオムでもある。詳しくは以下などを参照。

Prolog Programming for Artificial Intelligenceの第4版を読んでいる。演習問題の7.8をやってみた結果をメモしておく。

問題の概要は、N-Queen問題をFinite Domain上のCLP（Constraint Logic Programming）ライブラリを用いて解くときに、盤上の左の方のQueenから位置を試していくより、真ん中の方を優先した方が速くなるので試してみよというものだ。

通常、PrologのCLPにおいてはlabelingという述語を用いて、変数を確定する順番を定める。例えばlabeling([], [X1,X2, ... ,Xn])であれば、左のほうから（X1、X2、…の順に）値を試していく。labelingの第１引数は順番の決め方を指定するもので、[]だと左からになる。使用していた処理系SWI-Prologには、真ん中の変数からなどというオプションはなかった。仕方ないので、真ん中の方の変数が最初に来るように変数のリストを並べ替えて、labeling([], ...)を使うようにした。

並び替え部分のコードは次の通り。

% リスト Xs を（可能な限り）同じ長さで 2 つのリストに分けた時の、 %% 各リストの長さを求める。 %% 例: halflengths([1,2,3,4,5], 3, 2). halflengths(Xs, L1, L2) :- length(Xs, L), R is L mod 2, L1 is L // 2 + R, L2 is L // 2. %% リスト Xs を真ん中で２つのリストに分ける。 %% Xs の長さが奇数の場合は１つ目のリストが 1 つ分長くなる。 %% 例: half([1,2,3], [1,2], [3]). half(Xs, Y1, Y2) :- halflengths(Xs, L1, L2), length(Y1, L1), length(Y2, L2), append(Y1, Y2, Xs). %% ２つのリストを交互にマージしたようなリストを考える。 %% 例: merge([1,2], [3,4], [1,3,2,4]). merge([], X, X). merge(X, [], X). merge([X|Xs], [Y|Ys], [X,Y|Zs]) :- merge(Xs, Ys, Zs). %% リスト Xs から、真ん中の元が先に来るようなリスト Ys を構築 %% 例: middlefirst([1,2,3,4,5], [3,4,2,5,1]). middlefirst(Xs, Ys) :- half(Xs, H1, H2), reverse(H1, R1), merge(R1, H2, Ys).

labelingを行う前に上記middlefirstによる並べ替え処理を入れて実行した場合のベンチマークをとってみた結果が以下である。leftfirstというのが並べ替え処理を入れない場合、midfirstは入れた場合。さらに、labelingの第１引数に[ff]を指定した場合も載せてある。ここで、実行時間としては、最初の解を見つけるまでのCPU時間を述語statisticsで３回繰り返し測った平均を使用している。

google.load("visualization", "1", {packages:["corechart"]}); google.setOnLoadCallback(drawChart); function drawChart() { var data = google.visualization.arrayToDataTable([ ['N','leftfirst','midfirst','leftfirst ([ff])','midfirst ([ff])'], [4,0.666666666666667,0.666666666666667,0.666666666666667,0.666666666666667], [5,0.666666666666667,0.666666666666667,1,0.666666666666667], [6,2,3,3,2.33333333333333], [7,1.33333333333333,1,1.33333333333333,1], [8,6.33333333333333,6.33333333333333,6.66666666666667,6.66666666666667], [9,4,3.66666666666667,4.66666666666667,4], [10,9.33333333333333,9,6.33333333333333,5.33333333333333], [11,5.33333333333333,5.33333333333333,9,8], [12,20,21,20.6666666666667,19.3333333333333], [13,10.3333333333333,11.3333333333333,27.3333333333333,28.3333333333333], [14,135,133.333333333333,19,18.6666666666667], [15,106,105.333333333333,7.66666666666667,7.33333333333333], [16,738.333333333333,719.333333333333,10.3333333333333,10], [17,490.666666666667,476.666666666667,13,11.6666666666667], [18,3213.33333333333,3178.66666666667,13.6666666666667,14.3333333333333], [19,258.333333333333,258.333333333333,14.3333333333333,15], [20,16777,16789,28.3333333333333,30.3333333333333], [21,804.666666666667,790.333333333333,13.6666666666667,14.3333333333333], [22,138128.666666667,137760.333333333,15,13], [23,2343,2342.66666666667,21.6666666666667,22.3333333333333], [24,36592.6666666667,36504.3333333333,18.3333333333333,18.3333333333333], [25,4689.33333333333,4617,53.3333333333333,57.6666666666667] ]); var options = { title: 'Benchmark Results', vAxis: {title: "CPU Time (ms)", logScale: true}, hAxis: {title: "N"} }; var chart = new google.visualization.LineChart(document.getElementById('nqueen_chart_div')); chart.draw(data, options); }

Rには向いていない処理だと思うが、シミュレーション的処理をRでやることを考える。その際クラスオブジェクトのメンバを繰り返し書き換えるような処理が出てくる。

そのような処理において、R5クラスが速くなるかということが気になり、単純なプログラムでベンチマークを取ってみた。

コードは以下の通り。

## S3 class inc <- function(x) UseMethod("inc") inc.hoge <- function(h){ h$a = h$a + 1 return(h) } test.s3 = function(num=1000){ x = list(a=0) class(x) = c('hoge', 'list') for(i in seq(num)){ x = inc(x) } cat(x$a) } ## S4 class setClass("HogeS4", representation(a="numeric")) setMethod("inc", c("HogeS4"), function(x){ x@a = x@a + 1 return(x) }) test.s4 = function(num=1000){ x = new("HogeS4", a=0) for(i in seq(num)){ x = inc(x) } cat(x@a) } ## R5 class r5cls = setRefClass(Class="HogeR5", fields = list(a='numeric'), methods = list( inc = function(){ a <<- a + 1 })) test.r5 <- function(num=1000){ x = r5cls$new(a=0) for(i in seq(num)){ x$inc() } cat(x$a) }

ベンチマーク結果は次の通り。グラフにもしてみた。

> system.time(test.s3(1000000)) 1e+06 user system elapsed 27.581 0.012 27.628 > system.time(test.s4(1000000)) 1e+06 user system elapsed 34.882 0.000 34.924 > system.time(test.r5(1000000)) 1e+06 user system elapsed 27.698 0.000 27.728

google.load("visualization", "1", {packages:["corechart"]}); google.setOnLoadCallback(drawChart); function drawChart() { var data = google.visualization.arrayToDataTable([ ['Type', 'user', 'system', 'elapsed'], ['S3', 27.581, 0.012, 27.628], ['S4', 34.882, 0.000, 34.924], ['R5', 27.698, 0.000, 27.728] ]); var options = { title: 'Benchmark Results' }; var chart = new google.visualization.ColumnChart(document.getElementById('chart_div')); chart.draw(data, options); }

Rの可変長引数...を用いていると、以下のようにmissingがうまく働かないことがある。

f <- function(a){ if(missing(a)){ print("missing") }else{ print("not missing") } } g <- function(a){f(a)} g2 <- function(a){g(a)} h <- function(...){f(...)} h2 <- function(a){h(a)}

各関数を引数を与えずに実行すると次のような結果になる。

> g() [1] "missing" > h() [1] "missing" > g2() [1] "missing" > h2() [1] "not missing"

引数を与えていないのだから、h2()の結果は"missing"になってほしいところだ。

…を用いずに明示的に引数名を記述すれば上記の問題は避けられる。しかし、与えられた引数を他の関数にただ受け渡すだけの場合は、上記のように透過的に（呼び出す関数のインターフェイスに依存せず）書きたい。

引数が与えられていないことを正しく判定するには、例えば次のような関数をmissingの代わりに用いる方法がある。実際に評価してみてエラーとなるか否かにより判定している。

missing2 <- function(x){ t <- try(force(x), silent=TRUE) return(inherits(t, "try-error")) }

data.frame型を継承したS4クラスのオブジェクトに対し，そのまま代入操作を行うとエラーになる．

> setClass("hoge", contains = c("data.frame")) [1] "hoge" > ho <- new("hoge", data.frame(a=runif(3), b=runif(3))) > ho Object of class "hoge" a b 1 0.1446811 0.4779831 2 0.7206756 0.4498786 3 0.5258015 0.7642846 > ho[1,2] <- 100 警告メッセージ： In class(x) <- NULL : class(x) を NULL に設定しようとしています；結果はもはや S4 オブジェクトではなくなります > ho a b 1 0.1446811 0.4779831 2 0.7206756 0.4498786 3 0.5258015 0.7642846

data.frame型の代入操作（"[<-"函数）で，クラス属性を弄っていることが原因のようである．回避策としては，S3Partを用いる方法がある．

> ho <- new("hoge", data.frame(a=runif(3), b=runif(3))) > S3Part(ho, strict=T)[1,2] <- 100 > ho Object of class "hoge" a b 1 0.8765066 0.2053294 2 0.3706098 0.9298724 3 0.4628688 0.3079220 > S3Part(ho, strict=T)[1,2] <- 200 > ho Object of class "hoge" a b 1 0.8765066 200.0000000 2 0.3706098 0.9298724 3 0.4628688 0.3079220

代入函数をオーバーロードする方法もある．次のようにオーバーロードを行う．

setReplaceMethod( "[", "hoge", function(x, i, j, ..., value){ if(nargs() == 3){ if(missing(i)){ S3Part(x, strict=T)[j] <- value }else{ S3Part(x, strict=T)[i] <- value } }else{ if(missing(i)){ if(missing(j)){ S3Part(x, strict=T)[,] <- value }else{ S3Part(x, strict=T)[,j] <- value } }else{ if(missing(j)){ S3Part(x, strict=T)[i,] <- value }else{ S3Part(x, strict=T)[i,j] <- value } } } x })

するとエラーは出なくなる．

> ho <- new("hoge", data.frame(a=runif(3), b=runif(3))) > ho Object of class "hoge" a b 1 0.1446811 0.4779831 2 0.7206756 0.4498786 3 0.5258015 0.7642846 > ho[1,2] <- 200 > ho Object of class "hoge" a b 1 0.1446811 200.0000000 2 0.7206756 0.4498786 3 0.5258015 0.7642846

同様のエラーはho[["a"]] <- c(1,2,3)のように"[["でも起こる（"$"では大丈夫なようだ）．S3Partを用いるか，次のようにオーバーロードしたらよい，

setReplaceMethod( "[[", "hoge", function(x, i, j, ..., value){ if(missing(j)){ S3Part(x, strict=T)[[i]] <- value }else{ S3Part(x, strict=T)[[i,j]] <- value } x })

以下は，Rのreshape函数のHelpにあるexample．但し一部変更してある．

> head(Indometh) Subject time conc 1 1 0.25 1.50 2 1 0.50 0.94 3 1 0.75 0.78 4 1 1.00 0.48 5 1 1.25 0.37 6 1 2.00 0.19 > wide <- reshape(Indometh, v.names="conc", idvar="Subject", timevar="time", direction="wide") > head(wide) Subject conc.0.25 conc.0.5 conc.0.75 conc.1 conc.1.25 conc.2 1 1 1.50 0.94 0.78 0.48 0.37 0.19 12 2 2.03 1.63 0.71 0.70 0.64 0.36 23 3 2.72 1.49 1.16 0.80 0.80 0.39 34 4 1.85 1.39 1.02 0.89 0.59 0.40 45 5 2.05 1.04 0.81 0.39 0.30 0.23 56 6 2.31 1.44 1.03 0.84 0.64 0.42 conc.3 conc.4 conc.5 conc.6 conc.8 1 0.12 0.11 0.08 0.07 0.05 12 0.32 0.20 0.25 0.12 0.08 23 0.22 0.12 0.11 0.08 0.08 34 0.16 0.11 0.10 0.07 0.07 45 0.13 0.11 0.08 0.10 0.06 56 0.24 0.17 0.13 0.10 0.09 > reshape(wide, direction="long") Subject time conc 1.0.25 1 0.25 1.50 2.0.25 2 0.25 2.03 3.0.25 3 0.25 2.72 4.0.25 4 0.25 1.85 5.0.25 5 0.25 2.05 ... 2.8 2 8.00 0.08 3.8 3 8.00 0.08 4.8 4 8.00 0.07 5.8 5 8.00 0.06 6.8 6 8.00 0.09 > reshape(wide, idvar="Subject", varying=list(2:12), v.names="conc", direction="long") Subject time conc 1.1 1 1 1.50 2.1 2 1 2.03 3.1 3 1 2.72 4.1 4 1 1.85 5.1 5 1 2.05 ... 2.11 2 11 0.08 3.11 3 11 0.08 4.11 4 11 0.07 5.11 5 11 0.06 6.11 6 11 0.09 > ## times need not be numeric > df <- data.frame(id=rep(1:4,rep(2,4)), visit=I(rep(c("Before","After"),4)), x=rnorm(4), y=runif(4)) > df id visit x y 1 1 Before 1.00602114 0.4236643 2 1 After 1.74354714 0.6772507 3 2 Before 0.40998887 0.9493135 4 2 After 0.09952564 0.8352482 5 3 Before 1.00602114 0.4236643 6 3 After 1.74354714 0.6772507 7 4 Before 0.40998887 0.9493135 8 4 After 0.09952564 0.8352482 > reshape(df, timevar="visit", idvar="id", direction="wide") id x.Before y.Before x.After y.After 1 1 1.0060211 0.4236643 1.74354714 0.6772507 3 2 0.4099889 0.9493135 0.09952564 0.8352482 5 3 1.0060211 0.4236643 1.74354714 0.6772507 7 4 0.4099889 0.9493135 0.09952564 0.8352482 > ## warns that y is really varying > reshape(df, timevar="visit", idvar="id", direction="wide", v.names="x") id y x.Before x.After 1 1 0.4236643 1.0060211 1.74354714 3 2 0.9493135 0.4099889 0.09952564 5 3 0.4236643 1.0060211 1.74354714 7 4 0.9493135 0.4099889 0.09952564 > ## unbalanced 'long' data leads to NA fill in 'wide' form > df2 <- df[1:7,] > df2 id visit x y 1 1 Before 1.00602114 0.4236643 2 1 After 1.74354714 0.6772507 3 2 Before 0.40998887 0.9493135 4 2 After 0.09952564 0.8352482 5 3 Before 1.00602114 0.4236643 6 3 After 1.74354714 0.6772507 7 4 Before 0.40998887 0.9493135 > reshape(df2, timevar="visit", idvar="id", direction="wide") id x.Before y.Before x.After y.After 1 1 1.0060211 0.4236643 1.74354714 0.6772507 3 2 0.4099889 0.9493135 0.09952564 0.8352482 5 3 1.0060211 0.4236643 1.74354714 0.6772507 7 4 0.4099889 0.9493135 NA NA > ## Alternative regular expressions for guessing names > df3 <- data.frame(id=1:4, age=c(40,50,60,50), dose1=c(1,2,1,2), dose2=c(2,1,2,1), dose4=c(3,3,3,3)) > reshape(df3, direction="long", varying=3:5, sep="") id age time dose 1.1 1 40 1 1 2.1 2 50 1 2 3.1 3 60 1 1 4.1 4 50 1 2 1.2 1 40 2 2 2.2 2 50 2 1 3.2 3 60 2 2 4.2 4 50 2 1 1.4 1 40 4 3 2.4 2 50 4 3 3.4 3 60 4 3 4.4 4 50 4 3 > ## an example that isn't longitudinal data > state.x77 <- as.data.frame(state.x77) > long <- reshape(state.x77, idvar="state", ids=row.names(state.x77), times=names(state.x77), timevar="Characteristic", varying=list(names(state.x77)), direction="long") > reshape(long, direction="wide") state Population Income Alabama.Population Alabama 3615 3624 Alaska.Population Alaska 365 6315 Arizona.Population Arizona 2212 4530 Arkansas.Population Arkansas 2110 3378 California.Population California 21198 5114 ... Illiteracy Life Exp Murder HS Grad Alabama.Population 2.1 69.05 15.1 41.3 Alaska.Population 1.5 69.31 11.3 66.7 Arizona.Population 1.8 70.55 7.8 58.1 Arkansas.Population 1.9 70.66 10.1 39.9 California.Population 1.1 71.71 10.3 62.6 ... Frost Area Alabama.Population 20 50708 Alaska.Population 152 566432 Arizona.Population 15 113417 Arkansas.Population 65 51945 California.Population 20 156361 ... > reshape(long, direction="wide", new.row.names=unique(long$state)) state Population Income Illiteracy Alabama Alabama 3615 3624 2.1 Alaska Alaska 365 6315 1.5 Arizona Arizona 2212 4530 1.8 Arkansas Arkansas 2110 3378 1.9 California California 21198 5114 1.1 ... Life Exp Murder HS Grad Frost Area Alabama 69.05 15.1 41.3 20 50708 Alaska 69.31 11.3 66.7 152 566432 Arizona 70.55 7.8 58.1 15 113417 Arkansas 70.66 10.1 39.9 65 51945 California 71.71 10.3 62.6 20 156361 ... > ## multiple id variables > df3 <- data.frame(school=rep(1:3,each=4), class=rep(9:10,6), time=rep(c(1,1,2,2),3), score=rnorm(12)) > wide <- reshape(df3, idvar=c("school","class"), direction="wide") > wide school class score.1 score.2 1 1 9 0.0381327 0.7346048 2 1 10 0.2816448 -0.1762236 5 2 9 -2.4207748 -0.4001436 6 2 10 0.2742648 1.7715353 9 3 9 -0.6904898 0.5204333 10 3 10 0.2139066 0.3323449 > ## transform back > reshape(wide) school class time score.1 1.9.1 1 9 1 0.0381327 1.10.1 1 10 1 0.2816448 2.9.1 2 9 1 -2.4207748 2.10.1 2 10 1 0.2742648 3.9.1 3 9 1 -0.6904898 3.10.1 3 10 1 0.2139066 1.9.2 1 9 2 0.7346048 1.10.2 1 10 2 -0.1762236 2.9.2 2 9 2 -0.4001436 2.10.2 2 10 2 1.7715353 3.9.2 3 9 2 0.5204333 3.10.2 3 10 2 0.3323449 警告メッセージ： In reshapeWide(data, idvar = idvar, timevar = timevar, varying = varying, : some constant variables (y) are really varying