Go to Qiita Advent Calendar 2024 Top
LoginSignup
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@_R2_

Rcpp, RcppArmadilloの練習(1)

Posted at

0. Rcpp,RcppArmadilloの環境構築

実行環境はMac．MacでRcpp,RcppArmadilloをインストールするのはだいぶめんどくさい．

  1. Rでinstall.packages("Rcpp", dependencies=TRUE)
  2. XcodeとCommand Line Toolをインストールしてから，ターミナルで「sudo xcodebuild -license」と打ち込んでライセンス条項にagreeする．XcodeのインストールはApp Storeから．Command Line Toolはこちらを参照．
  3. Rでinstall.packages("RcppArmadillo", type = "source")
  4. 3.で「library not found for -lgfortran
    clang: error: linker command failed with exit code 1 (use -v to see invocation)」みたいなエラーが出る．gfortanが入っていない場合はこちらからインストールする
  5. 4.のエラーが出るのにgfortanが入っているはずの場合は，Finderからメニューの「移動」→「フォルダへ移動」と選んで，「~/.R」と検索する．ここに「Makevars」というファイルがあればそれを開き，無ければテキストエディタで新規作成する．Makevarsファイルに入れたら，「FLIBS=-L/usr/local/gfortran」という行を加える．
  6. もう一度install.packages("RcppArmadillo", type = "source")を試し，エラーが出なければ完了．

1.行列の基本演算

以下，こちらを参考にまとめる．

// MySample1.cpp
#include <RcppArmadillo.h> // 呪文
// [[Rcpp::depends(RcppArmadillo)]] // 呪文
// [[Rcpp::plugins("cpp11")]] // 呪文
// [[Rcpp::export]] // 呪文
void fSample1() // 関数fSample1を定義。戻り値なし。
{
  
  double da, db, dc; // 実数 
  arma::mat mA, mB, mC, mD; // RcppArmadilloに依存する行列
  
  
  /*---実数の定義---*/
  da = 2.0;
  db = arma::datum::pi; // piとかeとか特別な数をdatumで定義
  dc = da + db;
  
  Rprintf("da=%.4f, db=%.4f\n",da,db); // 表示にはRprintfかRcoutを使う。有効数字4桁
  Rprintf("da+db=%.4f\n",da+db);
  Rcpp::Rcout << "dc=" << dc <<"\n"; // オブジェクトを<<でつなげる
  
  
  /*---行列の定義---*/
  mA = { {1, 2}, 
         {3, 4} };
  
  mA.print("mA="); // 行列の表示
  
  mB = arma::zeros(2,2); //ゼロ行列 
  
  mB(0,0) = 10;
  mB(0,1) = 20;
  mB(1,0) = 30;
  mB(1,1) = 40;
  mB.print("mB=");
  
  mC ={ {100,200}, 
        {300,400}};
  mC.print("mC=");
  
  
  /*---簡単な行列演算---*/
  (da*mA).print("a*mA=");
  (mA+mB).print("mA+mB=");
  (mA*mB).print("mAmB="); // 行列積

  mD = arma::inv(mA); // 逆行列
  mD.print("mA^{-1}=");
  
  
  /*---ベクトルの定義---*/
  arma::vec vx = {1, 4, 3, 2, 5};
  vx.print("vx=");
  
  
  /*---簡単なベクトル演算---*/
  Rprintf("Mean     = %.4f \n", arma::mean(vx));
  Rprintf("Variance = %.4f \n", arma::var(vx));
  Rprintf("Minimum  = %.4f \n", arma::min(vx));
  Rprintf("Maximum  = %.4f \n", arma::max(vx));
  
  sort(vx, "ascend").print("Sort in ascending order");
  sort(vx, "descend").print("Sort in descending order");

  
  /*---簡単な行列演算part2---*/
  arma::mat mX = { {1, 30}, {4, 20}, {3, 50}, {2, 40}, {5, 10} };
  mX.print("mX=");
  arma::mean(mX).print("Mean");
  arma::var(mX).print("Variance");
  arma::min(mX).print("Minimum");
  arma::max(mX).print("Maximum");
  
  sort(mX, "descend").print("Sort in descending order");
}

/*** R
library(Rcpp)
library(RcppArmadillo)
fSample1() 
*/

###出力###
da=2.0000, db=3.1416
da+db=5.1416
dc=5.14159
mA=
   1.0000   2.0000
   3.0000   4.0000
mB=
   10.0000   20.0000
   30.0000   40.0000
mC=
   1.0000e+02   2.0000e+02
   3.0000e+02   4.0000e+02
a*mA=
   2.0000   4.0000
   6.0000   8.0000
mA+mB=
   11.0000   22.0000
   33.0000   44.0000
mAmB=
   7.0000e+01   1.0000e+02
   1.5000e+02   2.2000e+02
mA^{-1}=
  -2.0000   1.0000
   1.5000  -0.5000
