Hello, I had a problem this afternoon with an arma::mat constructor and the mat::insert_col method. Since I got the example from Dirk, I thought I would point out the problem I had and a solution. See
http://gallery.rcpp.org/articles/fast-linear-model-with-armadillo/ for background. Before I get to the example to reproduce the problem and the fix, let me say I'm posting from home as I can't post code from my work email; otherwise I'm fired. I do believe Rcpp::as should be used to marshall from NumericMatrix to arma::mat, not the arma::mat constructor that accesses raw memory, as Dirk uses. I read the Armadillo documentation this afternoon and do believe the problem with .insert_col is due to the handling of memory in the constructor mat(aux_mem*, n_rows, n_cols, copy_aux_mem = true, strict = true) see http://arma.sourceforge.net/docs.html Note that I've run similar examples on my Windows 7 box at work and on my OS X laptop (Mountain Lion) at home to confirm this is not os-dependent. Thanks, Dale Smith dtsm...@mindspring.com mat <- cbind(rnorm(50), rnorm(50), rnorm(50)) y <- rnorm(50) cppFunction(' List fastLm(NumericVector yr, NumericMatrix Xr) { int n = Xr.nrow(), k = Xr.ncol(); arma::mat X(Xr.begin(), n, k, false); arma::colvec y(yr.begin(), yr.size(), false); arma::colvec coef = arma::solve(X, y); arma::colvec resid = y - X*coef; arma::colvec ones = arma::ones<arma::colvec>(X.n_rows); X.insert_cols(1, ones); double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k)); arma::colvec stderrest = arma::sqrt( sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) ); return List::create(Named("coefficients") = coef, Named("stderr") = stderrest); }', depends = c("RcppArmadillo")) cppFunction(' List fastLmMod(NumericVector yr, NumericMatrix Xr) { int n = Xr.nrow(), k = Xr.ncol(); arma::mat X = Rcpp::as<arma::mat>(Xr); arma::colvec y(yr.begin(), yr.size(), false); arma::colvec ones = arma::ones<arma::colvec>(X.n_rows); X.insert_cols(0, ones); arma::colvec coef = arma::solve(X, y); arma::colvec resid = y - X*coef; double sig2 = arma::as_scalar(arma::trans(resid)*resid/(n-k)); arma::colvec stderrest = arma::sqrt( sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) ); return List::create(Named("coefficients") = coef, Named("stderr") = stderrest); }', depends = c("RcppArmadillo")) fastLm(y, mat) # error is "error: Mat::init(): mismatch between size of auxiliary memory and requested size" fastLmMod(y, mat) # coefficients are -0.15785722, -0.05668155, 0.15409366, -0.10079147 lm(y ~ mat) # confirm fastLmMod via lm function in R
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
