アルマジロその2. 演算編
前回と同様に適当にはっつけるので, 気になったらコピペしてみて結果を確認してみてください.
おまじない
library (Rcpp) library (inline) library (RcppArmadillo) library (Matrix) # for sp_mat
行列積
# ベクトル, 行列の掛け算
src <- 'arma::rowvec a = arma::randu(1, 5);
arma::colvec b = arma::randu(5, 1);
arma::mat X = arma::randu(5, 5);
arma::vec c = arma::randu(5);
double q = as_scalar(a * b);
double p = as_scalar(a * X * b);
double l = as_scalar(a * arma::diagmat (X) * b);
double m = as_scalar(c.t() * arma::inv(arma::diagmat(X)) * c);
arma::mat Y = b * a; //外積だよ
// double v = arma::cross(c, c); //使えない....
// double j = arma::dot(c, c); //使えない....
return Rcpp::List::create(Rcpp::Named("a * b") = q,
Rcpp::Named("a * X * b") = p,
Rcpp::Named("a * dX * b") = l,
Rcpp::Named("ct * idX * c") = m,
Rcpp::Named("b * a") = Y);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
四則演算
# 四則演算的 (element_wise)
src <- 'arma::mat A = arma::randu(5,5);
arma::mat B = arma::randu(5,5);
arma::mat C = A + B;
arma::mat D = A - B;
arma::mat E = A / B;
arma::mat F = A % B; // 要素積
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);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
数学関数
# 数学関数 (element_wise)
src <- 'arma::mat A = arma::randu(3,3);
arma::mat B = arma::exp (A); // eのa_ij乗, exp2で2のa_ij乗, exp10で10のa_ij乗
arma::mat C = arma::trunc_exp(A);
arma::mat D = arma::log(A); //log2, log10, trunc_logがある
arma::mat E = arma::pow(A, 2); // 要素の二乗
arma::mat F = arma::sqrt(A);
arma::mat G = arma::sqrt(A);
arma::mat H = arma::floor(A); // ceilもある
arma::mat I = arma::round(A);
arma::mat J = arma::sign(A);
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);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
三角関数
# 三角関数(element_wise)
# cos family: cos, acos, cosh, acosh
# sin family: sin, asin, sinh, asinh
# tan family: tan, atan, tanh, atanh
src <- 'arma::mat A = arma::randu(3,3);
arma::mat B = arma::cos(A);
return wrap(B);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
クロネッカー積
# みんな大好きクロネッカー積～～～～
src <- 'arma::mat A = arma::mat (2,2,arma::fill::randu);
arma::mat B = arma::mat (2,2,arma::fill::randu);
arma::mat C = arma::kron (A, B);
return wrap (C);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
相関係数行列, 分散共分散行列
# 相関係数行行列, 分散共分散行列
# 引数をマトリックスにすると, 列ごとに変数が保持されていると考える
# ベクトルはひとつの変数の観測値がならんでいると考える
src <- 'arma::mat A = arma::randn(10,5);
arma::mat B = arma::randn(10,5);
arma::vec c = arma::randn(10);
arma::vec d = arma::randn(10);
arma::mat E = arma::cor(A, B);
arma::mat F = arma::cov(A, B); // 第三引数に, 0だとN-1, 1はN
double g = as_scalar(arma::cor(c, d));
return Rcpp::List::create(Named("E") = E,
Named("F") = F,
Named("g") = g);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
ノルム
# p次ノルム
# 1, 2, 3, ...., "-inf", "inf", "fro"
src <- 'arma::vec a = arma::randu(5);
arma::mat A = arma::randu(5,5);
double a2norm = arma::norm(a, 2);
double AFnorm = arma::norm(A, "fro");
return Rcpp::List::create(Rcpp::Named("a2") = a2norm,
Rcpp::Named("a * a") = arma::sqrt(a.t() * a),
Rcpp::Named("Af") = AFnorm);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
ランク
# ランク
src <- 'arma::mat A = arma::randu(5,5);
arma::uword r = arma::rank (A);
return wrap(r);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
rowMaxsなど
# rowMaxsなど
# 1が行, 0が列ごとでござる.
# 型は, 結果に合わせて宣言しなければならない
src <- 'arma::mat A = arma::randu(5,5);
arma::mat B = arma::randu(5,5);
arma::vec q = arma::randu(5);
arma::rowvec a = arma::max(A, 0);
arma::colvec b = arma::min(A, 1);
double d = max(max(A)); // Aの最大値
double qmax = q.max();
arma::mat C = max(A, B); // 要素ごとの最大
return Rcpp::List::create(Rcpp::Named("a") = a,
Rcpp::Named("b") = b,
Rcpp::Named("qmax") = qmax,
Rcpp::Named("d") = d,
Rcpp::Named("C") = C);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
sum, prod
# sum関係(prodも)
# 1が行, 0が列ごとでござる.
# 型は, 結果に合わせて宣言しなければならない
src <- 'arma::mat A = arma::randu(5,5);
arma::mat B = arma::randu(5,5);
arma::vec q = arma::randu(5);
double x = accu(A); // 全要素の和
arma::mat C = A + B; // AとBの要素和
arma::rowvec D = sum(A, 0); //Aの列和
arma::mat E = cumsum(A, 0); //Aの列の累積和
arma::rowvec F = prod(A, 0);//Aの列積
arma::vec p = cumsum(q);
arma::vec l = q + q;
double i = prod(q);
return Rcpp::List::create(Rcpp::Named("x") = x,
Rcpp::Named("C") = C,
Rcpp::Named("D") = D,
Rcpp::Named("E") = E,
Rcpp::Named("F") = F,
Rcpp::Named("p") = p,
Rcpp::Named("l") = l);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
基本統計量
# 基本統計量
src <- 'arma::mat A = arma::randu(5,5);
arma::vec q = arma::randu(5);
arma::rowvec C = arma::mean(A, 0); //Aの列平均
arma::colvec D = arma::median(A, 1); //Aの行の中央値
arma::rowvec E = arma::stddev(A, 1, 0); // 第二引数は不偏推定量かどうか, 0はN-1, 1はN
arma::rowvec F = arma::var(A, 0, 0);
double qvar = arma::var(q, 0);
return Rcpp::List::create(Rcpp::Named("C") = C,
Rcpp::Named("D") = D,
Rcpp::Named("E") = E,
Rcpp::Named("F") = F,
Rcpp::Named("qvar") = qvar);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun()
all, any
# bool
# allとany
# anyも同じ
src <- 'arma::vec V = arma::randu(10);
arma::mat X = arma::randu(5,5);
bool status1 = all(V); // ゼロ要素がないときtrue
bool status2 = arma::all(V>0.5);
bool status3 = arma::all(vectorise(X)>0.6);
arma::umat A = arma::all(X>0.7); // 列ごとに適用されるので, rowvec
return Rcpp::List::create(Named("s1") = status1,
Named("s2") = status2,
Named("s3") = status3,
Named("A") = A);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun()
sort
src <- 'arma::mat A = arma::randu(10, 10);
arma::mat B = arma::sort(A, "ascend", 0); // 列ごとに, 昇順. decendで降順
arma::vec C = arma::randu(10);
arma::uvec d = sort_index(C);
C(d).print(); //ソートした順番で出力できる
return Rcpp::List::create(Named("A") = A,
Named("B") = B,
Named("C") = C,
Named("d") = d);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun()
unique
# unique
src <- 'arma::mat A = arma::ones(3,3);
arma::mat Y = unique(A);
return wrap(Y);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun()
コンボルージョン
# コンボルージョン
# armaはリサイクル則が適用されるだと....!?
src <- 'arma::vec A = arma::randu(128) - 0.5;
arma::vec B = arma::randu(128) - 0.5;
arma::vec C = arma::conv(A,B);
return wrap(C);
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun()
フーリエ変換
# 高速フーリエ変換
# 複素のときはifft
# 第二引数は, 変換の長さ
src <- 'arma::vec X = arma::randu(100);
arma::cx_vec Y = arma::fft(X, 128);
Y.print();
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
# 二次元高速フーリエ変換
src <- 'arma::mat A = arma::randu(100, 100);
arma::cx_mat B = arma::fft2(A);
arma::cx_mat C = arma::fft2(A, 128, 128);
B.print();
'
fun <- cxxfunction (body = src, plugin = "RcppArmadillo")
fun ()
定数
src <- 'return wrap(arma::datum::pi); ' fun <- cxxfunction (body = src, plugin = "RcppArmadillo") fun () # predefined constants # # datum::pi π, the ratio of any circle's circumference to its diameter # datum::inf infinity # datum::nan “not a number” (NaN); caveat: NaN is not equal to anything, even itself # datum::e base of the natural logarithm # datum::sqrt2 square root of 2 # datum::eps the difference between 1 and the least value greater than 1 that is representable (type and machine dependant) # datum::log_min log of minimum non-zero value (type and machine dependant) # datum::log_max log of maximum value (type and machine dependant) # datum::euler Euler's constant, aka Euler-Mascheroni constant # datum::gratio golden ratio # datum::m_u atomic mass constant (in kg) # datum::N_A Avogadro constant # datum::k Boltzmann constant (in joules per kelvin) # datum::k_evk Boltzmann constant (in eV/K) # datum::a_0 Bohr radius (in meters) # datum::mu_B Bohr magneton # datum::Z_0 characteristic impedance of vacuum (in ohms) # datum::G_0 conductance quantum (in siemens) # datum::k_e Coulomb's constant (in meters per farad) # datum::eps_0 electric constant (in farads per meter) # datum::m_e electron mass (in kg) # datum::eV electron volt (in joules) # datum::ec elementary charge (in coulombs) # datum::F Faraday constant (in coulombs) # datum::alpha fine-structure constant # datum::alpha_inv inverse fine-structure constant # datum::K_J Josephson constant # datum::mu_0 magnetic constant (in henries per meter) # datum::phi_0 magnetic flux quantum (in webers) # datum::R molar gas constant (in joules per mole kelvin) # datum::G Newtonian constant of gravitation (in newton square meters per kilogram squared) # datum::h Planck constant (in joule seconds) # datum::h_bar Planck constant over 2 pi, aka reduced Planck constant (in joule seconds) # datum::m_p proton mass (in kg) # datum::R_inf Rydberg constant (in reciprocal meters) # datum::c_0 speed of light in vacuum (in meters per second) # datum::sigma Stefan-Boltzmann constant # datum::R_k von Klitzing constant (in ohms) # datum::b Wien wavelength displacement law constant #
