アルマジロその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 #