vx=
   1.0000
   4.0000
   3.0000
   2.0000
   5.0000
Mean     =   3.0000 
Variance =   2.5000 
Minimum  =   1.0000 
Maximum  =   5.0000 
Sort in ascending order
   1.0000
   2.0000
   3.0000
   4.0000
   5.0000
Sort in descending order
   5.0000
   4.0000
   3.0000
   2.0000
   1.0000
mX=
    1.0000   30.0000
    4.0000   20.0000
    3.0000   50.0000
    2.0000   40.0000
    5.0000   10.0000
Mean
    3.0000   30.0000
Variance
   2.5000e+00   2.5000e+02
Minimum
    1.0000   10.0000
Maximum
    5.0000   50.0000
Sort in descending order
    5.0000   50.0000
    4.0000   40.0000
    3.0000   30.0000
    2.0000   20.0000
    1.0000   10.0000

2. 簡単な行列の作り方

// MySample2.cpp
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
void fSample2() // 関数fSample2を定義，戻り値なし
{
  arma::mat mX, mY;
  arma::vec va, vb;
  
  
  /*---ゼロ行列，単位行列の定義---*/
  mX = arma::zeros(2,2);
  mY = arma::eye(2,2);
  
  mX.print("mX=");
  mY.print("mY=");
  
  
  /*---特別なベクトルの定義---*/
  va = arma::ones(2); // 1のみを要素とする長さ2のベクトル
  vb = arma::regspace(1,0.5,5); // 1から5まで,0.5刻みの等差数列
  
  va.print("va=");
  vb.print("vb=");
  (mY*va).print("mY*va=");
  
  
  /*---行列同士の操作---*/
  arma::vec vp, vq, vr;

  vp = {1, 2};
  vq = arma::zeros(2);
  vq[0] = 3;
  vq[1] = 4;

  vp.print("vp=");
  vq.print("vq=");
  
  arma::mat mP, mQ;
  
  mP = join_horiz(vp, vq); // 行方向に連結(2行2列になる)
  mQ = join_vert(vp, vq);  // 列方向に連結(4行1列になる)
  
  mP.print("mP=(vp vq)=");
  mQ.print("mQ=(vp' vq')'=");
  
  (vp.t()*vq).print("vp'vq="); // ベクトル積を計算.vp.t()で転置．．
  (arma::as_scalar(vp.t()*vq)*mP).print("(vp'*vq)*mP="); // as.scalarで一次元ベクトルを数値化
  
  }

/*** R
 fSample2()
*/
  

mX=
        0        0
        0        0
mY=
   1.0000        0
        0   1.0000
va=
   1.0000
   1.0000
vb=
   1.0000
   1.5000
   2.0000
   2.5000
   3.0000
   3.5000
   4.0000
   4.5000
   5.0000
mY*va=
   1.0000
   1.0000
vp=
   1.0000
   2.0000
vq=
   3.0000
   4.0000
mP=(vp vq)=
   1.0000   3.0000
   2.0000   4.0000
mQ=(vp' vq')'=
   1.0000
   2.0000
   3.0000
   4.0000
vp'vq=
   11.0000
