2015-01-06

RstudioでC++を覚える

Rstudio R C++ Rcpp

2015-01-06

準備と動作確認 Preparation 〜RstudioでC++を覚える

Rstudio R C++ Rcpp

Rをインストール Install R
Rstudioをインストール Install Rstudio
RstudioのツールバーメニューからC++ファイルを新規作成、とすると、もしコンパイラとかの準備がまだだったらdevtoolsとかのインストールをするかどうか訊いてきてくれるので、それにしたがうのもよい
- コンパイラ(python用に入れたこれでよいの？)、Rtoolsかと迷ったらこちら、？ You need C++ compiler. Rtools is one way to get an appropriate one.
Rstudioを立ち上げる Start Rstudio
これも
これも
RコンソールでRcppパッケージをインストールして読み込み Install Rcpp package and load it in the R console panel of Rstudio.

install.packages("Rcpp")
library(Rcpp)

Rstudioでc++ファイルを作る Make a c++ file with Rstudio.
- RstudioのツールバーからC++ファイルを新しく作る、にするとパネルが開くので、最初はサンプルをそのまま貼る (こちら) Go to Rstudio's toolbar and select "make a new c++ file" and paste sample codes into it.
- "gibbs.cpp"という名前で保存する(ワーキングディレクトリに) Save it "gibbs.cpp" in your working directory.

#include <Rcpp.h>
using namespace Rcpp;
 
// [[Rcpp::export]]
NumericMatrix gibbs(int N, int thin) {
 
   NumericMatrix mat(N, 2);
   double x = 0, y = 0;
 
   RNGScope scope;
   for(int i = 0; i < N; i++) {
      for(int j = 0; j < thin; j++) {
         x = R::rgamma(3.0, 1.0 / (y * y + 4));
         y = R::rnorm(1.0 / (x + 1), 1.0 / sqrt(2 * x + 2));
      }
      mat(i, 0) = x;
      mat(i, 1) = y;
   }
 
   return(mat);
}

RstudioのRコンソールに戻って Go back to R console in Rstudio then,

# c++ファイルを読み込んで
sourceCpp("gibbs.cpp")
# c++ファイルに作ってあるgibbs()という関数を使う
gibbs(10,10)

こうなる You are seeing this

2015-01-06

練習 Practice 〜RstudioでC++を覚える

Rstudio R C++ Rcpp

こちらをいくつかやってみる Try some in this page.
Vector_Subsetting
- "vectorsubsetting.cpp"という名前で保存 Save below as "vectorsubsetting.cpp"

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector positives(NumericVector x) {
    return x[x > 0];
}

// [[Rcpp::export]]
List first_three(List x) {
    IntegerVector idx = IntegerVector::create(0, 1, 2);
    return x[idx];
}

// [[Rcpp::export]]
List with_names(List x, CharacterVector y) {
    return x[y];
}

// [[Rcpp::export]]
NumericVector in_range(NumericVector x, double low, double high) {
    return x[x > low & x < high];
}

// [[Rcpp::export]]
NumericVector no_na(NumericVector x) {
    return x[ !is_na(x) ];
}

bool is_character(SEXP x) {
    return TYPEOF(x) == STRSXP;
}

// [[Rcpp::export]]
List charvecs(List x) {
    return x[ sapply(x, is_character) ];
}

Rで読んで、実行する Read it into R and use it.

install.packages("microbenchmark")
library(microbenchmark)
sourceCpp("vectorsubsetting.cpp")
R_in_range <- function(x, low, high) {
    return(x[x > low & x < high])
}
x <- rnorm(1E5)
identical( R_in_range(x, -1, 1), in_range(x, -1, 1) )
microbenchmark( times=5, 
    R_in_range(x, -1, 1),
    in_range(x, -1, 1)
)

- うまく動いたら、何をしているか確認してみる。It runs ok? then read codes to understand how it works.
C++内でRの関数を使ってみる Use a R function in C++
- R function

x <- rnorm(1000)
fivenum(x) # min,max,quantiles(x,c(0.25,0.5,0.75))

- "callfunction.cpp"

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector callFunction(NumericVector x, Function f) {
    NumericVector res = f(x);
    return res;
}

- 使う Use it

sourceCpp("callfunction.cpp")
callFunction(x,fivenum)
fivenum(x)

