add more unit tests, change ash code

R/sufficient_stats_methods.R

@@ -352,26 +352,15 @@ update_variance_components.ss <- function(data, params, model, ...) {
       est_sa2 <- est_sa2 * (0.1 / max(est_sa2))
     }
 
-    # Simplify precision matrix for ash
-    diagXtOmegaX_mrash <- colSums(data$X^2) / mom_result$sigma2 
-    XtOmega_mrash <- t(data$X) / mom_result$sigma2
-
-    # Call mr.ash with residuals and simplified precision matrix
-    mrash_output <- mr.ash.alpha.mccreight::mr.ash(
+    # Call mr.ash with residuals
+    mrash_output <- mr.ash.alpha::mr.ash(
       X             = data$X,
       y             = residuals,
       sa2           = est_sa2,
       intercept     = FALSE,
       standardize   = FALSE,
       sigma2        = mom_result$sigma2,
       update.sigma2 = FALSE,
-      diagXtOmegaX  = diagXtOmegaX_mrash,
-      XtOmega       = XtOmega_mrash,
-      V             = data$eigen_vectors,
-      tausq         = 0,
-      sum_Dsq       = sum(data$eigen_values),
-      Dsq           = data$eigen_values,
-      VtXt          = data$VtXt,
       max.iter      = 3000
     )
 

R/susie_constructors.R

@@ -45,7 +45,7 @@ individual_data_constructor <- function(X, y, L = min(10, ncol(X)),
 
   # Validate input X
   if (!(is.double(X) & is.matrix(X)) &
-      !inherits(X, "CsparseMatrix") &
+      !inherits(X, "sparseMatrix") &
       is.null(attr(X, "matrix.type"))) {
     stop("Input X must be a double-precision matrix, or a sparse matrix, or ",
          "a trend filtering matrix.")
@@ -238,7 +238,7 @@ sufficient_stats_constructor <- function(XtX, Xty, yty, n,
   }
 
   if (!(is.double(XtX) && is.matrix(XtX)) &&
-      !inherits(XtX, "CsparseMatrix")) {
+      !inherits(XtX, "sparseMatrix")) {
     stop("XtX must be a numeric dense or sparse matrix.")
   }
 
@@ -276,7 +276,7 @@ sufficient_stats_constructor <- function(XtX, Xty, yty, n,
   if (any(is.infinite(Xty))) {
     stop("Input Xty contains infinite values.")
   }
-  if (!(is.double(XtX) & is.matrix(XtX)) & !inherits(XtX, "CsparseMatrix")) {
+  if (!(is.double(XtX) & is.matrix(XtX)) & !inherits(XtX, "sparseMatrix")) {
     stop("Input XtX must be a double-precision matrix, or a sparse matrix.")
   }
   if (anyNA(XtX)) {
@@ -761,7 +761,7 @@ rss_lambda_constructor <- function(z, R, n = NULL,
   if (!isSymmetric(R)) {
     stop("R is not a symmetric matrix.")
   }
-  if (!(is.double(R) & is.matrix(R)) & !inherits(R, "CsparseMatrix")) {
+  if (!(is.double(R) & is.matrix(R)) & !inherits(R, "sparseMatrix")) {
     stop("Input R must be a double-precision matrix or a sparse matrix.")
   }
 

R/susie_utils.R

@@ -38,7 +38,7 @@ muffled_corr <- function(x) {
 #' @keywords internal
 muffled_cov2cor <- function(x) {
   withCallingHandlers(cov2cor(x), warning = function(w) {
-    if (grepl("had 0 or NA entries; non-finite result is doubtful", w$message)){
+    if (grepl("had.*(0|non-positive).*NA entries.*result.*(dubious|doubtful)", w$message)){
       invokeRestart("muffleWarning")}
     })
 }

inst/code/sparse_matrix_strategy.Rmd

@@ -121,7 +121,7 @@ p <- 10000
 
 ```{r}
 X.dense  <- create_sparsity_mat(0.99,n,p)
-X.sparse <- as(X.dense,"dgCMatrix")
+X.sparse <- as(X.dense,"sparseMatrix")
 X.tilde  <- susieR:::set_X_attributes(X.dense) #returns a scaled X if input is a dense matrix
 X <- susieR:::set_X_attributes(X.sparse) #return an unsacled sparse X if input is a sparse matrix 
                                      #but computes column means and standard deviations

tests/testthat/test_rss_utils.R

@@ -32,7 +32,7 @@ test_that("compute_suff_stat with sparse matrix input", {
   base_data <- generate_base_data(n = 10, p = 5, seed = 3)
 
   # Sparse version
-  X_sparse <- as(base_data$X, "CsparseMatrix")
+  X_sparse <- as(base_data$X, "sparseMatrix")
 
   out_dense <- compute_suff_stat(base_data$X, base_data$y, standardize = FALSE)
   out_sparse <- compute_suff_stat(X_sparse, base_data$y, standardize = FALSE)