code and tests for 1.0.7 patch (#98)

zdk123 · Aug 6, 2019 · e4d7c0a · e4d7c0a
1 parent 0d48baf
commit e4d7c0a
Show file tree

Hide file tree

Showing 4 changed files with 47 additions and 32 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: SpiecEasi
 Title: Sparse Inverse Covariance for Ecological Statistical Inference
-Version: 1.0.6
+Version: 1.0.7
 Authors@R: c(
     person("Zachary", "Kurtz", role = c("aut", "cre"), email="zdkurtz@gmail.com"),
     person("Christian", "Mueller", role = "aut"),
@@ -16,7 +16,7 @@ Imports:
    methods,
    graphics,
    grDevices,
-   huge,
+   huge (>= 1.3.2),
    pulsar (>= 0.3.4),
    MASS,
    VGAM,

diff --git a/R/SparseICov.R b/R/SparseICov.R
@@ -51,20 +51,21 @@ sparseiCov <- function(data, method, npn=FALSE, verbose=FALSE, cov.output = TRUE
   method <- switch(method, glasso = "glasso", mb = "mb",
                    stop("Method not supported"))
 
-  if (is.null(args$lambda.min.ratio)) args$lambda.min.ratio <- 1e-3
+  if (is.null(args$lambda.min.ratio))
+    args$lambda.min.ratio <- 1e-3
 
-  if (method %in% c("glasso")) {
-    est <- do.call(huge::huge, c(args, list(x=data,
-                                            method=method,
-                                            verbose=verbose,
-                                            cov.output = cov.output)))
+  est <- do.call(huge::huge, c(args, list(x=data,
+                                          method=method,
+                                          verbose=verbose,
+                                          cov.output = cov.output)))
 
-  } else if (method %in% c('mb')) {
-    est <- do.call(utils::getFromNamespace('huge.mb', 'huge'),
-                   c(args, list(x=data, verbose=verbose)))
-    est$method <- 'mb'
-    est$data <- data
-    est$sym  <- ifelse(!is.null(args$sym), args$sym, 'or')
-  }
+  ## MB betas exported in huge>=1.3.2
+  # if (method %in% c('mb')) {
+  #   est <- do.call(utils::getFromNamespace('huge.mb', 'huge'),
+  #                  c(args, list(x=data, verbose=verbose)))
+  #   est$method <- 'mb'
+  #   est$data <- data
+  #   est$sym  <- ifelse(!is.null(args$sym), args$sym, 'or')
+  # }
   return(est)
 }
diff --git a/R/utilities.R b/R/utilities.R
@@ -66,10 +66,10 @@ getOptX.pulsar.refit <- function(est, getter='index') {
               est$select$stars$opt.index
                 },
         refit = {
-              est$refit$stars
+              Matrix::drop0(est$refit$stars)
         },
         merge = {
-              est$select$stars$merge[[getOptInd(est)]]
+              Matrix::drop0(est$select$stars$merge[[getOptInd(est)]])
         },
         stars = {
               est$select$stars$summary[getOptInd(est)]
@@ -79,19 +79,19 @@ getOptX.pulsar.refit <- function(est, getter='index') {
         },
         icov = {
           if (est$est$method == "glasso")
-            est$est$icov[[getOptInd(est)]]
+            Matrix::drop0(est$est$icov[[getOptInd(est)]])
           else
           stop("Run spiec-easi with method=\"glasso\"")
         },
         cov = {
           if (est$est$method == "glasso")
-            est$est$cov[[getOptInd(est)]]
+            Matrix::drop0(est$est$cov[[getOptInd(est)]])
           else
             stop("Run spiec-easi with method=\"glasso\"")
         },
         beta = {
           if (est$est$method == "mb")
-            est$est$beta[[getOptInd(est)]]
+            Matrix::drop0(est$est$beta[[getOptInd(est)]])
           else
             stop("Run spiec-easi with method=\"mb\"")
         }

diff --git a/tests/testthat/test_pulsar.R b/tests/testthat/test_pulsar.R
@@ -1,14 +1,14 @@
 context('setup')
 
-p <- 30
+p <- 20
 e <- p
 n <- 500
 set.seed(10010)
 g <- make_graph('erdos_renyi', p, e)
 S <- cov2cor(prec2cov(graph2prec(g)))
 X <- exp(rmvnorm(n, rep(0,p), S))
 
-pargs <- list(seed=10010, rep.num=15)
+pargs <- list(seed=10010, rep.num=10)
 
 context("SPIEC-EASI fit")
 lmx  <- .7
@@ -19,17 +19,26 @@ out <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx, lambda.min.rati
 
 ## StARs / B-StARS
 t1 <- system.time(
-out.stars <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx, lambda.min.ratio=lmr, nlambda=nlam, sel.criterion='stars', pulsar.select=TRUE, pulsar.params=pargs))
+out.stars <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx,
+  lambda.min.ratio=lmr, nlambda=nlam, sel.criterion='stars',
+  pulsar.select=TRUE, pulsar.params=pargs))
 t2 <- system.time(
-out.bstars <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx, lambda.min.ratio=lmr, nlambda=nlam, sel.criterion='bstars', pulsar.select=TRUE, pulsar.params=pargs))
+out.bstars <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx,
+  lambda.min.ratio=lmr, nlambda=nlam, sel.criterion='bstars',
+  pulsar.select=TRUE, pulsar.params=pargs))
 ## Batch Mode StARs / B-StARS
 options(batchtools.verbose=FALSE)
