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

#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]
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

> 

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