Merge pull request #15 from zwang2333/main

upgrade to 0.1.4
thebrisklab · Feb 5, 2025 · 6b53019 · 6b53019
2 parents 9c9dcda + f4c4e87
commit 6b53019
Show file tree

Hide file tree

Showing 10 changed files with 54 additions and 79 deletions.
diff --git a/.DS_Store b/.DS_Store
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: SparseICA
 Type: Package
 Title: Sparse Independent Component Analysis
-Version: 0.1.3
+Version: 0.1.4
 Date: 2025-01-15
 Author: Zihang Wang [aut, cre] (<https://orcid.org/0000-0001-9032-5412>),
   Irina Gaynanova [aut] (<https://orcid.org/0000-0002-4116-0268>),
@@ -21,7 +21,6 @@ LinkingTo: Rcpp, RcppArmadillo
 RoxygenNote: 7.3.1
 Encoding: UTF-8
 Imports:
-    RcppArmadillo (>= 14.2.2),
     Rcpp (>= 1.0.13),
     MASS (>= 7.3-58),
     irlba (>= 2.3.5),

diff --git a/NAMESPACE b/NAMESPACE
@@ -12,11 +12,10 @@ export(sparseICA)
 export(whitener)
 import(MASS)
 import(Rcpp)
-import(RcppArmadillo)
 import(ciftiTools)
 import(clue)
 import(irlba)
 import(parallel)
 importFrom("graphics", "abline")
 importFrom("stats", "coef", "cov", "lm", "rnorm", "runif", "sd")
-useDynLib(SparseICA, .registration = TRUE)
+useDynLib(SparseICA, .registration = TRUE)
diff --git a/R/SICAmain.R b/R/SICAmain.R
@@ -122,7 +122,7 @@ relax_and_split_ICA <- function(
   orth.method <- match.arg(orth.method)
   method <- match.arg(method)
 
-  if (verbose) cat("Centering/scaling data.\n")
+  if (verbose) message("Centering/scaling data.")
 
   xData_centered = scale(xData, center=TRUE, scale=FALSE)
   # Center and optionally scale data
@@ -141,7 +141,7 @@ relax_and_split_ICA <- function(
   if (d > p) stop("Number of components (n.comp) must not exceed number of columns in data.")
 
   # Whitening
-  if (verbose) cat("Whitening data.\n")
+  if (verbose) message("Whitening data.")
   whitener <- diag(d) # Default to identity
 
   if (lngca) {
@@ -164,7 +164,7 @@ relax_and_split_ICA <- function(
   }
   runs <- length(U.list)
 
-  if (verbose) cat("Starting Sparse ICA with", runs, "initializations.\n")
+  if (verbose) message("Starting Sparse ICA with ",runs," initialization(s).")
 
   ##############################################################################
   ############################ R implementation ################################
@@ -222,7 +222,7 @@ relax_and_split_ICA <- function(
 
     if(verbose){
       for (converge_i in 1:length(final_out_list$converge)) {
-        cat("Iteration ",converge_i,", tol=",final_out_list$converge[converge_i],".\n")}
+        message("Iteration ",converge_i,", tol = ",final_out_list$converge[converge_i],".")}
       }
 
     if (converge_plot) {
@@ -238,7 +238,7 @@ relax_and_split_ICA <- function(
     if(verbose){
       converge = final_Rcpp$converge[final_Rcpp$converge!=0]
       for (converge_i in 1:length(converge)) {
-        cat("Iteration ",converge_i,", tol=",converge[converge_i],".\n")
+        message("Iteration ",converge_i,", tol = ",converge[converge_i],".")
       }
     }
 
@@ -301,7 +301,7 @@ relax_and_split_ICA <- function(
 #' data(example_sim123)
 #'
 #' select_sparseICA = BIC_sparseICA(xData = example_sim123$xmat, n.comp = 3, 
-#'       method="C", BIC_plot = TRUE,verbose=FALSE, nu_list = seq(0.1,4,0.1))
+#'       method="C", BIC_plot = TRUE,verbose = TRUE, nu_list = seq(0.1,4,0.1))
 #'
 #' (my_nu = select_sparseICA$best_nu)
 #' }
@@ -322,14 +322,14 @@ BIC_sparseICA <- function(
   orth.method <- match.arg(orth.method)
   method <- match.arg(method)
 
-  if (verbose) cat("Centering and scaling data.\n")
+  if (verbose) message("Centering and scaling data.")
 
   # Center data
   xData_centered <- scale(xData, center = TRUE, scale = FALSE)
   v <- nrow(xData)
   t <- ncol(xData)
 
-  if (verbose) cat("Initializing reference Sparse ICA.\n")
+  if (verbose) message("Initializing reference Sparse ICA.")
 
   # Run reference Sparse ICA with very small `nu`
   ref_sparseICA <- relax_and_split_ICA(
@@ -345,7 +345,7 @@ BIC_sparseICA <- function(
   # Initialize storage for BIC values
   out_BIC <- numeric(length(nu_list))
 
-  if (verbose) cat("Starting BIC calculation over", length(nu_list), "nu values.\n")
+  if (verbose) message("Starting BIC calculation over ",length(nu_list)," nu values.")
 
   # Loop through `nu_list`
   for (i in seq_along(nu_list)) {
@@ -375,7 +375,7 @@ BIC_sparseICA <- function(
   # Select best `nu` based on minimum BIC
   best_nu <- nu_list[which.min(out_BIC)]
 
-  if (verbose) cat("The best nu selected by BIC is", best_nu, ".\n")
+  if (verbose) message("The best nu selected by BIC is ", best_nu, ".")
 
   # Optional BIC plot
   if (BIC_plot) {
@@ -390,7 +390,7 @@ BIC_sparseICA <- function(
     best_nu = best_nu
   )
 
-  if (verbose) print(Sys.time() - start.time)
+  if (verbose) message("BIC selection completed in ",round(difftime(Sys.time(),start.time,units = "secs"),2)," seconds.")
   return(result)
 }
 
@@ -447,10 +447,12 @@ BIC_sparseICA <- function(
 #' res_matched <- matchICA(my_sparseICA$estS,example_sim123$smat)
 #'
 #' # Visualize the estimated components
+#' oldpar <- par()$mfrow
 #' par(mfrow=c(1,3))
 #' image(matrix(res_matched[,1],33,33))
 #' image(matrix(res_matched[,2],33,33))
 #' image(matrix(res_matched[,3],33,33))
+#' par(mfrow=oldpar)
 #' }
 #' 
 #' @export
@@ -462,18 +464,18 @@ sparseICA <- function(
     verbose = TRUE, BIC_verbose = FALSE, converge_plot = FALSE, col.stand = TRUE,
     row.stand = FALSE, iter.stand = 5, positive_skewness = TRUE
 ) {
-  start.time <- Sys.time()
+  start.time1 <- Sys.time()
 
   if(nu == "BIC"){
-    if (verbose) cat("Selecting optimal nu using BIC-like criterion.\n")
+    if (verbose) message("Selecting the optimal nu using BIC-like criterion.")
     BIC_result <- BIC_sparseICA(
       xData = xData, n.comp = n.comp, nu_list = nu_list,
       whiten = whiten, lngca = lngca, orth.method = orth.method,
       method = method, use_irlba = use_irlba, eps = eps, maxit = maxit,
       verbose = BIC_verbose, col.stand = col.stand, row.stand = row.stand,
       iter.stand = iter.stand, BIC_plot = FALSE)
     nu <- BIC_result$best_nu
-    if (verbose) cat("The optimal nu is",nu,".\n")
+    if (verbose) message("The optimal nu is ",nu,".")
 
     sparseICA_result <- relax_and_split_ICA(
       xData = xData, n.comp = n.comp, U.list = U.list,
@@ -489,7 +491,7 @@ sparseICA <- function(
     sparseICA_result$BIC <- BIC_result$BIC
     sparseICA_result$nu_list <- BIC_result$nu_list
 
-    if (verbose) print(Sys.time() - start.time)
+    if (verbose) message("SparseICA is completed in ",round(difftime(Sys.time(),start.time1,units = "secs"),2)," seconds.")
 
     return(sparseICA_result)
   } else {
@@ -503,7 +505,7 @@ sparseICA <- function(
       positive_skewness = positive_skewness
     )
 
-    if (verbose) print(Sys.time() - start.time)
+    if (verbose) message("SparseICA is completed in ",round(difftime(Sys.time(),start.time1,units = "secs"),2)," seconds.")
 
     return(sparseICA_result)
   }

diff --git a/R/groupICA.R b/R/groupICA.R
@@ -60,33 +60,25 @@ create_group_list <- function(bids_path,pattern = "task-rest.*\\.dtseries\\.nii$
 #' @param npc An integer specifying the number of components to retain during subject-level PCA. Default is 85.
 #' @param iter_std An integer specifying the number of iterative standardization steps to apply to fMRI data. Default is 5.
 #' @param brainstructures A character vector specifying the brain structures to include in the analysis. Options are \code{"left"} (left cortex), \code{"right"} (right cortex), and/or \code{"subcortical"} (subcortex and cerebellum). Can also be \code{"all"} (obtain all three brain structures). Default is \code{c("left", "right")}.
+#' @param verbose A logical value indicating whether to print convergence information during execution. Default is \code{TRUE}.
 #'
 #' @return A numeric matrix containing the group-level principal components, with dimensions determined by the number of retained components (\code{n.comp}) and the concatenated data across all subjects.
 #'
 #' @details
 #' NOTE: This function requires the \code{ciftiTools} package to be installed, and set up the path to the Connectome Workbench folder by \code{ciftiTools.setOption()}. See the package \code{ciftiTools} documentation for more information.
 #'
-#' @examples
-#' # Example usage:
-#' # library(ciftiTools)
-#' # ciftiTools.setOption('wb_path', '/Applications/workbench')
-#' # Assuming `bids_dir` is the path to a BIDS dataset,
-#' # and `subject_list` is a named list of fMRI files:
-#' # groupPC <- gen_groupPC(bids_path = bids_dir, subj_list = subject_list, n.comp = 30, npc = 85)
-#' # print(dim(groupPC))
-#'
 #' @import ciftiTools
 #' @import irlba
 #' @import parallel
 #' @export
 gen_groupPC <- function(bids_path, subj_list, n.comp = 30, ncore=1,
-                        npc = 85, iter_std = 5, brainstructures = c("left", "right")) {
+                        npc = 85, iter_std = 5, brainstructures = c("left", "right"),verbose = TRUE) {
 
-  cat("# Start generating group PC of", length(subj_list), "subjects. #\n")
+  if (verbose) message("# Start generating group PC of", length(subj_list), "subjects. #")
 
   process_subject <- function(subject_name, subject_files) {
     nfile <- length(subject_files)
-    cat("##", subject_name, "has", nfile, "cortical surface fMRI data. ##\n")
+    if (verbose) message("##", subject_name, "has", nfile, "cortical surface fMRI data. ##")
 
     dat <- c()
     for (j in 1:nfile) {
@@ -107,14 +99,14 @@ gen_groupPC <- function(bids_path, subj_list, n.comp = 30, ncore=1,
       }
       dat <- cbind(dat, xii_mat)
     }
-    cat("## Total number of TRs:", dim(dat)[2], ".##\n")
+    if (verbose) message("## Total number of TRs:", dim(dat)[2], ".##")
 
     # Perform PCA
     subj_PCA <- prcomp_irlba(dat, npc)
     PC_subj <- subj_PCA$x
     dimnames(PC_subj) <- NULL
 
-    cat("## Retained", npc, "PCs. ##\n")
+    if (verbose) message("## Retained", npc, "PCs. ##")
     return(PC_subj)
   }
 
@@ -126,13 +118,13 @@ gen_groupPC <- function(bids_path, subj_list, n.comp = 30, ncore=1,
 
   # Combine all subject PCs
   groupPC <- do.call(cbind, groupPC_list)
-  cat("# Finish subject-level PCA! The concatenated matrix has dimension", dim(groupPC), ". #\n")
+  if (verbose) message("# Finish subject-level PCA! The concatenated matrix has dimension", dim(groupPC), ". #")
 
   # Perform group PCA
   temp <- whitener(X = groupPC, n.comp = n.comp, use_irlba = TRUE)
   groupPC <- temp$Z
 
-  cat("# Finish group PCA. #\n")
+  if (verbose) message("# Finish group PCA. #")
   return(groupPC)
 }
 
@@ -182,12 +174,6 @@ gen_groupPC <- function(bids_path, subj_list, n.comp = 30, ncore=1,
 #'   \item Executes Sparse ICA on the group-level PCs to estimate independent components.
 #' }
 #'
-#' @examples
-#' # Example usage:
-#' # Assuming `bids_dir` is the path to a BIDS dataset:
-#' # result <- group_sparseICA(bids_path = bids_dir, n.comp = 30, nu = "BIC")
-#' # str(result)
-#'
 #' @seealso \code{\link{create_group_list}}, \code{\link{gen_groupPC}}, \code{\link{BIC_sparseICA}}, \code{\link{sparseICA}}
 #' @export
 group_sparseICA <- function(bids_path, subj_list = NULL, nu = "BIC",
@@ -198,51 +184,51 @@ group_sparseICA <- function(bids_path, subj_list = NULL, nu = "BIC",
                             verbose = TRUE, BIC_verbose = FALSE, converge_plot = FALSE){
 
   if(verbose){
-    cat("##################################\n")
-    cat("##### Start group Sparse ICA #####\n")
-    cat("##################################\n\n")
-    cat("+++++ Step 1: Create subject list +++++\n")
+    message("##################################")
+    message("##### Start group Sparse ICA #####")
+    message("##################################")
+    message("+++++ Step 1: Create subject list +++++")
   }
 
   if(is.null(subj_list)){
     subj_list <- create_group_list(bids_path)
     if(verbose){
-      cat("# No input subject list, create using given path to BIDS. #\n")
-      cat("## Detect",length(subj_list),"subjects to be included. ##\n")
+      message("# No input subject list, create using given path to BIDS. #")
+      message("## Detect",length(subj_list),"subjects to be included. ##")
     }
   }else{
-    cat("# Use given subject list. #\n")
-    cat("## Detect",length(subj_list),"subjects to be included. ##\n")
+    message("# Use given subject list. #")
+    message("## Detect",length(subj_list),"subjects to be included. ##")
   }
 
   if(verbose){
-    cat("+++++ Step 2: Perform subject level PCA +++++\n")
+    message("+++++ Step 2: Perform subject level PCA +++++")
   }
 
   PC_group <- gen_groupPC(bids_path=bids_path,
                           subj_list = subj_list,
                           ncore = ncore,
                           n.comp = n.comp,
-                          npc = npc,iter_std = iter_std,brainstructures = brainstructures)
+                          npc = npc,iter_std = iter_std,brainstructures = brainstructures,verbose = verbose)
 
   if(verbose){
-    cat("+++++ Step 3: Select the tuning parameter +++++\n")
+    message("+++++ Step 3: Select the tuning parameter +++++")
   }
 
   if(nu == "BIC"){
-    cat("# Select nu by BIC. #\n")
+    if (verbose) message("# Select nu by BIC. #")
     nu_selection <- BIC_sparseICA(xData = PC_group, n.comp = n.comp, whiten = "none",method = method, 
                                   use_irlba = use_irlba,eps = eps,maxit = maxit,
                                   BIC_plot = BIC_plot,nu_list = nu_list,verbose=BIC_verbose)
     my_nu <- nu_selection$best_nu
-    cat("## The selected nu is",my_nu,". ##\n")
+    if (verbose) message("## The selected nu is",my_nu,". ##")
   }else{
     my_nu <- nu
-    cat("# Use given nu = ",my_nu,". #")
+    if (verbose) message("# Use given nu = ",my_nu,". #")
   }
 
   if(verbose){
-    cat("+++++ Step 4: Perform Sparse ICA on group PC +++++\n")
+    message("+++++ Step 4: Perform Sparse ICA on group PC +++++")
   }
 
   my_group_sparseICA <- sparseICA(xData = PC_group,
@@ -258,9 +244,9 @@ group_sparseICA <- function(bids_path, subj_list = NULL, nu = "BIC",
   my_group_sparseICA$nu_list <- nu_selection$nu_list
 
   if(verbose){
-    cat("###################################\n")
-    cat("##### Finish group Sparse ICA #####\n")
-    cat("###################################\n")
+    message("###################################")
+    message("##### Finish group Sparse ICA #####")
+    message("###################################")
   }
 
   return(my_group_sparseICA)

diff --git a/SparseICA.Rproj b/SparseICA.Rproj
@@ -15,3 +15,4 @@ LaTeX: pdfLaTeX
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
+PackageCheckArgs: --as-cran
diff --git a/man/BIC_sparseICA.Rd b/man/BIC_sparseICA.Rd
diff --git a/man/gen_groupPC.Rd b/man/gen_groupPC.Rd