読者です 読者をやめる 読者になる 読者になる

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