font size...
[iramuteq] / Rscripts / anacor.R
1 #################################################################################
2 #http://www.mail-archive.com/rcpp-devel@lists.r-forge.r-project.org/msg01513.html
3
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")
8   
9   wl <- function(...) cat(..., "\n", file=fh)
10   
11   ## header
12   wl(dim(m), length(m@x))
13   
14   globalCount <- 1
15   nper <- diff(m@p)
16   for(i in 1:ncol(m)) {
17     wl(nper[i])  ## Number of entries in this column
18     if (nper[i]==0) next
19     for(j in 1:nper[i]) {
20       wl(m@i[globalCount], m@x[m@p[i]+j])
21       globalCount <- globalCount+1
22     }
23   }
24 }
25
26 my.svd <- function(x, nu, nv, libsvdc.path=NULL, sparse.path=NULL) {
27   print('my.svd')
28   stopifnot(nu==nv)
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))
33   if (rc != 0)
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 
38   if (length(d)==1) {
39       nu <- nu + 1
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)
43   }
44   utfile <- paste(outfile, '-Ut', sep='')
45   ut <- matrix(scan(utfile,skip=1),nrow=nu,byrow=TRUE)
46   if (nrow(ut)==3) {
47       ut <- ut[-3,]
48   }
49   vt <- NULL
50   list(d=d, u=-t(ut), v=vt)
51 }
52 ###################################################################################
53
54 #from anacor package
55 boostana<-function (tab, ndim = 2, svd.method = 'svdR', libsvdc.path=NULL) 
56 {
57     #tab <- as.matrix(tab)
58     if (ndim > min(dim(tab)) - 1) 
59         stop("Too many dimensions!")
60     name <- deparse(substitute(tab))
61     if (any(is.na(tab))) 
62         print('YA NA')        
63         #tab <- reconstitute(tab, eps = eps)
64     n <- dim(tab)[1]
65     m <- dim(tab)[2]
66     N <- sum(tab)
67     #tab <- as.matrix(tab)
68     #prop <- as.vector(t(tab))/N
69     r <- rowSums(tab)
70     c <- colSums(tab)
71     qdim <- ndim + 1
72     r <- ifelse(r == 0, 1, r)
73     c <- ifelse(c == 0, 1, c)
74     print('make z')
75     z1 <- t(tab)/sqrt(c)
76     z2 <- tab/sqrt(r)
77     z <- t(z1) * z2
78     if (svd.method == 'svdlibc') {
79         #START NEW SVD
80         z <- as(z, "dgCMatrix")
81         tmpmat <- tempfile(pattern='sparse')
82         print('write sparse matrix')
83         write.sparse(z, tmpmat)
84         print('do svd')
85         sv <- my.svd(z, qdim, qdim, libsvdc.path=libsvdc.path, sparse.path=tmpmat)
86         #END NEW SVD
87     } else if (svd.method == 'svdR') {
88         print('start R svd')
89         sv <- svd(z, nu = qdim, nv = qdim) 
90         print('end svd')
91     } else if (svd.method == 'irlba') {
92         print('irlba')
93         sv <- irlba(z, qdim, qdim)
94         print('end irlba')
95     }
96     sigmavec <- (sv$d)[2:qdim]
97         x <- ((sv$u)/sqrt(r))[, -1]
98     x <- x * sqrt(N)
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"
106     result
107 }