2013-07-09

返り値をリストで統一してみる

Rcpp

Rの関数の返り値はリストとするのが基本になっているので、Rcppを使うときもそう決めてしまうのも手かもしれません。以下のように…

// [[Rcpp::export]]

List hogeReturningList(****) {
	****
	return List::create(
		Named("hogeInt") = x,
		Named("hogeNumericVector") = y,
		Named("hogeMatrix") = z
		) ;
}

#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp ;

// [[Rcpp::export]]

List list_int(int n,int d) {
	int ret = n+d;
	return List::create(
		Named("") = ret
		) ;
}

// [[Rcpp::export]]

List list_double(int n,int d) {
	double nd = n*1.0;
	double dd = d*1.0;
	double ret = pow(nd,1/dd);
	return List::create(
		Named("") = ret
		) ;

}
// [[Rcpp::export]]

List list_double2(double n,double d) {
	double ret = pow(n,1/d);
	return List::create(
		Named("") = ret
		) ;
}


// [[Rcpp::export]]

List list_intVec(int n,int d) {
	RNGScope scope;
	IntegerVector kv = seq_len(n);
	IntegerVector ret = RcppArmadillo::sample(kv, d, TRUE);
	return List::create(
		Named("") = ret
		) ;

}

// [[Rcpp::export]]

List list_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 List::create(
		Named("") = ret
		) ;

}

// [[Rcpp::export]]


List list_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 list_matVec(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
		) ;
}


// [[Rcpp::export]]

/*** R

list_int(3,2)
list_double(3,2)
list_double2(3,2)
list_intVec(6,10)
list_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 = list_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)

list_matVec(3,4)

*/

> sourceCpp("AllListReturn.cpp")

> list_int(3,2)
[[1]]
[1] 5


> list_double(3,2)
[[1]]
[1] 1.732051


> list_double2(3,2)
[[1]]
[1] 1.732051


> list_intVec(6,10)
[[1]]
 [1] 5 1 1 1 1 1 1 4 5 2


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


> 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( .... [TRUNCATED] 

> k <- 2

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

> test.matexp2 = list_numMat(A.eigen$values,A.eigen$vectors,k)

> range(unlist(test.matexp1)-unlist(test.matexp2))
[1] -3.312838e+71  1.863472e+71

> plot(unlist(test.matexp1),unlist(test.matexp2))

> abline(0,1,col=2)

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

$sequence
[1] 1 2 3

>