From 0a92a2708a4e6ecc405a1a102ebdf0653b6b82c0 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Fri, 1 Mar 2024 13:52:26 +0100 Subject: [PATCH 1/4] Fixed C++ version of `reg()` (closes #28) The translation of the original R definition of `reg()` was done in d4e3ef1cfad96b7dc7f567929e6397badacd8bc7, but that version incorrectly assigns the slope of `new1` instead of the intercept to `beta01(e)`. This is fixed here --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 13 +++++++++++++ src/reg.cpp | 16 ++++++++++++++-- 3 files changed, 31 insertions(+), 2 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index e22b1d2..017fb2e 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -53,6 +53,10 @@ multiples_of <- function(x, divisor, subset_out = FALSE) { .Call(`_MADMMplasso_multiples_of`, x, divisor, subset_out) } +lm_arma <- function(R, Z) { + .Call(`_MADMMplasso_lm_arma`, R, Z) +} + reg <- function(r, Z) { .Call(`_MADMMplasso_reg`, r, Z) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index dd8e899..5900628 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -101,6 +101,18 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// lm_arma +arma::vec lm_arma(const arma::vec& R, const arma::mat& Z); +RcppExport SEXP _MADMMplasso_lm_arma(SEXP RSEXP, SEXP ZSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::vec& >::type R(RSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Z(ZSEXP); + rcpp_result_gen = Rcpp::wrap(lm_arma(R, Z)); + return rcpp_result_gen; +END_RCPP +} // reg Rcpp::List reg(const arma::mat r, const arma::mat Z); RcppExport SEXP _MADMMplasso_reg(SEXP rSEXP, SEXP ZSEXP) { @@ -143,6 +155,7 @@ static const R_CallMethodDef CallEntries[] = { {"_MADMMplasso_model_p", (DL_FUNC) &_MADMMplasso_model_p, 6}, {"_MADMMplasso_modulo", (DL_FUNC) &_MADMMplasso_modulo, 2}, {"_MADMMplasso_multiples_of", (DL_FUNC) &_MADMMplasso_multiples_of, 3}, + {"_MADMMplasso_lm_arma", (DL_FUNC) &_MADMMplasso_lm_arma, 2}, {"_MADMMplasso_reg", (DL_FUNC) &_MADMMplasso_reg, 2}, {"_MADMMplasso_scale_cpp", (DL_FUNC) &_MADMMplasso_scale_cpp, 2}, {"_MADMMplasso_sqrt_sum_squared_rows", (DL_FUNC) &_MADMMplasso_sqrt_sum_squared_rows, 1}, diff --git a/src/reg.cpp b/src/reg.cpp index a3cbe47..27efc6c 100644 --- a/src/reg.cpp +++ b/src/reg.cpp @@ -1,5 +1,17 @@ #include // [[Rcpp::depends(RcppArmadillo)]] + +// [[Rcpp::export]] +arma::vec lm_arma(const arma::vec &R, const arma::mat &Z) { + // Add a column of ones to Z + arma::mat Z_intercept = arma::join_rows(arma::ones(Z.n_rows), Z); + + // Solve the system of linear equations + arma::vec coefficients = arma::solve(Z_intercept, R); + + return coefficients; +} + // [[Rcpp::export]] Rcpp::List reg( const arma::mat r, @@ -10,9 +22,9 @@ Rcpp::List reg( arma::mat theta01(Z.n_cols, r.n_cols, arma::fill::zeros); for (arma::uword e = 0; e < r.n_cols; e++) { - arma::vec new1 = arma::solve(Z, r.col(e)); + arma::vec new1 = lm_arma(r.col(e), Z); beta01(e) = new1(0); - theta01.col(e) = new1.tail(new1.n_elem); + theta01.col(e) = new1.tail(new1.n_elem - 1); } Rcpp::List out = Rcpp::List::create( From 122566f25964d87cb6219a867dd8de8bb81272d7 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Fri, 1 Mar 2024 14:18:24 +0100 Subject: [PATCH 2/4] Added unit tests for #28 --- tests/testthat/test-reg.R | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 tests/testthat/test-reg.R diff --git a/tests/testthat/test-reg.R b/tests/testthat/test-reg.R new file mode 100644 index 0000000..a05be49 --- /dev/null +++ b/tests/testthat/test-reg.R @@ -0,0 +1,25 @@ +# Original function ============================================================ +reg_R <- function(r, Z) { + K <- ncol(Z) + N <- nrow(Z) + beta01 <- matrix(0, 1, ncol(r)) + theta01 <- matrix(0, ncol(Z), ncol(r)) + for (e in 1:ncol(r)) { + new1 <- lm(r[, e] ~ Z, singular.ok = TRUE) + beta01[e] <- matrix(new1$coefficients[1]) + theta01[, e] <- as.vector(new1$coefficients[-1]) + } + return(list(beta0 = beta01, theta0 = theta01)) +} + +# Testing ====================================================================== +reps <- 10L +n_obs <- rpois(reps, lambda = 10L) +n_vars <- sample(2:10, reps, replace = TRUE) +test_that("reg() produces the correct output", { + for (rp in seq_len(reps)) { + r <- matrix(rnorm(n_obs[rp] * n_vars[rp]), n_obs[rp], n_vars[rp]) + z <- as.matrix(sample(0:1, n_obs[rp], replace = TRUE)) + expect_equal(reg(r, z), reg_R(r, z)) + } +}) From 05d3bd76f0694c679871035ec2860bda7f74efca Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Fri, 1 Mar 2024 14:25:38 +0100 Subject: [PATCH 3/4] Replaced `1:ncol(r)` (#28) The usage is discouraged and problematic when `r` has one column. Thankfully caught by CodeFactor. --- tests/testthat/test-reg.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-reg.R b/tests/testthat/test-reg.R index a05be49..225a5e8 100644 --- a/tests/testthat/test-reg.R +++ b/tests/testthat/test-reg.R @@ -4,7 +4,7 @@ reg_R <- function(r, Z) { N <- nrow(Z) beta01 <- matrix(0, 1, ncol(r)) theta01 <- matrix(0, ncol(Z), ncol(r)) - for (e in 1:ncol(r)) { + for (e in seq_len(ncol(r))) { new1 <- lm(r[, e] ~ Z, singular.ok = TRUE) beta01[e] <- matrix(new1$coefficients[1]) theta01[, e] <- as.vector(new1$coefficients[-1]) From 0c77491e14906e4c8848b6ec3178ed39d2025ae5 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Fri, 1 Mar 2024 14:32:06 +0100 Subject: [PATCH 4/4] Removed lints in test for #28 --- tests/testthat/test-reg.R | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/testthat/test-reg.R b/tests/testthat/test-reg.R index 225a5e8..0d78536 100644 --- a/tests/testthat/test-reg.R +++ b/tests/testthat/test-reg.R @@ -1,7 +1,5 @@ # Original function ============================================================ reg_R <- function(r, Z) { - K <- ncol(Z) - N <- nrow(Z) beta01 <- matrix(0, 1, ncol(r)) theta01 <- matrix(0, ncol(Z), ncol(r)) for (e in seq_len(ncol(r))) { @@ -20,6 +18,6 @@ test_that("reg() produces the correct output", { for (rp in seq_len(reps)) { r <- matrix(rnorm(n_obs[rp] * n_vars[rp]), n_obs[rp], n_vars[rp]) z <- as.matrix(sample(0:1, n_obs[rp], replace = TRUE)) - expect_equal(reg(r, z), reg_R(r, z)) + expect_identical(reg(r, z), reg_R(r, z), tolerance = 1e-10) } })