On 11 March 2013 at 22:42, Jo?o Daniel Nunes Duarte wrote:
| Hello,
|
| I'm trying to calculate principal component using 'princomp' function
from
| RcppArmadillo. Here's the cpp code:
|
| #include <RcppArmadillo.h>
|
| RcppExport SEXP pca(SEXP mats) {
| try {
|
| Rcpp::NumericMatrix matr(mats);
| int n = matr.nrow(), k = matr.ncol();
|
| arma::mat mat(matr.begin(), n, k, false);
| arma::colvec pca;
|
| pca = arma::princomp(mat);
|
| return Rcpp::wrap(mat);
|
| } catch(...) {
| ::Rf_error("c++ error");
| }
| }
|
| However, when I "R CMD check" the package, I get the following error:
|
| ** testing if installed package can be loaded
| Error in dyn.load(file, DLLpath = DLLpath, ...) :
| unable to load shared object
'/home/tecto/cpp/Rcpp/amora.Rcheck/amora/libs/
| amora.so':
| /home/tecto/cpp/Rcpp/amora.Rcheck/amora/libs/amora.so: undefined
symbol:
| dgesvd_
| Error: loading failed
| Execution halted
|
| I've read the Armadillo documentation for that function (http://
| arma.sourceforge.net/docs.html#princomp), however it was not clear to
me how to
| use it correctly.
|
| Does my code have any mistake?
Yes, several in fact. princomp() returns a matrix, not a vector. You were
trying to return mat when you probably meant to return pca.
As for your linking error: should not happen. On all systems,
RcppArmadillo
should get these BLAS functions from R by linking to R. Unless your R is
built in a weird way.
Here is what I get in the slightly corrected (and rewritten to use
sourceCpp() instead) example below:
R> sourceCpp("/tmp/joao.cpp")
R> M <- matrix(1:9, 3, 3)
R> localpca(M)
1.0000 4.0000 7.0000
2.0000 5.0000 8.0000
3.0000 6.0000 9.0000
0.5774 0.8165 0
0.5774 -0.4082 -0.7071
0.5774 -0.4082 0.7071
[,1] [,2] [,3]
[1,] 0.57735 0.816497 0.000000
[2,] 0.57735 -0.408248 -0.707107
[3,] 0.57735 -0.408248 0.707107
R>
The program follows:
-----------------------------------------------------------------------------
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::NumericMatrix localpca(Rcpp::NumericMatrix matr) {
int n = matr.nrow(), k = matr.ncol();
arma::mat mat(matr.begin(), n, k, false);
Rcpp::Rcout << mat << std::endl;
arma::mat pca = arma::princomp(mat);
Rcpp::Rcout << pca << std::endl;
return Rcpp::wrap(pca);
}
/*** R
M <- matrix(1:9, 3, 3)
localpca(M)
*/
-----------------------------------------------------------------------------
And for what it is worth, R returns the same (modulo ordering):
R> prcomp(M)
Standard deviations:
[1] 1.73205 0.00000 0.00000
Rotation:
PC1 PC2 PC3
[1,] 0.57735 0.000000 0.816497
[2,] 0.57735 -0.707107 -0.408248
[3,] 0.57735 0.707107 -0.408248
Hth, Dirk
|
| Cheers,
|
| Jo?o Daniel
|
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
--
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com