[Rcpp-devel] Principal Component in Armadillo

[João Daniel Nunes Duarte]

Thanks for your explanation, Dirk! Indeed, I looked at Armadillo
documentation, but it was not clear what was the kind of return.

Regarding the linking error, I had created the package using
Rcpp.package.skeleton (because I didn't know I was going to use
RcppArmadillo), and then I made the necessary changes according to
documentation (
http://romainfrancois.blog.free.fr/index.php?post/2010/05/19/RcppArmadillo-0.2.1).
So I created the package again using "RcppArmadillo.package.skeleton" and
everything worked fine!

Actually, what "RcppArmadillo.package.skeleton" creates is slightly
different from what is described on documentation (maybe it is system
dependent):

- Documentation

# Makevars
PKG_LIBS = $(shell $(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()" )
$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

# Makevars.win
PKG_LIBS = $(shell Rscript.exe -e "Rcpp:::LdFlags()") $(LAPACK_LIBS)
$(BLAS_LIBS) $(FLIBS)

- Documentation

# Makevars
PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` $(LAPACK_LIBS)
$(BLAS_LIBS) $(FLIBS)

# Makevars.win
PKG_LIBS = $(shell $(R_HOME)/bin/Rscript.exe -e "Rcpp:::LdFlags()")
$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)


So, it seems that the R libraries was not being imported during the package
compilation; problem solved!

Regards,

João Daniel


2013/3/11 Dirk Eddelbuettel <edd at debian.org>

>
>
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130312/c52fc95d/attachment.html>