<div>First of all, great job with Rcpp - it's an awesome tool, and finally gave me a reason to go back to my C++ books and get beyond "Hello World"!</div>
<div> </div>
<div>I did have a question about RcppEigen, though: specifically, when I use RcppEigen through inline, is bounds checking disabled? If not, how could I go about disabling it? I did try adding the eigen_no_debug flag to my includes, but it wouldn't compile. (I know the preferred way might be to build a package, but inline is so much faster for quickly testing changes to my code, I hate to give it up right now.)</div>
<div> </div>
<div>As a separate but related question: is there agreement about good dynamic (not fixed-size) C++ analogues to R's array? Especially anaologues that are Rcpp-friendly. :) I really want to be able to quickly loop over and access elements using something like "array(i,j,k,l,m)".</div>
<div> </div>
<div>BTW, sorry if these are simple questions, but I'm learning C/C++ for this project as I go along.</div>
<div> </div>
<div>Thanks,</div>
<div>Andy</div>
<div> </div>
<div> </div>
<div>Background for the curious: </div>
<div> </div>
<div>The vast majority of my code is looping over elements of multidimensional numeric arrays (not fixed at compile time) and performing simple multiplication, addition, and subtraction on those elements. Basically, the time required to access elements of a numeric array is pretty important. A much much smaller set of my code requires some matrix algebra and linear algebra routines (matrix multiplication, cholesky decomp, etc.) on small vectors and matrices. Right now I use BLAS/LAPACK, but it's harder for me to read and experiment with, so I thought I'd give Eigen a try. </div>
<div> </div>
<div>Since RcppArmadillo was much slower than using STL vectors last time I tried it, I decided to benchmark Eigen before rewriting any code (see code below comparing STL vector, NumericMatrix, and MatrixXd). It seems a good bit slower, which I suspect might be due to bounds checking when I access elements of the matrix.</div>
<div> </div>
<div>Ultimately, I guess I could mix and match if I have to (use STL vectors for the simple operations, and other classes for the matrix algebra stuff), but for some reason that "feels" wrong..... but I'm no C++ person, so I'm probably wrong.</div>
<div> </div>
<div> </div>
<div># Here's a copy of the R syntax for the benchmark I mention above:</div>
<div> </div>
<div>library(RcppEigen)<br>library(inline)<br>library(rbenchmark)</div>
<div><br>mat1<-matrix(rbinom(200*200, 1, 0.5), ncol=200)</div>
<div> </div>
<div>header<-'<br>using namespace std;<br>using namespace Eigen;<br>'</div>
<div>RcppMatrixtest<-'<br>const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);</div>
<div>int p = X.nrow();<br>std::vector<double> V(p*p);<br>Eigen::MatrixXd EM(p,p);</div>
<div> </div>
<div>// Copy data into STL vector and Eigen dynamic matrix<br>for (int i=0; i<p; ++i){<br> for (int j=0; j<p; ++j){<br> EM(i,j) = V[i + j*p] = X(i,j);<br> }<br>}</div>
<div> </div>
<div>double tran = 0.0;</div>
<div> </div>
<div>for (int i=0; i<p; ++i){<br> for (int j=0; j<p; ++j){<br> for (int k=0; k<p; ++k){<br> if ((j != i) && (k != i) && (j != k)) {<br> tran += X(i,j)*X(j,k)*X(i,k);<br>
}<br> }<br> }<br>}</div>
<div> </div>
<div>return Rcpp::wrap(tran);<br>'</div>
<div>rcpp<-cxxfunction(<br> signature(<br> mat="numeric"<br> ),<br> body=RcppMatrixtest,<br> includes=header,<br> plugin="RcppEigen",<br> verbose=T<br> )</div>
<div><br>STLvectortest<-'<br>const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);</div>
<div>int p = X.nrow();<br>std::vector<double> V(p*p);<br>Eigen::MatrixXd EM(p,p);</div>
<div> </div>
<div>// Copy data into STL vector and Eigen dynamic matrix<br>for (int i=0; i<p; ++i){<br> for (int j=0; j<p; ++j){<br> EM(i,j) = V[i + j*p] = X(i,j);<br> }<br>}</div>
<div> </div>
<div>double tran = 0.0;</div>
<div> </div>
<div>for (int i=0; i<p; ++i){<br> for (int j=0; j<p; ++j){<br> for (int k=0; k<p; ++k){<br> if ((j != i) && (k != i) && (j != k)) {<br> tran += V[i + j*p]*V[j + k*p]*V[i + k*p];<br>
}<br> }<br> }<br>}</div>
<div> </div>
<div>return Rcpp::wrap(tran);<br>'</div>
<div>stl<-cxxfunction(<br> signature(<br> mat="numeric"<br> ),<br> body=STLvectortest,<br> includes=header,<br> plugin="RcppEigen",<br> verbose=T<br> )</div>
<div> </div>
<div>Reigentest<-'<br>const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);</div>
<div>int p = X.nrow();<br>std::vector<double> V(p*p);<br>Eigen::MatrixXd EM(p,p);</div>
<div> </div>
<div>// Copy data into STL vector and Eigen dynamic matrix<br>for (int i=0; i<p; ++i){<br> for (int j=0; j<p; ++j){<br> EM(i,j) = V[i + j*p] = X(i,j);<br> }<br>}</div>
<div> </div>
<div>double tran = 0.0;</div>
<div> </div>
<div>for (int i=0; i<p; ++i){<br> for (int j=0; j<p; ++j){<br> for (int k=0; k<p; ++k){<br> if ((j != i) && (k != i) && (j != k)) {<br> tran += EM(i,j)*EM(j,k)*EM(i,k);<br>
}<br> }<br> }<br>}</div>
<div> </div>
<div>return Rcpp::wrap(tran);<br>'</div>
<div>eig<-cxxfunction(<br> signature(<br> mat="numeric"<br> ),<br> body=Reigentest,<br> includes=header,<br> plugin="RcppEigen",<br> verbose=T<br> )</div>
<div> </div>
<div>Rcppbase<-'<br>const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);</div>
<div>int p = X.nrow();<br>std::vector<double> V(p*p);<br>Eigen::MatrixXd EM(p,p);</div>
<div>// Copy data into STL vector and Eigen dynamic matrix<br>for (int i=0; i<p; ++i){<br> for (int j=0; j<p; ++j){<br> EM(i,j) = V[i + j*p] = X(i,j);<br> }<br>}</div>
<div> </div>
<div>double tran = 0.0;</div>
<div> </div>
<div>return Rcpp::wrap(tran);<br>'</div>
<div>rcppbase<-cxxfunction(<br> signature(<br> mat="numeric"<br> ),<br> body=Rcppbase,<br> includes=header,<br> plugin="RcppEigen",<br> settings=settings,<br> verbose=T<br> )</div>
<div> </div>
<div># rcppbase is there to provide an estimate of how much time it takes to copy, create </div>
<div># necessary data structures, assuming it's not optimized away by the compiler....</div>
<div> </div>
<div>m1<-benchmark(<br> rcpp(mat1),<br> stl(mat1),<br> eig(mat1),<br> rcppbase(mat1),<br> replications=500)</div>
<div> </div>
<div>m1<br></div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>