1 #################################################################################
2 #http://www.mail-archive.com/rcpp-devel@lists.r-forge.r-project.org/msg01513.html
4 write.sparse <- function (m, to) {
5 ## Writes in a format that SVDLIBC can read
6 stopifnot(inherits(m, "dgCMatrix"))
7 fh <- file(to, open="w")
9 wl <- function(...) cat(..., "\n", file=fh)
12 wl(dim(m), length(m@x))
17 wl(nper[i]) ## Number of entries in this column
20 wl(m@i[globalCount], m@x[m@p[i]+j])
21 globalCount <- globalCount+1
26 my.svd <- function(x, nu, nv, libsvdc.path=NULL, sparse.path=NULL) {
29 outfile <- file.path(tempdir(),'sout')
30 cmd <- paste(libsvdc.path, '-o', outfile, '-d')
31 #rc <- system(paste("/usr/bin/svd -o /tmp/sout -d", nu, "/tmp/sparse.m"))
32 rc <- system(paste(cmd, nu, sparse.path))
34 stop("Couldn't run external svd code")
35 res1 <- paste(outfile,'-S', sep='')
36 d <- scan(res1, skip=1)
37 #FIXME : sometimes, libsvdc doesn't find solution with 2 dimensions, but does with 3
40 #rc <- system(paste("/usr/bin/svd -o /tmp/sout -d", nu, "/tmp/sparse.m"))
41 rc <- system(paste(cmd, nu, sparse.path))
42 d <- scan(res1, skip=1)
44 utfile <- paste(outfile, '-Ut', sep='')
45 ut <- matrix(scan(utfile,skip=1),nrow=nu,byrow=TRUE)
50 list(d=d, u=-t(ut), v=vt)
52 ###################################################################################
55 boostana<-function (tab, ndim = 2, svd.method = 'svdR', libsvdc.path=NULL)
57 #tab <- as.matrix(tab)
58 if (ndim > min(dim(tab)) - 1)
59 stop("Too many dimensions!")
60 name <- deparse(substitute(tab))
63 #tab <- reconstitute(tab, eps = eps)
67 #tab <- as.matrix(tab)
68 #prop <- as.vector(t(tab))/N
72 r <- ifelse(r == 0, 1, r)
73 c <- ifelse(c == 0, 1, c)
78 if (svd.method == 'svdlibc') {
80 z <- as(z, "dgCMatrix")
81 tmpmat <- tempfile(pattern='sparse')
82 print('write sparse matrix')
83 write.sparse(z, tmpmat)
85 sv <- my.svd(z, qdim, qdim, libsvdc.path=libsvdc.path, sparse.path=tmpmat)
87 } else if (svd.method == 'svdR') {
89 sv <- svd(z, nu = qdim, nv = qdim)
91 } else if (svd.method == 'irlba') {
93 sv <- irlba(z, qdim, qdim)
96 sigmavec <- (sv$d)[2:qdim]
97 x <- ((sv$u)/sqrt(r))[, -1]
99 x <- x * outer(rep(1, n), sigmavec)
100 dimlab <- paste("D", 1:ndim, sep = "")
101 colnames(x) <- dimlab# <- colnames(y) <- dimlab
102 rownames(x) <- rownames(tab)
103 result <- list(ndim = ndim, row.scores = x,
104 singular.values = sigmavec, eigen.values = sigmavec^2)
105 class(result) <- "boostanacor"