裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

cxxfunction array の実装

2012年01月09日 | ブログラミング

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

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村