(vp'*vq)*mP=
   11.0000   33.0000
   22.0000   44.0000

3.行列の要素を取り出す

// MySample3.cpp
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// using namespace arma;
// [[Rcpp::export]]

void fSample3()
{
  arma::mat mA, mB;
  double da;

  
  /*---行列の要素を参照---*/
  mA = { {1,2,3}, 
         {4,5,6},
         {7,8,9}};
  mA.print("mA=");
  ( mA.t() ).print("mA'=");
  
  da = mA(0, 0); // 行列mAの[0,0]要素をdaとおく．
  Rprintf("mA[0][0]= %5.3f\n", da);
  Rprintf("mA[0][1]= %5.3f\n", mA(0,1));
  Rprintf("mA[1][1]= %5.3f\n", mA(1,1));

  
  /*---スライス---*/
  ( mA(arma::span(0,0), arma::span(0,2)) ).print("mA[0][0:2]");
  ( mA(arma::span(0,2), arma::span(0,0)) ).print("mA[0:2][0]");
  ( mA(arma::span(0,1), arma::span(1,2)) ).print("mA[0:1][1:2]");
  ( mA(arma::span(1,2), arma::span(1,2)) ).print("mA[1:2][1:2]");
  
}  


/*** R
fSample3()
*/

###出力###
mA=
   1.0000   2.0000   3.0000
   4.0000   5.0000   6.0000
   7.0000   8.0000   9.0000
mA'=
   1.0000   4.0000   7.0000
   2.0000   5.0000   8.0000
   3.0000   6.0000   9.0000
mA[0][0]= 1.000
mA[0][1]= 2.000
mA[1][1]= 5.000
mA[0][0:2]
   1.0000   2.0000   3.0000
mA[0:2][0]
   1.0000
   4.0000
   7.0000
mA[0:1][1:2]
   2.0000   3.0000
   5.0000   6.0000
mA[1:2][1:2]
   5.0000   6.0000
   8.0000   9.0000

行列のインクリーメント，クロネッカー積

// MySample4.cpp
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// using namespace Rcpp
// using namespace arma;
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
void fSample4a()
{
 arma:: mat mA, mB, vc;
  
  mA = {{1,2},
        {3,4}};
  mA.print("mA=");
  mB = mA; // 代入．mAをいじってもmBは変わらない．

  mA++;  // mA = mA + 1; or mA += 1;と同値．各要素に1を加える．
  mA.print("mA=");
  mB.print("mB=");
  
  vc = {10, 20}; // 1 × 2 vector (縦ベクトル)
  vc = vc.t();  // 2 x 1 vector (横ベクトル)
  vc.print("vc=");
  
  mA--;  // mA = mA - 1; or mA -= 1;と同値．各要素から1引く．
  mA.print("mA=");
  
  (mA*mA).print("mA^2="); // 行列積
  (mA%mA).print("mA.^2="); // 要素積
  
  mB++;
  mB.print("mB=");
  ( mA*arma::inv(mB) ).print("mA*mB^{-1}="); // 行列商
  ( mA/mB).print("mA./mB="); // 各要素の割り算
  
  (mA*mB).print("mA*mB="); // 行列積
  (mA%mB).print("mA.*mB="); // 要素積
  ( arma::kron(mA, vc.t()) ).print("mA**vc'"); // クロネッカー積
}

// [[Rcpp::export]]
void fSample4b()
{
  double dx1, dx2, dx3, dx4;
  dx1 = 0; dx2= 1; dx3 = 2; dx4= 3;
  Rprintf("x1=%5.3f, x2=%5.3f, x3=%5.3f, x4=%5.3f\n", dx1, dx2, dx3, dx4);
  
  dx1 +=1; dx2 -=1; dx3 *=2, dx4 /=2;
  Rprintf("x1=%5.3f, x2=%5.3f, x3=%5.3f, x4=%5.3f\n", dx1, dx2, dx3, dx4);
}

/*** R
fSample4a() 
fSample4b() 
*/


###出力###
mA=
   1.0000   2.0000
   3.0000   4.0000
mA=
   2.0000   3.0000
   4.0000   5.0000
mB=
   1.0000   2.0000
   3.0000   4.0000
vc=
   10.0000
   20.0000
mA=
   1.0000   2.0000
   3.0000   4.0000
mA^2=
    7.0000   10.0000
   15.0000   22.0000
mA.^2=
    1.0000    4.0000
    9.0000   16.0000
mB=
   2.0000   3.0000
   4.0000   5.0000
mA*mB^{-1}=
   1.5000  -0.5000
   0.5000   0.5000
mA./mB=
   0.5000   0.6667
   0.7500   0.8000
mA*mB=
   10.0000   13.0000
   22.0000   29.0000
mA.*mB=
    2.0000    6.0000
   12.0000   20.0000
mA**vc'
   10.0000   20.0000   20.0000   40.0000
   30.0000   60.0000   40.0000   80.0000

> fSample4b() 
x1=0.000, x2=1.000, x3=2.000, x4=3.000
x1=1.000, x2=0.000, x3=4.000, x4=1.500

5. Loop

// MySample5.cpp
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// using namespace arma;

// [[Rcpp::export]]
void fSample5()
{
  int i;
  
  
  /*---for文---*/
  Rprintf("loop: increasing\n");
  for( i=0; i < 3; ++i)
    {
      Rprintf("i=%i\n",i);
    }
  
  Rprintf("loop: decreasing\n");
  for( i=3; i > 0; --i)
    {
      Rprintf("i=%i\n",i);
    }
  
  
  /*---while文---*/
  i = 0;
  Rprintf("loop: increasing\n");
  while(i < 3)
    {
      Rprintf("i=%i\n",i);
      ++i;
    }
  
  i = 3;
  Rprintf("loop: decreasing\n");
  while(i > 0)
    {
      Rprintf("i=%i\n",i);
      --i;
    }
  
  
  
  /*---do-while文*/
  i = 0;
  Rprintf("loop: increasing\n");
  do
    {
      Rprintf("i=%i\n",i);
      ++i;
    }while(i < 3);
  
  i = 3;
  Rprintf("loop: decreasing\n");
  do
    {
      Rprintf("i=%i\n",i);
      --i;
    }while(i > 0);
}

/*** R
fSample5() 
*/

loop: increasing
i=0
i=1
i=2
loop: decreasing
i=3
i=2
i=1
loop: increasing
i=0
i=1
i=2
loop: decreasing
i=3
i=2
i=1
loop: increasing
i=0
i=1
i=2
loop: decreasing
i=3
i=2
i=1
0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?