様々な言語によるコードをインラインで R の関数にする

cfunction は，様々な言語によるコードをインラインで R の関数にしてくれる。
対応している言語は "C++", "C", "Fortran", "F95", "ObjectiveC", "ObjectiveC++"。
FORTRAN も含まれているのが私にとっては懐かしい。

code <- "
      integer i
      do 10 i = 2, n(1)
        x(i) = x(i-1)+x(i)
   10 continue
"

fun <- cfunction(signature(n="integer", x="numeric"), code, convention=".Fortran")

ただ，計算内容によっては，R の関数の方が速い場合があるので，注意が必要。

実行結果

> n <- 10000000
> x <- rnorm(n)

> system.time(a <- cumsum(x))
   ユーザ   システム       経過 
     0.104      0.000      0.104

> system.time(b <- fun(n, x)$x)
   ユーザ   システム       経過 
     0.180      0.003      0.182

> all.equal(a, b)
[1] TRUE

前の記事のプログラミング例を FORTRAN でやってみると，計算時間は同じくらいになる。

> code <- "
+       integer i, j
+       do 10 i = 1, nx(1)
+           do 20 j = 1, ny(1)
+               z(i+j-1) = z(i+j-1)+x(i)*x(j)
+    20   continue
+    10 continue
+ "
 
> fun3 <- cfunction(signature(nx="integer", x="integer", ny="integer", y="integer", nz="integer", z="integer"), code, convention=".Fortran")

> system.time(ans <- fun3(1000, x, 1000, y, 1999, z))
   ユーザ   システム       経過 
     0.001      0.000      0.002

> all.equal(ans$z, a)
[1] TRUE

C++ コードをインラインで R の関数にする

inline パッケージ（関連：Rcpp パッケージ）を使って，C++ で書いたプログラムをインラインで R の関数にすることができる。
速度は素晴らしい（他のやり方は試していないので比較できないが，R で書いたプログラムよりは格段に速い）。

以下のように使う。

# C++ プログラムを cxxfunction で，R の関数にする
library(inline)
src <- '
    Rcpp::NumericVector xa(a);
    Rcpp::NumericVector xb(b);
    int n_xa = xa.size(), n_xb = xb.size();
    Rcpp::NumericVector xab(n_xa + n_xb - 1);
    for(int i = 0; i < n_xa; i++)
        for (int j = 0; j < n_xb; j++)
             xab[i + j] += xa[i] * xb[j];
    return xab;
'
fun <- cxxfunction(signature(a = "numeric", b = "numeric"), src, plugin = "Rcpp")

# R で関数を定義してみる
fun2 <- function(xa, xb) {
    xab <- integer(length(c(xa, xb))-1)
    for (i in seq.int(xa))
        for (j in seq.int(xb))
          xab[i+j-1] <- xab[i+j-1]+xa[i]*xb[j]
    return(xab)
}

実行結果は以下の通り

> x <- 1:1000
> y <- 1:1000

> system.time(a <- fun(x, y))
   ユーザ   システム       経過 
     0.003      0.000      0.003

> system.time(b <- fun2(x, y))
   ユーザ   システム       経過 
     6.298      0.053      6.249

> all.equal(a, b)
[1] TRUE

ダメ出し：Rで円周率計算

entertainment-lab
Rで円周率計算
http://entertainment-lab.blogspot.com/2011/06/r.html

> system.time({
+ set.seed(1)
+ points.num <- 100000
+ coords <- matrix(runif(2*points.num),ncol=2)
+ dist.from.origin <- apply(coords,1,function(x){sqrt(sum(x^2))})
+ myPI <- sum(dist.from.origin<1)/points.num*4
+ myPI
+ })
   ユーザ   システム       経過 
     0.944      0.026      0.966

rowSums, mean を使う。sqrt は不要。余計なことをしなければ速度は上がる。

> system.time({
+ set.seed(1)
+ points.num <- 100000
+ coords <- matrix(runif(2*points.num),ncol=2)
+ myPI <- mean(rowSums(coords^2) < 1)*4
+ myPI
+ })
   ユーザ   システム       経過 
     0.011      0.001      0.014

ダメ出し：ryamadaのRソース管理帳

http://d.hatena.ne.jp/ryamada22222/ の中の MakeOTable 関数は以下のようになっているが

MakeOTable<-function(p,k=c(0,1)){
    n<-length(k)
    m<-matrix(0,n,length(p))
    for(i in 1:n){
        m[i,which(p==k[i])]<-1
    }
    m
}

こんな風にした方がよい。なんと，1行でできる。

MakeOTable2 <- function(p, k=0:1) {
    sapply(p, "==", k)+0
}

R では，== も「関数」なのだ。

>"=="(3, 4)
[1] FALSE
> "=="(3, 3)
[1] TRUE