> sourceCpp("callfunction.cpp")
> callFunction(x,fivenum)
[1] -2.96348996 -0.70554582 -0.01415115  0.66923375  3.05126692
> fivenum(x)
[1] -2.96348996 -0.70554582 -0.01415115  0.66923375  3.05126692

2015-01-06

自分のためにRcppを使ってみる〜RstudioでC++を覚える

Rstudio R C++ Rcpp

だいたいわかった。C++は速いが、ごつい。便利な関数もあるにはあるが、Rやpythonのようにきめ細かくない。…ということが
もう一つ知っていることがある。自分はC++に慣れていないということ。→速くできるのかもしれないけれど、自分ができないなら、それは「速くできないこと」だということ。
さて、今、手元にある「重くて遅いRでの処理」のどこにC++を入れて速くするか、試行錯誤してみることにする。
ここでの課題は「統計遺伝学」というよりも「離散数学」なので、こちらに書いておくことにする。

2015-01-06

Armadilloを使う〜RstudioでC++を覚える

Rstudio R C++ Rcpp

Rcppパッケージの基本、それをRstudioで扱うことの基本はわかった
けれど、「素のC++」だけでは、ごつすぎて使いにくい。特に線形代数とか、いわゆる数値計算のためのベクトル・行列・アレイなどが使いやすくて、それの演算ライブラリがないとだめなので、Armadilloを入れて、それを使えるようになることにする
そもそもArmadilloはC++用の線形代数ライブラリ(こちら)
それをRcppで使うようにしたのがRcppArmadilloパッケージなので、以下をして使えるようにする

install.packages("RcppArmadillo")
library(RcppArmadillo)

RcppArmadilloのチュートリアルとかもあるけれど、手っ取り早く、Rstudioで動かしているサンプルを見つけよう→こちら
- 結局、以下を"LinearRegression.cpp"という名前で保存して(これより上にヘッダーは不要(#include
  も不要)
  - 逆行列計算ができていたり(solve)、行列の転置ができていたり(mX.t())、行列の積ができていたり(mX.t() * mX)するところがArmadilloライブラリの恩恵。List型でRに返したりしているところはRcppArmadilloパッケージの実装の恩恵

// LinearRegression.cpp
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace arma;  // use the Armadillo library for matrix computations
using namespace Rcpp;

// [[Rcpp::export]]
List fnLinRegRcpp(vec vY, mat mX) {

  // compute the OLS estimator & model residuals
  vec vBeta =  solve(mX.t()*mX, mX.t()*vY);  
  vec vResid = vY - mX * vBeta;

  // construct the return object
  List ret;
  ret["beta"] = vBeta;
  ret["resid"] = vResid;

  return ret;
}

// END

- R側でこのファイルを読み込んでコンパイルして("sourceCpp("hoge.cpp")"の代りに"Rcpp::sourceCpp("hoge.cpp")"

library(devtools)
library(Rcpp)
library(RcppArmadillo)
Rcpp::sourceCpp("LinearRegression.cpp", showOutput = TRUE, rebuild = FALSE) # もちろんファイルの置き場所はこれでアクセスできるように…

- 関数を使う

# generate some sample data
iK = 4
iN = 100
mX = cbind(1, matrix(rnorm(iK*iN), iN, iK))
vBeta0 = c(2, 3.5, 0.11, 6.33, 23)
vY = rnorm(iN, mean = mX %*% vBeta0)

# test the function
linReg1 = fnLinRegRcpp(vY, mX)
linReg1$beta  # coefficient estimates

# compare the results to the built-in lm.fit function
lm.fit(y = vY, x = mX)$coef  # coefficient estimates

# END

2015-01-06

C++を？

R C++

「C++かCを覚えてみる？」という話になった
一からやってもよいけれど、多分、飽きてしまうし、Rと連携して使うことが多そうだから、Rcppを使うことを前提にC++を覚えてみる、例題は自分たちがRで書いているコードをC++化する、ということにする、というようなやり方ではどうだろうか…
こんな文書があった。これをなぞるかどうか考えてみる
Rcpp Galleryも検討する
Rstudioを使うならRcpp with Rstudio(En),RstudioでRcpp...(Jp)→RstudioでC++を覚える(考える暇があったらやってみよう)
その他資料