save some memory

DESCRIPTION

@@ -45,5 +45,5 @@ Suggests:
 LazyData: yes
 LazyDataCompression: xz
 NeedsCompilation: no
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.2
 VignetteBuilder: knitr

NAMESPACE

@@ -41,6 +41,7 @@ importFrom(Matrix,rowSums)
 importFrom(Matrix,sparseMatrix)
 importFrom(Matrix,t)
 importFrom(Matrix,tcrossprod)
+importFrom(Rfast,med)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,aes_string)
 importFrom(ggplot2,annotate)
@@ -69,7 +70,6 @@ importFrom(stats,cor)
 importFrom(stats,cov2cor)
 importFrom(stats,dnorm)
 importFrom(stats,lm)
-importFrom(stats,median)
 importFrom(stats,optim)
 importFrom(stats,optimize)
 importFrom(stats,pnorm)

R/susie_rss.R

@@ -173,23 +173,6 @@ susie_rss = function (z, R, z_ld_weight = 0, prior_variance = 50,
     stop(paste0("The dimension of correlation matrix (", nrow(R)," by ",
                 ncol(R),") does not agree with expected (",length(z)," by ",
                 length(z),")"))
-  if (!is_symmetric_matrix(R))
-    stop("R is not a symmetric matrix")
-  if (!(is.double(R) & is.matrix(R)) & !inherits(R,"CsparseMatrix"))
-    stop("Input R must be a double-precision matrix, or a sparse matrix")
-
-  if (any(is.infinite(z)))
-    stop("z contains infinite value")
-
-  # Check for NAs in R.
-  if (any(is.na(R)))
-    stop("R matrix contains missing values")
-
-  # Replace NAs in z with zeros.
-  if (any(is.na(z))) {
-    warning("NA values in z-scores are replaced with 0")
-    z[is.na(z)] = 0
-  }
 
   # Modify R by z_ld_weight; this modification was designed to ensure
   # the column space of R contained z, but susie_suff_stat does not

R/susie_ss.R

@@ -141,7 +141,7 @@ susie_suff_stat = function (bhat, shat, R, n, var_y, XtX, Xty, yty,
                 ") does not agree with expected (",length(Xty)," by ",
                 length(Xty),")"))
   if (!is_symmetric_matrix(XtX))
-    stop("XtX is not a symmetric matrix")
+    stop("Input XtX or R is not a symmetric matrix")
 
   # MAF filter.
   if (!is.null(maf)) {
@@ -153,11 +153,17 @@ susie_suff_stat = function (bhat, shat, R, n, var_y, XtX, Xty, yty,
   }
 
   if (any(is.infinite(Xty)))
-    stop("Xty contains infinite values")
+    stop("Input Xty or zscores contains infinite values")
   if (!(is.double(XtX) & is.matrix(XtX)) & !inherits(XtX,"CsparseMatrix"))
-    stop("Input X must be a double-precision matrix, or a sparse matrix")
+    stop("Input XtX or R must be a double-precision matrix, or a sparse matrix")
   if (any(is.na(XtX)))
-    stop("XtX matrix contains NAs")
+    stop("Input XtX or R matrix contains NAs")
+  # Replace NAs in z with zeros.
+  if (any(is.na(Xty))) {
+    warning("NA values in Xty or z-scores are replaced with 0")
+    Xty[is.na(Xty)] = 0
+  }
+
 
   if (check_input) {
 
@@ -310,9 +316,15 @@ susie_suff_stat = function (bhat, shat, R, n, var_y, XtX, Xty, yty,
 
   # SuSiE CS and PIP.
   if (!is.null(coverage) && !is.null(min_abs_corr)) {
-    R = muffled_cov2cor(XtX)
-    s$sets = susie_get_cs(s,coverage = coverage,Xcorr = R,
-                          min_abs_corr = min_abs_corr)
+    if(any(!(diag(XtX) %in% c(0,1)))){
+      s$sets = susie_get_cs(s,coverage = coverage,Xcorr = muffled_cov2cor(XtX),
+                            min_abs_corr = min_abs_corr,
+                            check_symmetric = FALSE)
+    }else{
+      s$sets = susie_get_cs(s,coverage = coverage,Xcorr = XtX,
+                            min_abs_corr = min_abs_corr,
+                            check_symmetric = FALSE)
+    }
     s$pip = susie_get_pip(s,prune_by_cs = FALSE,prior_tol = prior_tol)
   }
 

R/susie_utils.R

@@ -246,11 +246,14 @@ susie_get_posterior_samples = function (susie_fit, num_samples) {
 #' @export
 #'
 susie_get_cs = function (res, X = NULL, Xcorr = NULL, coverage = 0.95,
-                         min_abs_corr = 0.5, dedup = TRUE, squared = FALSE) {
+                         min_abs_corr = 0.5, dedup = TRUE, squared = FALSE,
+                         check_symmetric = TRUE) {
   if (!is.null(X) && !is.null(Xcorr))
     stop("Only one of X or Xcorr should be specified")
-  if (!is.null(Xcorr) && !is_symmetric_matrix(Xcorr))
-    stop("Xcorr matrix must be symmetric")
+  if (check_symmetric){
+    if (!is.null(Xcorr) && !is_symmetric_matrix(Xcorr))
+      stop("Xcorr matrix must be symmetric")
+  }
   null_index = 0
   include_idx = rep(TRUE,nrow(res$alpha))
   if (!is.null(res$null_index)) null_index = res$null_index
@@ -442,7 +445,7 @@ n_in_CS = function(res, coverage = 0.9) {
 
 # Subsample and compute min, mean, median and max abs corr.
 #
-#' @importFrom stats median
+#' @importFrom Rfast med
 get_purity = function(pos, X, Xcorr, squared = FALSE, n = 100) {
   if (length(pos) == 1)
     c(1,1,1)
@@ -453,22 +456,26 @@ get_purity = function(pos, X, Xcorr, squared = FALSE, n = 100) {
 
       X_sub = X[,pos]
       if (length(pos) > n) {
-
         # Remove identical columns.
         pos_rm = sapply(1:ncol(X_sub),
                        function(i) all(abs(X_sub[,i] - mean(X_sub[,i])) <
                                        .Machine$double.eps^0.5))
         if (length(pos_rm))
           X_sub = X_sub[,-pos_rm]
       }
-      value = abs(muffled_corr(as.matrix(X_sub)))
-    } else
-      value = abs(Xcorr[pos,pos])
+      Rsub = muffled_corr(as.matrix(X_sub))
+      value = abs(Rsub[upper.tri(Rsub, diag = TRUE)])
+      rm(Rsub)
+    } else{
+      Xcorr = Xcorr[pos, pos]
+      value = abs(Xcorr[upper.tri(Xcorr, diag = TRUE)])
+      rm(Xcorr)
+    }
     if (squared)
       value = value^2
     return(c(min(value,na.rm = TRUE),
-             mean(value,na.rm = TRUE),
-             median(value,na.rm = TRUE)))
+             sum(value,na.rm = TRUE)/sum(!is.na(value)),
+             med(value,na.rm = TRUE)))
   }
 }