array とはいっても，C 言語では単なる一次元ベクトルを三個以上の添え字でアクセスするだけ。
C 言語と R では，行と列の優先順序が正反対という点のみ。
「Rcpp::Arrayがないので」 では，本当の意味での array を操作していることにはなっていない。
つまり，添え字で array をアクセスするためには，以下のようにすべしということ。
# 必要なパッケージ
library(Rcpp)
library(inline)
# C++ ソース本体
src<-'
int Nx=as<int>(N1);
int Ny=as<int>(N2);
int Nz=as<int>(N3);
Rcpp::NumericVector a(Nx*Ny*Nz);
for (int k = 0; k < Nz; k++) {
for (int j = 0; j < Ny; j++) {
for (int i = 0; i < Nx; i++) {
a[i+Nx*j+Nx*Ny*k]= (i + 1)*(j + 1)*(k + 1)*1.0;
}
}
}
return a;
'
# R 内で、関数をコンパイル
fx.test1 <- cxxfunction(signature(N1 = "integer", N2 = "integer", N3 = "integer"),
include = c(""), src, plugin="Rcpp")
# 使う
library(rgl)
N1 <- 10
N2 <- 15
N3 <- 20
test.out <- fx.test1(N1, N2, N3)
test.out2 <- array(1, c(N1, N2, N3))
addresses <- which(test.out2 > 0, arr.ind = TRUE)
plot3d(addresses, col = gray((test.out - min(test.out)) /
(max(test.out) - min(test.out))))
A <- array(test.out, c(N1, N2, N3))
plot3d(slice.index(A, 1), slice.index(A, 2), slice.index(A, 3),
col = rainbow(256)[A / (max(A)) * 256 + 1], alpha = ifelse(A != 0, 1, 0))
以下に，基本的例を
> library(inline)
> src <- '
+ int n1 = as<int>(N1);
+ int n2 = as<int>(N2);
+ int n3 = as<int>(N3);
+ Rcpp::NumericVector x(n1*n2*n3);
+ int i, j, k;
+ for (k = 0; k < n3; k++) {
+ for (j = 0; j < n2; j++) {
+ for (i = 0; i < n1; i++) {
+ x[i+n1*j+n1*n2*k] = (i+1)*100+(j+1)*10+k+1;
+ }
+ }
+ }
+ return wrap(x);
+ '
>
> arry <- cxxfunction(signature(N1="integer", N2="integer", N3="integer"), src, plugin="Rcpp")
> a <- arry(2, 3, 4)
> array(a, dim=c(2, 3, 4))
, , 1
[,1] [,2] [,3]
[1,] 111 121 131
[2,] 211 221 231
, , 2
[,1] [,2] [,3]
[1,] 112 122 132
[2,] 212 222 232
, , 3
[,1] [,2] [,3]
[1,] 113 123 133
[2,] 213 223 233
, , 4
[,1] [,2] [,3]
[1,] 114 124 134
[2,] 214 224 234
配列の次元は，Rcpp の内部で .attr("dim") = Dimension(n1, n2, n3); によって定義することも出来る。
> library(inline)
> src <- '
+ int n1 = as<int>(N1);
+ int n2 = as<int>(N2);
+ int n3 = as<int>(N3);
+ Rcpp::NumericVector x(n1*n2*n3);
+ int i, j, k;
+ for (k = 0; k < n3; k++) {
+ for (j = 0; j < n2; j++) {
+ for (i = 0; i < n1; i++) {
+ x[i+n1*j+n1*n2*k] = (i+1)*100+(j+1)*10+k+1;
+ }
+ }
+ }
+ x.attr("dim") = Dimension(n1, n2, n3);
+ return wrap(x);
+ '
>
> arry <- cxxfunction(signature(N1="integer", N2="integer", N3="integer"), src, plugin="Rcpp")
> ( a <- arry(2, 3, 4) )
, , 1
[,1] [,2] [,3]
[1,] 111 121 131
[2,] 211 221 231
, , 2
[,1] [,2] [,3]
[1,] 112 122 132
[2,] 212 222 232
, , 3
[,1] [,2] [,3]
[1,] 113 123 133
[2,] 213 223 233
, , 4
[,1] [,2] [,3]
[1,] 114 124 134
[2,] 214 224 234