-bout.stars <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx, lambda.min.ratio=lmr, nlambda=nlam, sel.criterion='stars', pulsar.select='batch', pulsar.params=pargs)
-bout.bstars <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx, lambda.min.ratio=lmr, nlambda=nlam, sel.criterion='bstars', pulsar.select='batch', pulsar.params=pargs)
+bout.stars <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx,
+  lambda.min.ratio=lmr, nlambda=nlam, sel.criterion='stars',
+  pulsar.select='batch', pulsar.params=pargs)
+bout.bstars <- spiec.easi(X, method='mb', verbose=FALSE, lambda.max=lmx,
+  lambda.min.ratio=lmr, nlambda=nlam, sel.criterion='bstars',
+  pulsar.select='batch', pulsar.params=pargs)
 
 
 test_that("no pulsar has same output", {
-  expect_equal(as.matrix(out$est$path[[out.stars$select$stars$opt.index]]),
+  tmp <- out.stars$select$stars$opt.index
+  expect_equal(as.matrix(Matrix::drop0(out$est$path[[tmp]])),
                as.matrix(out.stars$refit$stars))
 })
 
@@ -67,12 +76,17 @@ test_that("Getter API throws errors if no pulsar selection", {
 })
 
 runtests <- function(out) {
-  expect_equal(getOptInd(out), (i<-out$select$stars$opt.index))
-  expect_equal(getOptNet(out), out$refit$stars)
-  expect_equal(getRefit(out), out$refit$stars)
+  expect_equal(getOptInd(out),
+              (i<-out$select$stars$opt.index))
+  expect_equal(getOptNet(out),
+               Matrix::drop0(out$refit$stars))
+  expect_equal(getRefit(out),
+               Matrix::drop0(out$refit$stars))
   expect_equal(getOptLambda(out), out$lambda[i])
-  expect_equal(getOptMerge(out), out$select$stars$merge[[i]])
-  expect_equal(getOptBeta(out), out$est$beta[[i]])
+  expect_equal(getOptMerge(out),
+               Matrix::drop0(out$select$stars$merge[[i]]))
+  expect_equal(getOptBeta(out),
+               Matrix::drop0(out$est$beta[[i]]))
 }
 
 test_that("Getter API, pulsar / stars ", {