I had a bash at this as I'm currently learning some C++ / Rcpp. I made a function where you pass in the row start index and row end index, and the desired matrix column (which is what I assumed you would want). <br><br>
Here is the code I came up with including some benchmarks I ran, interestingly, using a bigger matrix meant that the speedup from using RcppArmadillo was less, I guess this is to do with the question I asked earlier on the mailing list (<a href="http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2011-September/002876.html">http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2011-September/002876.html</a>) but I've not yet had time to fully understand the solutions.<br>
<br>Jeff<br><br>library(inline)<br>library(Rcpp)<br><br>src <- '<br> arma::mat mat = as<arma::mat>(sMat);<br> int rowFrom = as<int>(sRowFrom) - 1;<br> int rowTo = as<int>(sRowTo) - 1;<br>
int col = as<int>(sCol) - 1;<br><br> arma::colvec column = mat(arma::span(rowFrom, rowTo), arma::span(col, col));<br><br> return wrap(arma::var(column));<br> '<br><br>cppFunction <- cxxfunction(signature(sMat = "NumericMatrix", sRowFrom = "integer", sRowTo = "integer", sCol = "integer"), <br>
src, "RcppArmadillo")<br><br>rFunction <- function(mat, rowFrom, rowTo, col) {<br> return(var(mat[rowFrom:rowTo, col]))<br>}<br><br># test<br>mat <- matrix(rnorm(100), 10, 10)<br>all.equal(cppFunction(mat, 2, 5, 2), rFunction(mat, 2, 5, 2))<br>
<br># benchmark<br>library(rbenchmark)<br><br>benchmark(cppFunction(mat, 2, 5, 2), rFunction(mat, 2, 5, 2),<br> columns = c("test", "replications", "elapsed", "relative", "user.self", "sys.self"),<br>
order = "relative", replications = 1e5)<br><br> test replications elapsed relative user.self sys.self<br>1 cppFunction(mat, 2, 5, 2) 100000 0.82 1.000000 0.81 0<br>
2 rFunction(mat, 2, 5, 2) 100000 5.41 6.597561 5.41 0<br><br># big matrix benchmark<br>bigMat <- matrix(rnorm(10000), 100, 100)<br><br>benchmark(cppFunction(bigMat, 2, 5, 2), rFunction(bigMat, 2, 5, 2),<br>
columns = c("test", "replications", "elapsed", "relative", "user.self", "sys.self"),<br> order = "relative", replications = 1e5)<br><br>
test replications elapsed relative user.self sys.self<br>1 cppFunction(bigMat, 2, 5, 2) 100000 2.19 1.000000 2.19 0<br>2 rFunction(bigMat, 2, 5, 2) 100000 5.45 2.488584 5.43 0<br>
<br><br><div class="gmail_quote">On Wed, Sep 21, 2011 at 1:35 PM, Dirk Eddelbuettel <span dir="ltr"><<a href="mailto:edd@debian.org">edd@debian.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div class="im"><br>
On 20 September 2011 at 23:00, Noah Silverman wrote:<br>
| Nice suggestion, but it may not work for what I'm trying to do.<br>
|<br>
| I'm building up a matrix of values over a time series as part of a big loop.<br>
| At certain iterations, I need to calculate some summary statistics on a few<br>
| things. In R, it is trivial to subset any part of a matrix. Not sure to do<br>
| that in Rcpp.<br>
<br>
</div>RcppArmardillo makes Armadillo and Rcpp (and hence R) interoperable.<br>
<br>
Which means Zarrar has just solved your problem for you. I don't quite see<br>
where that 'may not work' as you claim. Care to expand?<br>
<br>
You also may want to look at things like the Rcpp::SubMatrix class. Here is a<br>
quick cut and paste from one of the unit test:<br>
<br>
NumericMatrix xx(4, 5);<br>
xx(0,0) = 3;<br>
xx(0,1) = 4;<br>
xx(0,2) = 5;<br>
xx(1,_) = xx(0,_);<br>
xx(_,3) = xx(_,2);<br>
SubMatrix<REALSXP> yy = xx( Range(0,2), Range(0,3) ) ;<br>
NumericMatrix res = yy ;<br>
return res;<br>
<br>
Supply a matrix and two ranges, get back a submatrix. So it's there. Also<br>
note the sugar syntax used here accessing entire rows and colums, just like<br>
in R.<br>
<br>
Dirk<br>
<div><div></div><div class="h5"><br>
|<br>
|<br>
| --<br>
| Noah Silverman<br>
| UCLA Department of Statistics<br>
| 8117 Math Sciences Building<br>
| Los Angeles, CA 90095<br>
|<br>
| On Sep 20, 2011, at 10:09 PM, Zarrar Shehzad wrote:<br>
|<br>
|<br>
| I am not sure if there are native functions in Rcpp but you could use<br>
| RcppArmadillo to solve your problem.<br>
|<br>
| So say Xs = x:<br>
|<br>
| // Convert from SEXP => Rcpp => Arma<br>
| Rcpp::NumericMatrix Xr(Xs);<br>
| arma::mat X(Xr.begin(), Xr.nrow(), Xr.ncol(), false);<br>
|<br>
| // Get subset of matrix and calculate variance<br>
| // (i.e., do var(x[3:10,5]))<br>
| arma::var( X.submat(3, 10, 5, 5) );<br>
|<br>
| The armadillo docs are really nice so you can see those for more info on<br>
| the var and submat functions.<br>
|<br>
| Cheers<br>
| Zarrar<br>
|<br>
| On Wed, Sep 21, 2011 at 12:16 AM, Noah Silverman <<a href="mailto:noahsilverman@ucla.edu">noahsilverman@ucla.edu</a>><br>
| wrote:<br>
|<br>
| Hello,<br>
|<br>
| I want to calculate the variance of a subset of a matrix column.<br>
|<br>
| For example, if I wanted the variance of items 3-10 in column 5.<br>
|<br>
| In R, this would be:<br>
|<br>
| x <- matrix(rnorm(100), nrow=10, ncol=10)<br>
| varx <- var(x[3:10,5])<br>
|<br>
| In Rcpp, I can construct a matrix object Rcpp::NumericMatrix x<br>
| The var function is available thanks to the great Sugar implementation<br>
|<br>
| But, how can I easily reference a subset of a column to calculate the<br>
| variance?<br>
|<br>
| Ideas?<br>
|<br>
|<br>
| --<br>
| Noah Silverman<br>
| UCLA Department of Statistics<br>
| 8117 Math Sciences Building #8208<br>
| Los Angeles, CA 90095<br>
|<br>
|<br>
| _______________________________________________<br>
| Rcpp-devel mailing list<br>
| <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
| <a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
|<br>
|<br>
|<br>
|<br>
|<br>
</div></div>| ----------------------------------------------------------------------<br>
<div class="im">| _______________________________________________<br>
| Rcpp-devel mailing list<br>
| <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
| <a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</div><font color="#888888">--<br>
New Rcpp master class for R and C++ integration is scheduled for<br>
San Francisco (Oct 8), more details / <a href="http://reg.info" target="_blank">reg.info</a> available at<br>
<a href="http://www.revolutionanalytics.com/products/training/public/rcpp-master-class.php" target="_blank">http://www.revolutionanalytics.com/products/training/public/rcpp-master-class.php</a><br>
</font><div><div></div><div class="h5">_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</div></div></blockquote></div><br>