Mm 87 rcpp progress bar (#106)

* MM-87 progress bar and interruption of c++ code * MM-87 check as cran runs * MM-87 logging that explains progress bar * MM-87 wrap long lines * MM-87 remaining long lines * MM-87 clang-format
lcrawlab · Aug 19, 2023 · 3547c63 · 3547c63
1 parent 008b4ac
commit 3547c63
Show file tree

Hide file tree

Showing 6 changed files with 127 additions and 115 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: mvMAPIT
 Title: Multivariate Genome Wide Marginal Epistasis Test
-Version: 2.0.1
+Version: 2.0.2
 URL: https://github.com/lcrawlab/mvMAPIT, https://lcrawlab.github.io/mvMAPIT/
 Authors@R: c(
     person("Julian", "Stamp", email = "julian_stamp@brown.edu",
@@ -62,6 +62,7 @@ LinkingTo:
     Rcpp,
     RcppArmadillo,
     RcppParallel,
+    RcppProgress,
     RcppSpdlog,
     testthat
 VignetteBuilder:

diff --git a/NEWS.md b/NEWS.md
@@ -2,11 +2,12 @@
 
 * Added Cauchy combination test `cauchy_combined` including vignette that compares combination methods
 * `simulate_traits` now returns genotype matrix with causal epistatic variants named according to the trait they affect
+* Added progress bar and possibility to interrupt C++ routine using `RcppProgress`
 
 # mvMAPIT 2.0.1 release
 
 * Fix LTO issues when submitting to CRAN. The testthat issue https://github.com/r-lib/testthat/issues/1230
-describes the solution chosen. Created this github gist to reproduce LTO errors: https://gist.github.com/jdstamp/056475683110aacdb1e6761872ab1e05.
+describes the solution chosen. Created this GitHub gist to reproduce LTO errors: https://gist.github.com/jdstamp/056475683110aacdb1e6761872ab1e05.
 
 # mvMAPIT 2.0.0.1 pre-release
 

diff --git a/R/MAPIT.R b/R/MAPIT.R
@@ -120,6 +120,7 @@ mvmapit <- function(
 
     variance_components_ind <- get_variance_components_ind(Y)
     if (test == "hybrid") {
+        log$info("Running normal C++ routine.")
         vc.mod <- MAPITCpp(X, Y, Z, C, variantIndex, "normal", cores = cores, NULL)  # Normal Z-Test
         pvals <- vc.mod$pvalues
         # row.names(pvals) <- rownames(X)
@@ -143,17 +144,19 @@ mvmapit <- function(
             threshold
         )
 
-        log$info("Running davies method on selected SNPs.")
+        log$info("Running davies C++ routine on selected SNPs.")
         vc.mod <- MAPITCpp(X, Y, Z, C, ind, "davies", cores = cores, NULL)
         davies.pvals <- mvmapit_pvalues(vc.mod, X, accuracy)
         pvals[, variance_components_ind][ind_matrix] <- davies.pvals[, variance_components_ind][ind_matrix]
     } else if (test == "normal") {
+        log$info("Running normal C++ routine.")
         vc.mod <- MAPITCpp(X, Y, Z, C, variantIndex, "normal", cores = cores, NULL)
         pvals <- vc.mod$pvalues
         pves <- vc.mod$PVE
         timings <- vc.mod$timings
     } else {
         ind <- ifelse(variantIndex, variantIndex, 1:nrow(X))
+        log$info("Running davies C++ routine.")
         vc.mod <- MAPITCpp(X, Y, Z, C, ind, "davies", cores = cores, NULL)
         pvals <- mvmapit_pvalues(vc.mod, X, accuracy)
         pvals <- set_covariance_components(variance_components_ind, pvals)

diff --git a/src/MAPIT.cpp b/src/MAPIT.cpp
@@ -11,12 +11,15 @@
 #include "mqs/mqs.h"
 
 using std::chrono::duration_cast;
+
 using std::chrono::milliseconds;
 using std::chrono::steady_clock;
 
-// [[Rcpp::depends(RcppArmadillo)]]
 // [[Rcpp::plugins(cpp11)]]
 // [[Rcpp::plugins(openmp)]]
+// [[Rcpp::depends(RcppArmadillo)]]
+// [[Rcpp::depends(RcppProgress)]]
+#include <progress.hpp>
 
 ////////////////////////////////////////////////////////////////////////////
 
@@ -116,144 +119,148 @@ Rcpp::List MAPITCpp(
 #ifdef _OPENMP
   omp_set_num_threads(cores);
 #endif
+  Progress prgrs(p, true);
 #pragma omp parallel for schedule(dynamic)
   for (i = 0; i < p; i++) {
+    if (!Progress::check_abort()) {
+      prgrs.increment();
 #ifdef WITH_LOGGER
-    logger->info("Variant {}/{}.", i + 1, p);
+      logger->info("Variant {}/{}.", i + 1, p);
 #endif
 
-    if (skip_variant(ind, i)) {
+      if (skip_variant(ind, i)) {
 #ifdef WITH_LOGGER
-      logger->info("Variant {}/{} not of interest. Skip to next.", i + 1, p);
+        logger->info("Variant {}/{} not of interest. Skip to next.", i + 1, p);
 #endif
-      continue;
-    }
+        continue;
+      }
 
-    std::vector<arma::mat> matrices;
-    matrices.resize(4);
-    arma::mat &G = matrices[0];
-    arma::mat &K = matrices[1];
-    arma::mat &M = matrices[2];
-    arma::mat &Cc = matrices[3];
-    if (C.isNotNull()) {
-      Cc = Rcpp::as<arma::mat>(C.get());
-    } else {
-      matrices.resize(3);
-    }
+      std::vector<arma::mat> matrices;
+      matrices.resize(4);
+      arma::mat &G = matrices[0];
+      arma::mat &K = matrices[1];
+      arma::mat &M = matrices[2];
+      arma::mat &Cc = matrices[3];
+      if (C.isNotNull()) {
+        Cc = Rcpp::as<arma::mat>(C.get());
+      } else {
+        matrices.resize(3);
+      }
 
-    // Compute K and G covariance matrices
-    auto start = steady_clock::now();
-    // Create the linear kernel
-    const arma::rowvec x_k = X(arma::span(i), arma::span::all);
-    K = compute_k_matrix(GSM, x_k, p);
-    G = compute_g_matrix(K, x_k);
+      // Compute K and G covariance matrices
+      auto start = steady_clock::now();
+      // Create the linear kernel
+      const arma::rowvec x_k = X(arma::span(i), arma::span::all);
+      K = compute_k_matrix(GSM, x_k, p);
+      G = compute_g_matrix(K, x_k);
 
-    auto end = steady_clock::now();
-    execution_t(i, 0) = duration_cast<milliseconds>(end - start).count();
+      auto end = steady_clock::now();
+      execution_t(i, 0) = duration_cast<milliseconds>(end - start).count();
 
 #ifdef WITH_LOGGER_FINE
-    logger->info("Dimensions of polygenic background: {} x {}.", K.n_cols,
-                 K.n_rows);
+      logger->info("Dimensions of polygenic background: {} x {}.", K.n_cols,
+                   K.n_rows);
 #endif
 
-    // Transform K and G using projection M
-    start = steady_clock::now();
-    arma::mat b = arma::zeros(n, z + 2);
-    b.col(0) = arma::ones<arma::vec>(n);
-    if (z > 0) {
-      b.cols(1, z) = Zz.t();
-    }
-    b.col(z + 1) = arma::trans(x_k);
+      // Transform K and G using projection M
+      start = steady_clock::now();
+      arma::mat b = arma::zeros(n, z + 2);
+      b.col(0) = arma::ones<arma::vec>(n);
+      if (z > 0) {
+        b.cols(1, z) = Zz.t();
+      }
+      b.col(z + 1) = arma::trans(x_k);
 
-    M = compute_projection_matrix(n, b);
-    K = project_matrix(K, b);
-    G = project_matrix(G, b);
+      M = compute_projection_matrix(n, b);
+      K = project_matrix(K, b);
+      G = project_matrix(G, b);
 
-    if (C.isNotNull()) {
-      Cc = project_matrix(Cc, b);
-    }
-    b.reset();
+      if (C.isNotNull()) {
+        Cc = project_matrix(Cc, b);
+      }
+      b.reset();
 
-    const arma::mat Yc = Y * M;
-    end = steady_clock::now();
-    execution_t(i, 1) = duration_cast<milliseconds>(end - start).count();
+      const arma::mat Yc = Y * M;
+      end = steady_clock::now();
+      execution_t(i, 1) = duration_cast<milliseconds>(end - start).count();
 
-    arma::mat q;
-    std::vector<arma::vec> phenotypes;
-    start = steady_clock::now();
-    if (combinatorial) {
-      phenotypes = matrix_to_vector_of_rows(Yc);
+      arma::mat q;
+      std::vector<arma::vec> phenotypes;
+      start = steady_clock::now();
+      if (combinatorial) {
+        phenotypes = matrix_to_vector_of_rows(Yc);
 
-    } else {
-      arma::vec yc = vectorise(Yc);
-      phenotypes = matrix_to_vector_of_rows(yc.as_row());
-    }
-    end = steady_clock::now();
-    execution_t(i, 2) = duration_cast<milliseconds>(end - start).count();
-
-    start = steady_clock::now();
-    q = compute_q_matrix(phenotypes, matrices);
-    end = steady_clock::now();
-    execution_t(i, 3) = duration_cast<milliseconds>(end - start).count();
-
-    start = steady_clock::now();
-    arma::mat S = compute_s_matrix(matrices);
-    end = steady_clock::now();
-    execution_t(i, 4) = duration_cast<milliseconds>(end - start).count();
-
-    // Compute delta and Sinv
-    const float det_S = arma::det(S);
-    if (det_S == 0) {
+      } else {
+        arma::vec yc = vectorise(Yc);
+        phenotypes = matrix_to_vector_of_rows(yc.as_row());
+      }
+      end = steady_clock::now();
+      execution_t(i, 2) = duration_cast<milliseconds>(end - start).count();
+
+      start = steady_clock::now();
+      q = compute_q_matrix(phenotypes, matrices);
+      end = steady_clock::now();
+      execution_t(i, 3) = duration_cast<milliseconds>(end - start).count();
+
+      start = steady_clock::now();
+      arma::mat S = compute_s_matrix(matrices);
+      end = steady_clock::now();
+      execution_t(i, 4) = duration_cast<milliseconds>(end - start).count();
+
+      // Compute delta and Sinv
+      const float det_S = arma::det(S);
+      if (det_S == 0) {
 #ifdef WITH_LOGGER
-      logger->warn("The determinant of the S matrix is {}.", det_S);
-      logger->info("Skip variant {}.", i + 1);
+        logger->warn("The determinant of the S matrix is {}.", det_S);
+        logger->info("Skip variant {}.", i + 1);
 #endif
-      continue;
-    }
-    arma::mat Sinv = arma::inv(S);
-    arma::mat delta = Sinv * q;
+        continue;
+      }
+      arma::mat Sinv = arma::inv(S);
+      arma::mat delta = Sinv * q;
 
-    // Save point estimates of the epistasis component
-    sigma_est.row(i) = delta.row(0);
+      // Save point estimates of the epistasis component
+      sigma_est.row(i) = delta.row(0);
 #ifdef WITH_LOGGER_FINE
-    logger->info("S-1({}):\n {}.", i + 1, matrix_to_string(Sinv));
-    logger->info("delta({}):\n {}.", i + 1, matrix_to_string(delta));
-    logger->info("q({}):\n {}.", i + 1, matrix_to_string(q));
+      logger->info("S-1({}):\n {}.", i + 1, matrix_to_string(Sinv));
+      logger->info("delta({}):\n {}.", i + 1, matrix_to_string(delta));
+      logger->info("q({}):\n {}.", i + 1, matrix_to_string(q));
 #endif
 
-    start = steady_clock::now();
-    if (testMethod == "normal") {
-      arma::vec var_delta =
-          compute_variance_delta(phenotypes, Sinv, delta, matrices);
-      // Save SE of the epistasis component
-      sigma_se.row(i) = arma::trans(sqrt(var_delta));
+      start = steady_clock::now();
+      if (testMethod == "normal") {
+        arma::vec var_delta =
+            compute_variance_delta(phenotypes, Sinv, delta, matrices);
+        // Save SE of the epistasis component
+        sigma_se.row(i) = arma::trans(sqrt(var_delta));
 #ifdef WITH_LOGGER_FINE
-      logger->info("var_delta({}):\n {}.", i + 1,
-                   vector_to_string(sqrt(var_delta)));
+        logger->info("var_delta({}):\n {}.", i + 1,
+                     vector_to_string(sqrt(var_delta)));
 #endif
-    } else if (testMethod == "davies") {
-      try {
-        Lambda.slice(i) = davies_routine(S, Sinv, q, matrices);
+      } else if (testMethod == "davies") {
+        try {
+          Lambda.slice(i) = davies_routine(S, Sinv, q, matrices);
 
 #ifdef WITH_LOGGER_FINE
-        logger->info("Lambda.slice({}). {}", i + 1,
-                     matrix_to_string(Lambda.slice(i)));
+          logger->info("Lambda.slice({}). {}", i + 1,
+                       matrix_to_string(Lambda.slice(i)));
 #endif
-      } catch (std::exception &e) {
+        } catch (std::exception &e) {
 #ifdef WITH_LOGGER
-        logger->error("Error: {}.", e.what());
-        logger->info("Skip davies method for variant {}.", i + 1);
+          logger->error("Error: {}.", e.what());
+          logger->info("Skip davies method for variant {}.", i + 1);
 #endif
-        continue;
+          continue;
+        }
       }
-    }
-    end = ::steady_clock::now();
-    execution_t(i, 5) = duration_cast<milliseconds>(end - start).count();
-    // Compute the PVE
-    pve.row(i) = compute_pve(delta, 0);
+      end = ::steady_clock::now();
+      execution_t(i, 5) = duration_cast<milliseconds>(end - start).count();
+      // Compute the PVE
+      pve.row(i) = compute_pve(delta, 0);
 #ifdef WITH_LOGGER_FINE
-    logger->info("PVE({}): \n{}.", i + 1, matrix_to_string(pve.row(i)));
+      logger->info("PVE({}): \n{}.", i + 1, matrix_to_string(pve.row(i)));
 #endif
+    }
   }
 #ifdef WITH_LOGGER
   logger->info("Elapsed time: {}", sw);

diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -12,29 +12,29 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
 // MAPITCpp
-Rcpp::List MAPITCpp(const arma::mat X, const arma::mat Y, Rcpp::Nullable<Rcpp::NumericMatrix> Z, Rcpp::Nullable<Rcpp::NumericMatrix> C, Rcpp::Nullable<Rcpp::NumericVector> variantIndices, std::string testMethod, int cores, Rcpp::Nullable<Rcpp::NumericMatrix> GeneticSimilarityMatrix);
+Rcpp::List MAPITCpp(const arma::mat X, const arma::mat Y, Rcpp::Nullable <Rcpp::NumericMatrix> Z, Rcpp::Nullable <Rcpp::NumericMatrix> C, Rcpp::Nullable <Rcpp::NumericVector> variantIndices, std::string testMethod, int cores, Rcpp::Nullable <Rcpp::NumericMatrix> GeneticSimilarityMatrix);
 RcppExport SEXP _mvMAPIT_MAPITCpp(SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP CSEXP, SEXP variantIndicesSEXP, SEXP testMethodSEXP, SEXP coresSEXP, SEXP GeneticSimilarityMatrixSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const arma::mat >::type X(XSEXP);
     Rcpp::traits::input_parameter< const arma::mat >::type Y(YSEXP);
-    Rcpp::traits::input_parameter< Rcpp::Nullable<Rcpp::NumericMatrix> >::type Z(ZSEXP);
-    Rcpp::traits::input_parameter< Rcpp::Nullable<Rcpp::NumericMatrix> >::type C(CSEXP);
-    Rcpp::traits::input_parameter< Rcpp::Nullable<Rcpp::NumericVector> >::type variantIndices(variantIndicesSEXP);
+    Rcpp::traits::input_parameter< Rcpp::Nullable <Rcpp::NumericMatrix> >::type Z(ZSEXP);
+    Rcpp::traits::input_parameter< Rcpp::Nullable <Rcpp::NumericMatrix> >::type C(CSEXP);
+    Rcpp::traits::input_parameter< Rcpp::Nullable <Rcpp::NumericVector> >::type variantIndices(variantIndicesSEXP);
     Rcpp::traits::input_parameter< std::string >::type testMethod(testMethodSEXP);
     Rcpp::traits::input_parameter< int >::type cores(coresSEXP);
-    Rcpp::traits::input_parameter< Rcpp::Nullable<Rcpp::NumericMatrix> >::type GeneticSimilarityMatrix(GeneticSimilarityMatrixSEXP);
+    Rcpp::traits::input_parameter< Rcpp::Nullable <Rcpp::NumericMatrix> >::type GeneticSimilarityMatrix(GeneticSimilarityMatrixSEXP);
     rcpp_result_gen = Rcpp::wrap(MAPITCpp(X, Y, Z, C, variantIndices, testMethod, cores, GeneticSimilarityMatrix));
     return rcpp_result_gen;
 END_RCPP
 }
 
-RcppExport SEXP run_testthat_tests(SEXP use_xml_sxp);
+RcppExport SEXP run_testthat_tests(void);
 
 static const R_CallMethodDef CallEntries[] = {
     {"_mvMAPIT_MAPITCpp", (DL_FUNC) &_mvMAPIT_MAPITCpp, 8},
-    {"run_testthat_tests", (DL_FUNC) &run_testthat_tests, 1},
+    {"run_testthat_tests", (DL_FUNC) &run_testthat_tests, 0},
     {NULL, NULL, 0}
 };
 

diff --git a/vignettes/study-wtccc-mice.Rmd b/vignettes/study-wtccc-mice.Rmd
@@ -216,7 +216,7 @@ In the first two columns, we list SNPs and their genetic location according to t
 
 11: Juan Carlos Souto, Laura Almasy, Montserrat Borrell, Francisco Blanco-Vaca, José Mateo, José Manuel Soria, Inma Coll, Rosa Felices, William Stone, Jordi Fontcuberta, and John Blangero. Genetic susceptibility to thrombosis and its relationship to physiological risk factors: The GAIT study. The American Journal of Human Genetics, 67(6):1452–1459, 2000. ISSN 0002-9297. doi: 10.1086/316903. URL https://www.sciencedirect.com/science/article/pii/S0002929707632145.
 
-12: Joachim Weischenfeldt, Inge Damgaard, David Bryder, Kim Theilgaard-Mönch, Lina A. Thoren, Finn Cilius Nielsen, Sten Eirik W. Jacobsen, Claus Nerlov, and Bo Torben Porse. NMD is essential for hematopoietic stem and progenitor cells and for eliminating by-products of programmed DNA rearrangements. Genes & Development, 22(10):1381–1396, 2008. ISSN 0890-9369, 1549-5477. doi: 10.1101/gad.468808. URL http://genesdev.cshlp.org/content/22/10/1381. Company: Cold Spring Harbor Laboratory Press Distributor: Cold Spring Harbor Laboratory Press Institution: Cold Spring Harbor Laboratory Press Label: Cold Spring Harbor Laboratory Press Publisher: Cold Spring Harbor Lab.
+12: Joachim Weischenfeldt, Inge Damgaard, David Bryder, Kim Theilgaard-Mönch, Lina A. Thoren, Finn Cilius Nielsen, Sten Eirik W. Jacobsen, Claus Nerlov, and Bo Torben Porse. NMD is essential for hematopoietic stem and progenitor cells and for eliminating by-products of programmed DNA rearrangements. Genes & Development, 22(10):1381–1396, 2008. ISSN 0890-9369, 1549-5477. doi: 10.1101/gad.468808. URL https://genesdev.cshlp.org/content/22/10/1381. Company: Cold Spring Harbor Laboratory Press Distributor: Cold Spring Harbor Laboratory Press Institution: Cold Spring Harbor Laboratory Press Label: Cold Spring Harbor Laboratory Press Publisher: Cold Spring Harbor Lab.
 
 13: Kevin Berendse, Maxim Boek, Marion Gijbels, Nicole N. Van der Wel, Femke C. Klouwer, Marius A. van den Bergh-Weerman, Abhijit Babaji Shinde, Rob Ofman, Bwee Tien Poll-The, Sander M.
 Houten, Myriam Baes, Ronald J. A. Wanders, and Hans R. Waterham. Liver disease predominates in a mouse model for mild human zellweger spectrum disorder. Biochimica et Biophysica Acta (BBA) - Molecular Basis of Disease, 1865(10):2774–2787, 2019. ISSN 0925-4439.