2 #################################################################################
3 #http://www.mail-archive.com/rcpp-devel@lists.r-forge.r-project.org/msg01513.html
5 write.sparse <- function (m, to) {
6 ## Writes in a format that SVDLIBC can read
7 stopifnot(inherits(m, "dgCMatrix"))
8 fh <- file(to, open="w")
10 wl <- function(...) cat(..., "\n", file=fh)
13 wl(dim(m), length(m@x))
18 wl(nper[i]) ## Number of entries in this column
21 wl(m@i[globalCount], m@x[m@p[i]+j])
22 globalCount <- globalCount+1
27 my.svd <- function(x, nu, nv, libsvdc.path=NULL, sparse.path=NULL) {
30 outfile <- file.path(tempdir(),'sout')
31 cmd <- paste(libsvdc.path, '-o', outfile, '-d')
32 #rc <- system(paste("/usr/bin/svd -o /tmp/sout -d", nu, "/tmp/sparse.m"))
33 rc <- system(paste(cmd, nu, sparse.path))
35 stop("Couldn't run external svd code")
36 res1 <- paste(outfile,'-S', sep='')
37 d <- scan(res1, skip=1)
38 #FIXME : sometimes, libsvdc doesn't find solution with 2 dimensions, but does with 3
41 #rc <- system(paste("/usr/bin/svd -o /tmp/sout -d", nu, "/tmp/sparse.m"))
42 rc <- system(paste(cmd, nu, sparse.path))
43 d <- scan(res1, skip=1)
45 utfile <- paste(outfile, '-Ut', sep='')
46 ut <- matrix(scan(utfile,skip=1),nrow=nu,byrow=TRUE)
51 list(d=d, u=-t(ut), v=vt)
53 ###################################################################################
56 boostana<-function (tab, ndim = 2, libsvdc = FALSE, libsvdc.path=NULL)
58 #tab <- as.matrix(tab)
59 if (ndim > min(dim(tab)) - 1)
60 stop("Too many dimensions!")
61 name <- deparse(substitute(tab))
64 #tab <- reconstitute(tab, eps = eps)
68 #tab <- as.matrix(tab)
69 #prop <- as.vector(t(tab))/N
73 r <- ifelse(r == 0, 1, r)
74 c <- ifelse(c == 0, 1, c)
81 z <- as(z, "dgCMatrix")
82 tmpmat <- tempfile(pattern='sparse')
83 print('write sparse matrix')
84 write.sparse(z, tmpmat)
86 sv <- my.svd(z, qdim, qdim, libsvdc.path=libsvdc.path, sparse.path=tmpmat)
90 sv <- svd(z, nu = qdim, nv = qdim)
93 sigmavec <- (sv$d)[2:qdim]
94 x <- ((sv$u)/sqrt(r))[, -1]
96 x <- x * outer(rep(1, n), sigmavec)
97 dimlab <- paste("D", 1:ndim, sep = "")
98 colnames(x) <- dimlab# <- colnames(y) <- dimlab
99 rownames(x) <- rownames(tab)
100 result <- list(ndim = ndim, row.scores = x,
101 singular.values = sigmavec, eigen.values = sigmavec^2)
102 class(result) <- "boostanacor"