iramuteq.org Git - iramuteq/blob - Rscripts/anacor.R

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