compiler はなかなか優れもの

簡単に使えて，効果もそれなりに見込める

> func <- 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)
+ }
> system.time(a <- func(1:100, 10:200))
   ユーザ   システム       経過 
     0.160      0.001      0.159
> library(compiler)
> func2 <- cmpfun(func)
> system.time(b <- func2(1:100, 10:200))
   ユーザ   システム       経過 
     0.054      0.002      0.065 # 3 倍速い
> all.equal(a, b)
[1] TRUE

cxxfunction 使えそうなものを作ってみた

ColMeans, ColSums はあるが，colMins はないので，簡単だし作ってみた。

> src <- '
+ Rcpp::NumericMatrix x(ax);
+ int nc = x.ncol();
+ int nr = x.nrow();
+ Rcpp::NumericVector res(nc);
+ for (int j = 0; j < nc; j++) {
+     double res1 = x(0, j);
+     for (int i = 1; i < nr; i++) {
+         if (res1 > x(i, j)) res1 = x(i, j);
+     }
+     res[j] = res1;
+ }
+ return res;
+ '
>
> colMins <- cxxfunction(signature(ax="numeric"), src, plugin="Rcpp")
> colMins(matrix(1:24, 4))
[1]  1  5  9 13 17 21
> x <- matrix(rnorm(100000000), 10000)
> system.time(a <- apply(x, 2, min))
   ユーザ   システム       経過 
     2.776      1.205      3.966
> system.time(b <- colMins(x))
   ユーザ   システム       経過 
     0.265      0.000      0.265
> all.equal(a, b)
[1] TRUE

まあ，apply で十分速いという評価で。

cxxfunction std::transform を使う例

> src <- '
+     Rcpp::List input(data);
+     Rcpp::Function f(fun);
+     Rcpp::List output(input.size());
+     std::transform(input.begin(), input.end(), output.begin(), f);
+     output.names() = input.names();
+     return output;
+ '
> cpp_lapply <- cxxfunction(signature(data="list", fun="function"), src, plugin="Rcpp")
> (a <- as.data.frame(matrix(1:12, 3, 4)))
  V1 V2 V3 V4
1  1  4  7 10
2  2  5  8 11
3  3  6  9 12
> cpp_lapply(a, summary)
$V1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    1.0     1.5     2.0     2.0     2.5     3.0

$V2
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    4.0     4.5     5.0     5.0     5.5     6.0

$V3
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    7.0     7.5     8.0     8.0     8.5     9.0

$V4
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   10.0    10.5    11.0    11.0    11.5    12.0

> cpp_lapply(a, min)
$V1
[1] 1

$V2
[1] 4

$V3
[1] 7

$V4
[1] 10

cxxfunction 関数を呼ぶ例

inc <- '
    int func2(int n) {
        return n*2;
    }
'
src <- '
    Rcpp::IntegerVector N(n);
    Rcpp::IntegerVector ANS(1);
    ANS = func2(N(0)); /* [ ] でも ( ) でもよい */
    return ANS;
'
func <- cxxfunction(signature(n = "numeric"), body = src, inc=inc, plugin="Rcpp")

func(5)

Rcpp におけるクラス

スカラー変数というのはない（長さ 1 のベクトル）
Integer*, Numeric*, Logical*, Character*
　　　　　　　* は，Vector, Matrix

他に，List など

> src <- '
+ Rcpp::NumericVector x(xs);
+ Rcpp::NumericVector y(ys);
+ return Rcpp::wrap(ifelse(x < y, x*x, -(y*y)));
+ '
> a <- cxxfunction(signature(xs="numeric", ys="numeric"), src, plugin="Rcpp")
> a(rnorm(10), rnorm(10))
 [1]  1.04121737  1.72515331  0.24451134 -0.51313099  0.00310971  0.10984876 -0.37575399  0.22136341  0.90569395 -0.02640077

Rcpp: でなくて，こんな風にも書けるようだ。

> func <- cxxfunction(signature(xs="numeric", ns="integer"),
+     'int n = as<int>(ns);
+     double x = as<double>(xs);
+     double res;
+     res = pow(x, n);
+     return wrap(res);
+     ',
+     plugin="Rcpp")
> func(8.2, 3)
[1] 551.368

cxxfunction 行列を引数とする関数の例

行列を引数とする関数の例

> library(inline)
> src <- '
+ Rcpp::NumericMatrix xcpp(x);
+ int nrows = xcpp.nrow();
+ int ncolumns = xcpp.ncol();
+ for (int i = 0; i < nrows; i++) {
+     for (int j = 0; j < ncolumns; j++) {
+         xcpp[nrows * j + i] *= 2;
+     }
+ }
+ return xcpp;
+ '
> doublematrix <- cxxfunction(signature(x = "numeric"), body = src, plugin="Rcpp")
> print(doublematrix(matrix(1:10, nrow = 2)))
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    6   10   14   18
[2,]    4    8   12   16   20

