Rとのデータのやりとりをさらに復習する

```#include <RcppArmadilloExtensions/sample.h>
using namespace Rcpp ;

// [[Rcpp::export]]

int int_int(int n,int d) {
int ret = n+d;
return ret;
}

// [[Rcpp::export]]

double int_double(int n,int d) {
double nd = n*1.0;
double dd = d*1.0;
double ret = pow(nd,1/dd);
return ret;
}
// [[Rcpp::export]]

double double_double(double n,double d) {
double ret = pow(n,1/d);
return ret;
}

// [[Rcpp::export]]

IntegerVector int_intVec(int n,int d) {
RNGScope scope;
IntegerVector kv = seq_len(n);
IntegerVector r = RcppArmadillo::sample(kv, d, TRUE);
return r;
}

// [[Rcpp::export]]

IntegerMatrix int_intMat(int n,int d1,int d2) {
RNGScope scope;
IntegerVector kv = seq_len(n);
IntegerMatrix ret(d1,d2);
for(int i=0;i<d1;++i){
IntegerVector r = RcppArmadillo::sample(kv, d2, FALSE);
for(int j=0;j<d2;++j){
ret(i,j) = r(j);
}
}

return ret;
}

// [[Rcpp::export]]

List numVecNumMat_NumMat(NumericVector a,NumericMatrix b,double k) {
int n = b.nrow(), m = b.ncol();
NumericVector a2(a.size());
for(int i=0;i<n;++i){
a2[i] = a[i] * k;
}
arma::mat eigen_vectors(b.begin(),n,m,false);
arma::colvec eigen_values(a2.begin(),a2.size(),false);
return List::create(
Named("m") = eigen_vectors*(arma::diagmat(arma::exp(eigen_values)))*(arma::inv(eigen_vectors))
);

}
// [[Rcpp::export]]

List int_intList(int n,int d) {
IntegerMatrix x(n,d);
RNGScope scope;
IntegerVector kv = seq_len(n);

for(int i=0;i<n;++i){
IntegerVector r = RcppArmadillo::sample(kv, d, TRUE);
for(int j=0;j<d;++j){
x(i,j) = r[j];
}
}
return List::create(
Named("matrix") = x,
Named("sequence") = kv
) ;
}

/*** R

int_int(3,2)
int_double(3,2)
double_double(3,2)
int_intVec(6,10)
int_intMat(8,4,8)
A = matrix(runif(10000),100,100)
## Make it symmetric for stability of eigen decomposition
A = A + t(A)

A.eigen = eigen(A)

## Here is an R function to compute the matrix exponential

matexp.eigen = function(eigen.list,k){
return(eigen.list\$vectors%*%
diag(exp(eigen.list\$values*k))%*%solve(eigen.list\$vectors))
}
k <- 2
test.matexp1 = matexp.eigen(A.eigen,k)
test.matexp2 = numVecNumMat_NumMat(A.eigen\$values,A.eigen\$vectors,k)
range(unlist(test.matexp1)-unlist(test.matexp2))
plot(unlist(test.matexp1),unlist(test.matexp2))
abline(0,1,col=2)

int_intList(3,4)

*/
```
```> sourceCpp("InOutTemplate.cpp")

> int_int(3,2)
[1] 5

> int_double(3,2)
[1] 1.732051

> double_double(3,2)
[1] 1.732051

> int_intVec(6,10)
[1] 1 3 6 1 1 5 4 5 5 3

> int_intMat(8,4,8)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    5    8    7    3    1    4    2    6
[2,]    6    8    5    3    2    4    1    7
[3,]    5    4    8    6    7    2    3    1
[4,]    7    3    1    6    8    4    5    2

> A = matrix(runif(4),2,2)

> ## Make it symmetric for stability of eigen decomposition
> A = A + t(A)

> A.eigen = eigen(A)

> ## Here is an R function to compute the matrix exponential
>
> matexp.eigen = function(eigen.list,k){
+    return(eigen.list\$vectors%*%
+     diag( .... [TRUNCATED]

> test.matexp1 = matexp.eigen(A.eigen,1)

> test.matexp2 = numVecNumMat_NumMat(A.eigen\$values,A.eigen\$vectors,1)

> range(unlist(test.matexp1)-unlist(test.matexp2))
[1] 0 0

> int_intList(3,4)
\$matrix
[,1] [,2] [,3] [,4]
[1,]    1    3    3    2
[2,]    2    1    3    2
[3,]    1    3    2    1

\$sequence
[1] 1 2 3

>
```
• 指数行列の部分はこちらを参考にしました