MM-58: refactoring should save a lot of memory (#54) [rinstall]

* MM-58: refactoring should save a lot of memory * MM-58: improved davies method with C * MM-58: not sure whether tests are working
lcrawlab · Apr 5, 2022 · 8e53b19 · 8e53b19
1 parent f5b9524
commit 8e53b19
Show file tree

Hide file tree

Showing 4 changed files with 46 additions and 46 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,13 +1,7 @@
-# Generated by using Rcpp::compileAttributes() -> do not edit by hand Generator
-# token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-MAPITCpp <- function(
-    X, Y, Z = NULL, C = NULL, variantIndices = NULL, testMethod = "normal", cores = 1L,
-    GeneticSimilarityMatrix = NULL
-) {
-    .Call(
-        "_mvMAPIT_MAPITCpp", PACKAGE = "mvMAPIT", X, Y, Z, C, variantIndices, testMethod,
-        cores, GeneticSimilarityMatrix
-    )
+MAPITCpp <- function(X, Y, Z = NULL, C = NULL, variantIndices = NULL, testMethod = "normal", cores = 1L, GeneticSimilarityMatrix = NULL) {
+    .Call('_mvMAPIT_MAPITCpp', PACKAGE = 'mvMAPIT', X, Y, Z, C, variantIndices, testMethod, cores, GeneticSimilarityMatrix)
 }
 
diff --git a/src/MAPIT.cpp b/src/MAPIT.cpp
@@ -129,12 +129,24 @@ Rcpp::List MAPITCpp(
       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);
+    }
+
     // 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);
-    arma::mat K = compute_k_matrix(GSM, x_k, p);
-    arma::mat G = compute_g_matrix(K, x_k);
+    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<microseconds>(end - start).count();
@@ -153,27 +165,16 @@ Rcpp::List MAPITCpp(
     }
     b.col(z + 1) = arma::trans(x_k);
 
-    arma::mat M = compute_projection_matrix(n, b);
+    M = compute_projection_matrix(n, b);
     b.reset();
-    arma::mat Kc = M * K * M;
-    K.reset();
-    arma::mat Gc = M * G * M;
-    G.reset();
-    arma::mat Cc;
-    std::vector<arma::mat> matrices;
+    K = M * K * M;
+    G = M * G * M;
 
     if (C.isNotNull()) {
-      Cc = M * Rcpp::as<arma::mat>(C.get()) * M;
-      matrices = {Gc, Kc, Cc, M};
-    } else {
-      matrices = {Gc, Kc, M};
+      Cc = M * Cc * M;
     }
-    Kc.reset();
-    Gc.reset();
-    Cc.reset();
 
     const arma::mat Yc = Y * M;
-    M.reset();
     end = steady_clock::now();
     execution_t(i, 1) = duration_cast<microseconds>(end - start).count();
 

diff --git a/src/mapit/davies.cpp b/src/mapit/davies.cpp
@@ -27,22 +27,16 @@ arma::vec davies_routine_vec(const arma::mat &S, const arma::mat &Sinv,
   arma::vec eigval;
   arma::mat eigvec;
 
-  if (num_variance_components == 3) {
-    // C is NULL
-    arma::vec q_sub = q.subvec(1, 2);
-    arma::mat S_sub = S.submat(1, 1, 2, 2);
-    arma::vec delta_null = arma::inv(S_sub) * q_sub;
-
-    arma::eig_sym(eigval, eigvec,
-                  delta_null(0) * matrices[1] + delta_null(1) * matrices[2]);
-  } else {
-    // C is not NULL
-    arma::vec delta_null = Sinv * q;
-
-    arma::eig_sym(eigval, eigvec,
-                  delta_null(0) * matrices[1] + delta_null(1) * matrices[2] +
-                      delta_null(2) * matrices[3]);
+  arma::vec q_sub = q.subvec(1, num_variance_components - 1);
+  arma::mat S_sub =
+      S.submat(1, 1, num_variance_components - 1, num_variance_components - 1);
+  arma::vec delta_null = arma::inv(S_sub) * q_sub;
+
+  arma::mat A = arma::zeros(arma::size(matrices[0]));
+  for (int i = 0; i < q_sub.n_elem; i++) {
+    A = A + delta_null(i) * matrices[i + 1];
   }
+  arma::eig_sym(eigval, eigvec, A);
   arma::mat EV_mat = compute_positive_ev_matrix(eigvec, eigval);
   evals = arma::eig_sym(EV_mat * compute_h_matrix(Sinv, matrices) * EV_mat);
   return evals;

diff --git a/tests/testthat/test-mvmapit-pariwise-executes-without-error.R b/tests/testthat/test-mvmapit-pariwise-executes-without-error.R
@@ -190,10 +190,13 @@ test_that(
         p <- 4
         n <- 10
         d <- 3
+
         pvalues <- matrix(
             c(
-                0.13395441, NA, 0, NA, NA, 0.009014612, NA, 0, NA, 0, NA, NA, 0,
-                NA, 0.01039137, NA, 0, NA, NA, 0, NA, 0, NA, 0, NA, NA, 0, NA
+                0.6977266, NA,0.4000627, NA, NA,0.3035032, NA,
+                0.4784944, NA,0.2936015, NA, NA,0.4665585, NA,
+                0.9274069, NA,0.1260672, NA, NA,0.1060084, NA,
+                0.4216691, NA,0.1615168, NA, NA,0.1526920, NA
             ),
             nrow = p, ncol = 7, byrow = TRUE
         )
@@ -230,8 +233,16 @@ test_that(
         d <- 1
         pvalues <- matrix(
             c(
-                0, 0.7422639, 0.7635732, 0, 0.1059719, 0.5984818, 0.7906796, 0, 0,
-                0.487107
+                0.7080633,
+                0.0000000,
+                0.0000000,
+                0.1740376,
+                0.2393295,
+                0.2031174,
+                0.0000000,
+                0.1372735,
+                0.1017353,
+                0.2053481
             ),
             nrow = p, ncol = 1, byrow = TRUE
         )