From 8e53b1945982326425e81dcee12e640cc87be70c Mon Sep 17 00:00:00 2001 From: jdstamp Date: Mon, 4 Apr 2022 20:33:24 -0400 Subject: [PATCH] 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 --- R/RcppExports.R | 14 +++----- src/MAPIT.cpp | 35 ++++++++++--------- src/mapit/davies.cpp | 24 +++++-------- ...-mvmapit-pariwise-executes-without-error.R | 19 +++++++--- 4 files changed, 46 insertions(+), 46 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 14a6a49e..e03abfa0 100644 --- 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 index b0b251be..f651cdec 100644 --- a/src/MAPIT.cpp +++ b/src/MAPIT.cpp @@ -129,12 +129,24 @@ Rcpp::List MAPITCpp( continue; } + std::vector 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(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(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 matrices; + K = M * K * M; + G = M * G * M; if (C.isNotNull()) { - Cc = M * Rcpp::as(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(end - start).count(); diff --git a/src/mapit/davies.cpp b/src/mapit/davies.cpp index 8bb49da8..06b8d803 100644 --- 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 index 9defce0c..82c8eecb 100644 --- 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 )