4

私は次の機能を持っています

fun <- cxxfunction(
    signature(x="numeric", y="numeric",N="interger", w="numeric", p="numeric"), 
    plugin="RcppArmadillo", 
    includes=c("#include <stdlib.h>", "#include <cmath>","#include <numeric>","#include <algorithm>","#include <vector>"), 
    body='
        using namespace Rcpp;
   RNGScope scope;

        NumericVector xa(x);
        NumericVector ya(y);
        int Na = as<int>(N);
        int n_xa = xa.size(), n_ya = ya.size();
        arma::mat wa = Rcpp::as<arma::mat>(w);
        Rcpp::NumericMatrix n(p); //one column to each x[i]/y
        NumericVector limite(n_xa);
        arma::mat z(n_ya, n_xa);
        arma::mat beta(n_ya, n_xa);
        arma::colvec one(xa.begin(), xa.size());
        arma::colvec betaM(ya.begin(), ya.size());
        arma::colvec betaP(ya.begin(), ya.size());
        arma::colvec Prob(ya.begin(), ya.size());
        arma::colvec betaC(ya.begin(), ya.size());
        z.col(0) = arma::zeros<arma::mat>(n_ya, 1);
        NumericVector nV(n_ya);
        NumericVector nT(n_ya);

        int i, j, randLR;

            for (i=0; i<n_xa; i++) {
                for (j=0; j<n_ya; j++) {
                   z(j,i) = 0;  
                }
                 randLR = rand() % n_ya;
         // randonly initialize z
                 z(randLR,i)=1;
                 limite(i)  = i;
                 one(i) = 1;
            }
    // the number of ways can be produced
            beta = wa * z;
            one = beta * one;


        int k, l,l2, pos;
        NumericMatrix prob_table(n_xa, Na);
        NumericMatrix class_table(n_xa, Na);

    int oldval;
            for(k=0; k<Na; k++){

               // limite = f(limite, n_xa);
       std::random_shuffle(limite.begin(),limite.end());
               for (l=0; l<n_xa; l++) {
                   l2=limite(l);     
        // trick, where the z.col came from 
                   betaM = one - wa * z.col(l2);
           // sum delta and normalize
           Prob = betaM;
                   std::transform(Prob.begin(), Prob.end(), Prob.begin(), std::bind2nd(std::plus<double>(),1));
                   std::transform(Prob.begin(), Prob.end(), Prob.begin(), std::bind2nd(std::divides<double>(),std::accumulate(Prob.begin(), Prob.end(), 0.0)));

                   nV = n(_, l2);
                   std::transform(nV.begin(), nV.end(), nV.begin(), std::bind2nd(std::divides<double>(),std::accumulate(nV.begin(), nV.end(), 0.0)));
           // multipling the the likelihoods
                   int nP;
                   for (nP=0; nP<betaP.size(); nP++) {
                            betaP(nP) = Prob(nP)*nV(nP);
            if(z(nP, l2)==1){
                oldval = nP;
            }

        }

                       // continuing control
                       if (std::accumulate(betaP.begin(), betaP.end(), 0.0) > 0) {
                            std::transform(betaP.begin(), betaP.end(), betaP.begin(), std::bind2nd(std::divides<double>(),std::accumulate(betaP.begin(), betaP.end(), 0.0)));
                       }
                       else {
                           // std::transform(betaP.begin(), betaP.end(), betaP.begin(), std::bind2nd(std::plus<double>(),1));
                       }
                        // a way to sample according a distribution, based on rogers
          double base = ::Rf_runif(0,1);

                       std::partial_sum (betaP.begin(), betaP.end(), betaC.begin());
                      int sz = betaC.size();
                      int idx;
                      for (idx=0; idx<sz; idx++) {
                          if(betaC[idx] > base) {
                              pos = idx;
                              break;
                          }
                      }


           if(pos!=oldval){
                        one = betaM - wa * z.col(l2);
               int nP;
                           for (nP=0; nP<betaP.size(); nP++) {
                    z(nP, l2) = 0;
                        }
                        z(pos, l2) = 1;
                        one = betaM + wa * z.col(l2);
        }       

                   prob_table(l2,k) = betaP(pos);
                   pos = pos + 1;
                   class_table(l2,k) = pos;
               }
}

//      betaP = z.col(l2);
        return Rcpp::List::create( Rcpp::Named("beta")=beta, Rcpp::Named("l2")=l2,Rcpp::Named("z")=z, Rcpp::Named("oldval")=oldval, Rcpp::Named("pos")=pos, Rcpp::Named("limite")=limite, Rcpp::Named("n")=nV, Rcpp::Named("betaM")=betaM, Rcpp::Named("betaP")=betaP, Rcpp::Named("Prob")=Prob,Rcpp::Named("prob_table")=prob_table, Rcpp::Named("class_table")=class_table );
        '
)

すべてのオペレーティング システム、64 ビット、8 GB RAM で、コードを次のようにロードします。

library(inline)
source("fun.R")

小さな例では、関数はすべてのシステムで完全に実行されますが、問題のサイズが大きくなると、windows と mac の増加が桁違いに大きくなるようです。

動作する小さな例:

varx <- 1:100
vary <- 1:200
w <- matrix(sample(0:1, 200^2, replace=TRUE), 200, 200)
wm <- matrix(abs(rnorm(200*100, 0.5, 0.5)), 200, 100)
system.time(conn <- fun(varx, vary, 5000, w, wm))

Linux でのみ動作する大きな例:

https://dl.dropboxusercontent.com/u/10712588/var1_teste.rda https://dl.dropboxusercontent.com/u/10712588/var2_teste.rda

load("var1_teste.rda")
load("var2_teste.rda")
varx <- 1:ncol(wl$wm)
vary <- 1:nrow(wl$wm)
system.time(conn <- fun(varx, vary, 5000, w, wl$wm))

ubuntu では 30 分、windows または mac では 2 日かかり、終了しません。運用システムに違いはありますか？

4

1 に答える 1

6

最初に小さな問題でテストすることを検討しましたか?

テストしたすべてのシステムでの実行可能性を確認するために、最初に既存の RcppArmadillo のサンプルの 1 つを再構築することを検討しましたか?

また、新しい Rcpp リリースでは、'Rcpp Attributes' 機能により、より簡単に行うことができます。検討：

R> cppFunction('arma::mat op(arma::colvec x) { return x*x.t(); }',
+              depends="RcppArmadillo")
R> op(1:3)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    4    6
[3,]    3    6    9
R>

ここでは、1 つの呼び出しと 1 つの (ラップされた) 行を使用して、Armadillo C++ 関数を作成し、列ベクトルの外積を生成します。

于 2013-09-22T23:05:10.963 に答える