> src1 <- '
+     Rcpp::NumericMatrix xbem(xbe);
+     int nrows = xbem.nrow();
+     Rcpp::NumericVector gv(g);
+     for (int i = 1; i < nrows; i++) {
+         xbem(i,_) = xbem(i-1,_) * gv[0] + xbem(i,_);
+     }
+     return xbem;
+ '
> funCPP1 <- cxxfunction(signature(xbe = "numeric", g="numeric"), body = src1, plugin="Rcpp")
> A <- matrix(1:12, 3, 4)
> funCPP1(A, 0.5)
     [,1] [,2]  [,3] [,4]
[1,] 1.00  4.0  7.00   10
[2,] 2.50  7.0 11.50   16
[3,] 4.25  9.5 14.75   20

> src <- '
+     Rcpp::NumericMatrix input(x);
+     int nc = input.ncol();
+     Rcpp::NumericMatrix output = clone<NumericMatrix>(input);
+     for (int i = 1; i < nc; i++) {
+         output.column(i) = output.column(i-1)+input.column(i);
+     }
+     return output;
+ '
> row.cumsum <- cxxfunction(signature(x="numeric"), src, plugin="Rcpp")
> (a <- matrix(1:12, 3, 4))
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> row.cumsum(a)
     [,1] [,2] [,3] [,4]
[1,]    1    5   12   22
[2,]    2    7   15   26
[3,]    3    9   18   30

様々な言語によるコードをインラインで 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

ダメ出し：行列・配列を無駄に使ってはいけない

http://d.hatena.ne.jp/KABIRA/

# 平行移動と拡大縮小で(0,1)におさめる
U<-V
U[1,]<-(V[1,]-min(V[1,]))/(max(V[1,])-min(V[1,]))
U[2,]<-(V[2,]-min(V[2,]))/(max(V[2,])-min(V[2,]))
U[3,]<-(V[3,]-min(V[3,]))/(max(V[3,])-min(V[3,]))

というのがあるが， U <- V は必要な大きさの行列を確保するためだけで，べつに V を U に代入しているわけではない。その後の3行も，無駄な記述である。添え字が違うだけだ。もし行数が100行あっても，100行のプログラムを書くわけではあるまい。3行だから，for 文で書いても行数はたいして変わらないということでこんな風に書いたのだろう。適切に書くと以下のようになる。

U <- t(apply(V, 1, function(x) (x-min(x)) / diff(range(x))))

U のメモリー確保も必要ないし，添え字を書く必要もないので，実にすっきり書ける。

ポイント：「明示的に」繰り返しを記述しなければならないということは，無駄なことの証明

ダメ出し：三角ダイアグラムを描く（Drawing ternary diagrams）

http://staff.aist.go.jp/a.noda/programs/ternary/ternary.html

ベクトル化されている関数はその機能を使おう。短いプログラムは分かりやすいとは限らないが，長いプログラムが分かりづらいのは確か。

####################
#
# Ternary.plot
#
# Drawing ternary plot of compositional data
# (modified from triangle.plot in ade4)
#
# Noda Atsushi
# Last updated: 2008/09/17 16:27:46.

# 0.1.2: axes. 軸（三角形の辺）を描くかどうかを選択できるようにした．
# 0.1.1: プロットされる点の型 (pch) と色(col) と大きさ (cex) を選択できるようにした．
# 0.1: just beggining.

ternary.plot <- function (dat, label=rownames(dat), axes=TRUE, clabel=0,
    cpoint=1, cline=0, lwd=1, addline=TRUE, label.vertices=TRUE,
    cline.col="black", pch=par("pch"), col=par("col"), cex=par("cex")) {
 
    par.orig <- par(mar=par("mar"))
    on.exit(par(par.orig))

    nam <- colnames(dat)
    dat <- dat/rowSums(dat)
   
    xy <- cbind((dat[,1]+2*dat[,3])/2, sqrt(3)*dat[,1]/2)

    plot(0, 0, type="n", xlim=c(0, 1), ylim=c(-0.1, 1), xlab="",
         ylab="", xaxt="n", yaxt="n", asp=1, frame.plot=FALSE)
   
    if (addline) {
        int <- 10
        i <- 1:(int-1)/int
        i.2 <- i/2
        sqrt3.2.i <- sqrt(3)*i.2
        segments(1-i.2,  sqrt3.2.i,  1-i, 0, col="lightgrey")
        segments(i.2,  sqrt3.2.i,  i, 0, col="lightgrey")
        segments(1-i.2,  sqrt3.2.i,  i.2, sqrt3.2.i, col="lightgrey")
    }

    if (axes)
        polygon(c(0.5, 0, 1), c(sqrt3.2, 0, 0), lwd=1)
    if (label.vertices)
        text(c(0.5, 0, 1), c(sqrt3.2, 0, 0), labels=nam, cex=1.5, pos=c(3, 1, 1))
    if (cpoint > 0)
        points(xy, pch=pch, cex=cex*cpoint, col=col)
    if (cline > 0)
        lines(xy, lwd=1, col=cline.col)
    if (clabel > 0)
        text(xy[, 1], xy[, 2], label, clabel)

    return(invisible(xy))
}