アルマジロその1. 行列操作編
ここでは, RcppArmadilloにおける行列の初期化などを紹介します. 特に説明もせず, バンバン例を載せるので気になった所をコピペして結果を確認してください.
おまじない
library (Rcpp) library (inline) library (RcppArmadillo) library (Matrix) # for sp_mat （疎行列に必要）
行列の初期化
# 配列の初期化
# 複素行列(cx) 疎行列(sp)は, 必要になったらまた調べる
# rowvec, colvecをつくるときは, arma::randn(1, 5)とかにする
# arma::randn<arma::mat>(1,5)とかにすれば, そのまま使える
src <- 'arma::mat A = arma::mat(5,5); // アロケート 第三引数にfill値をいれてもよい
A.randn(); // .fill(0)で全て0になったりもする
arma::mat B = arma::mat(5,5, arma::fill::randu); // 一様乱数で初期化
arma::mat C = arma::zeros(5,5); // すべて0
arma::mat D = arma::ones(5,5); // すべて１
arma::mat E = arma::eye(5,5); // 単位行列
arma::mat F = arma::randn(5,5); // 標準正規乱数
arma::cube G = arma::randu(5,5,2); // 三次元配列と, 一様乱数
arma::vec H = arma::randu(5); // ベクトル（列）
arma::ivec I = arma::randi(5, arma::distr_param(-10, 10)); // 1から10までの整数から5つ, imat, icubeもある
arma::vec J = arma::linspace(10, 14, 10); // ivecはダメだよ, linspaceはseqと同じなので, seq(from, to, length.out)
arma::sp_mat K = arma::sp_mat(5,5); // RのdgCMatrixクラスで出力. Matrixが必要だよ
arma::mat Q = "1 2 3; 4 5 6;";
return Rcpp::List::create(Rcpp::Named("A") = A,
Rcpp::Named("B") = B,
Rcpp::Named("C") = C,
Rcpp::Named("D") = D,
Rcpp::Named("E") = E,
Rcpp::Named("F") = F,
Rcpp::Named("G") = G,
Rcpp::Named("H") = H,
Rcpp::Named("I") = I,
Rcpp::Named("J") = J,
Rcpp::Named("K") = K,
Rcpp::Named("Q") = Q);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
行列初期化２（疎行列）
src <- 'arma::sp_mat A = arma::sprandu(10, 10, 0.1); // 密度.1で一様乱数をいれた疎行列, sprandnもある
arma::sp_mat B = arma::spones(A); // Aの非ゼロ要素を１に変換
arma::sp_mat C = arma::speye(10, 10); // 疎行列の対角行列
return Rcpp::List::create(A, B, C);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
次元いじり
# 次元変え
src <- 'arma::mat A = Rcpp::rnorm(10, 0, 1);
arma::mat B;
arma::mat C;
arma::mat E;
A.reshape(5,2); // 5×2
B.copy_size(A); //Aと同じ大きさの行列つくる
C = repmat(A, 4, 5); //matlabと同じ
arma::mat D = A; //このコピーはあんましよくなさそう
arma::vec v = vectorise(A, 0); // Aの要素で一列のベクトルをつくる, 第二引数を1にすると行ベクトルになる
v.t().print();
E = arma::fliplr(A); //列を入れ替えたAをこぴ. flipudは行を入れ替える
A.swap_rows(0, 4); //Aの一行目と五行目を入れ替える. colsもある
//A.reset(); //<0,0>matrix
//A.set_size(2, 5);
//A.resize(10,10);
return wrap (0);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
行列の基本情報
# 行列の特性にアクセス
# 他にも, スパース行列の非0要素を把握するn_nonzeroがある
# ３次現配列で, ３次目を指定するときは .slices(k)や, .slices(t,u)など
src <- 'arma::mat A = Rcpp::rnorm(10, 0, 1); // 次元がなく, 要素だけ
A.reshape(5,2);
int n = A.n_cols, k = A.n_rows, nk = A.n_elem, nk2 = A.size();
A.clear(); //A のエレメントを消す
A.empty(); //Aが要素を持たないときtrue
bool a, b, c, d;
a = A.is_square();
b = A.is_vec();
c = A.is_rowvec();
d = A.is_colvec();
LogicalVector D = Rcpp::LogicalVector::create(a, b, c, d);
// return wrap(A.empty());
return Rcpp::List::create(Rcpp::Named("ncol") = n,
Rcpp::Named("nrow") = k,
Rcpp::Named("elem") = nk,
Rcpp::Named("nk2") = nk2,
Rcpp::Named("D") = D);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
対角成分
# 対角成分関係
src <- 'arma::mat A = Rcpp::rnorm(25, 0, 1); // 次元がなく, 要素だけ
A.reshape(5,5);
arma::vec a = A.diag(); // 列ベクトルになっているよ
arma::vec b = A.diag(1);
arma::vec c = A.diag(-2);
arma::vec d = arma::diagvec(A); // aと同じ
double e = trace (A);
arma::mat B = diagmat (A); // Aの対角成分だけの行列
arma::mat C = diagmat (a); // Bと同じになる
return Rcpp::List::create(Rcpp::Named("A") = A,
Rcpp::Named("a") = a,
Rcpp::Named("b") = b,
Rcpp::Named("c") = c,
Rcpp::Named("d") = d,
Rcpp::Named("e") = e,
Rcpp::Named("B") = B,
Rcpp::Named("C") = C,
Rcpp::Named("v") = v);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
転置とか
# 転置, 三角, 対称とか
src <- 'arma::mat A = Rcpp::rnorm(25, 0, 1); // 次元がなく, 要素だけ
A.reshape(5,5);
arma::mat At = A.t(); // もしくは, trans(A)
arma::mat U = trimatu(A); // 上三角を抽出
arma::mat L = trimatl(A);
arma::mat B = symmatu(A); // 上三角を使って対称行列をつくる
arma::mat C = symmatl(A);
return Rcpp::List::create(Rcpp::Named("At") = At,
Rcpp::Named("A") = A,
Rcpp::Named("U") = U,
Rcpp::Named("L") = L,
Rcpp::Named("B") = B,
Rcpp::Named("C") = C); // Rcpp::はなくてもよい？
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
交換
# 交換
# swap
src <- 'arma::mat A = arma::zeros(5,5);
arma::mat B = arma::ones(7, 6);
A.swap(B);
return Rcpp::List::create(A, B);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
行列の結合
# 行列の接合
# .insert_rowsもある
# .insert_slicesもある(for cuve)
src <- 'arma::mat A = arma::mat (5,10,arma::fill::randu);
arma::mat B = arma::mat (5,2,arma::fill::randu);
arma::mat C = arma::ones (5, 5);
arma::mat D, E;
A.insert_cols(2, B); // 3列目にBを入れる
B.insert_cols(1, 5); // ２列目に, 5×5の行列を入れる
C.insert_cols(5, A);
D = join_rows (C, A); // Cと同じ, join_rowsはRのcbindと同じ
E = join_horiz(A, B); // cbind(A, B), join_vertがrbind
return Rcpp::List::create(Rcpp::Named("A") = A,
Rcpp::Named("B") = B,
Rcpp::Named("C") = C,
Rcpp::Named("D") = D,
Rcpp::Named("E") = E);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
列ごとにhogehoge
# apply的な.
# .cols(0, 2)とかで, 1, 2, 3列にアクセス
# .col(0)で, 1列目だけにアクセス
# indicesを使うともっと便利
# each_row()もある
# cubeはslices
src <- 'arma::mat A = Rcpp::rnorm(25, 0, 1); // 次元がなく, 要素だけ
A.reshape(5,5);
arma::vec b = arma::linspace(10, 14, 5); // 10から14の列ベクトル
A.cols(0, 3).each_col() += b; // 列ごとにbを足す, 当然, .colsはなくて良い
arma::uvec indices(1); // インデックスの大きさ
indices(0) = 1;
A.each_col(indices) = b; //代入
arma::mat B = A.col(0);
return wrap(B);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
サブマトにアクセス
# 部分的操作
src <- 'arma::mat A = arma::zeros(5,5);
A.submat(0,1, 2,3) = arma::randn(3,3); // サブマトの左上と右下の指定. ここでは, A[1,2]からA[3,4]
A(arma::span(0,2), arma::span(1,3)) = arma::randn(3,3); // A[1:3, 2:4] in R
arma::mat B = A.submat(0, 1, 2, 3);
arma::mat C = A(arma::span(0,2), arma::span(1,3));
arma::mat D = A(0, 0, arma::size(3,3)); // A[1,1]を左上にして, 3×3の行列をとる
arma::mat X = arma::randn(5,5);
arma::uvec indices; // 参照インデックス用
// Xの要素のうち, 正要素を取り出す
//return wrap(X>0);
//return wrap(find(X>0));
indices = find(X>0);
X.elem(indices).print();
// Xの要素のうち, 0以上に10を足す
indices = find(X>0);
X.elem(indices) += 10;
X.print();
// 負の値をゼロにする
// くどいので, もっといい方法ありそう
indices = find(X<0);
arma::vec tem;
tem.copy_size(indices);
tem.fill(0);
X.elem(indices) = tem;
X.print();
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
# printデバック
src <- 'arma::mat A = arma::mat (2,2,arma::fill::randu);
A.print("A:"); // Aの出力
A.raw_print(); // 丸めずに出力
return wrap(0);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
Rcout
# Rcout
src <- 'arma::mat A = arma::mat (2,2,arma::fill::randu);
Rcout << "The value is " << A.row(1)(0) << std::endl;
return wrap(0);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
