www.iramuteq.org Git - iramuteq/blob - 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 }