Translated the hh-lambda loop

.gitignore

@@ -5,3 +5,4 @@
 src/*.o
 src/*.so
 src/*.dll
+aux/

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: MADMMplasso
 Title: Multi Variate Multi Response 'ADMM' with Interaction Effects
-Version: 0.0.0.9013
+Version: 0.0.0.9015
 Authors@R:
     c(
         person(

R/MADMMplasso.R

@@ -123,12 +123,6 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = 0.001, ma
     gg <- gg1
   }
 
-  lam_list <- list()
-  obj <- NULL
-  n_main_terms <- NULL
-  non_zero_theta <- NULL
-  my_obj <- list()
-
   my_W_hat <- generate_my_w(X = X, Z = Z)
 
   svd.w <- svd(my_W_hat)
@@ -202,14 +196,19 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = 0.001, ma
     doParallel::registerDoParallel(cl = cl)
     foreach::getDoParRegistered()
 
-    my_values <- foreach(i = 1:nlambda, .packages = "MADMMplasso", .combine = rbind) %dopar% {
+    my_values_matrix <- foreach(i = 1:nlambda, .packages = "MADMMplasso", .combine = rbind) %dopar% {
       admm_MADMMplasso(
         beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, my_W_hat, XtY,
         y, N, e.abs, e.rel, alpha, lam[i, ], alph, svd.w, tree, my_print,
         invmat, gg[i, ], legacy
       )
     }
     parallel::stopCluster(cl)
+
+    # Converting to list so hh_nlambda_loop_cpp can handle it
+    for (hh in seq_len(nrow(my_values_matrix))) {
+      my_values[[hh]] <- my_values_matrix[hh, ]
+    }
   } else if (!parallel && !pal) {
     my_values <- lapply(
       seq_len(nlambda),
@@ -221,87 +220,42 @@ MADMMplasso <- function(X, Z, y, alpha, my_lambda = NULL, lambda_min = 0.001, ma
         )
       }
     )
+  } else {
+    # This is triggered when parallel is FALSE and pal is 1
+    my_values <- list()
   }
 
-  hh <- 1
-  while (hh <= nlambda) {
-    lambda <- lam[hh, ]
-
-    start_time <- Sys.time()
-    if (pal) {
-      my_values <- admm_MADMMplasso(
-        beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, my_W_hat, XtY,
-        y, N, e.abs, e.rel, alpha, lambda, alph, svd.w, tree, my_print, invmat,
-        gg[hh, ], legacy
-      )
-
-      beta <- my_values$beta
-      theta <- my_values$theta
-      my_obj[[hh]] <- list(my_values$obj)
-      beta0 <- my_values$beta0
-      theta0 <- my_values$theta0 ### iteration
-      beta_hat <- my_values$beta_hat
-      y_hat <- my_values$y_hat
-    }
-    cost_time <- Sys.time() - start_time
-    print(cost_time)
-    if (parallel && !pal) {
-      beta <- my_values[hh, ]$beta
-      theta <- my_values[hh, ]$theta
-      my_obj[[hh]] <- list(my_values[hh, ]$obj)
-      beta0 <- my_values[hh, ]$beta0
-      theta0 <- my_values[hh, ]$theta0 ### iteration
-      beta_hat <- my_values[hh, ]$beta_hat
-      y_hat <- my_values[hh, ]$y_hat
-    } else if (!parallel && !pal) {
-      beta <- my_values[[hh]]$beta
-      theta <- my_values[[hh]]$theta
-      my_obj[[hh]] <- list(my_values[[hh]]$obj)
-      beta0 <- my_values[[hh]]$beta0
-      theta0 <- my_values[[hh]]$theta0 ### iteration
-      beta_hat <- my_values[[hh]]$beta_hat
-      y_hat <- my_values[[hh]]$y_hat
-    }
-
-    beta1 <- as(beta * (abs(beta) > tol), "sparseMatrix")
-    theta1 <- as.sparse3Darray(theta * (abs(theta) > tol))
-    beta_hat1 <- as(beta_hat * (abs(beta_hat) > tol), "sparseMatrix")
-
-    n_interaction_terms <- count_nonzero_a((theta1))
-
-    n_main_terms <- (c(n_main_terms, count_nonzero_a((beta1))))
-
-    obj1 <- (sum(as.vector((y - y_hat)^2))) / (D * N)
-    obj <- c(obj, obj1)
-
-    non_zero_theta <- (c(non_zero_theta, n_interaction_terms))
-    lam_list <- (c(lam_list, lambda))
-
-    BETA0[[hh]] <- beta0
-    THETA0[[hh]] <- theta0
-    BETA[[hh]] <- as(beta1, "sparseMatrix")
-    BETA_hat[[hh]] <- as(beta_hat1, "sparseMatrix")
-
-    Y_HAT[[hh]] <- y_hat
-    THETA[[hh]] <- as.sparse3Darray(theta1)
-
-    if (hh == 1) {
-      print(c(hh, (n_main_terms[hh]), non_zero_theta[hh], obj1))
-    } else {
-      print(c(hh, (n_main_terms[hh]), non_zero_theta[hh], obj[hh - 1], obj1))
-    }
-
-    hh <- hh + 1
-  } ### lambda
+  loop_output <- hh_nlambda_loop(
+    lam, nlambda, beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it,
+    my_W_hat, XtY, y, N, e.abs, e.rel, alpha, alph, svd.w, tree, my_print,
+    invmat, gg, tol, parallel, pal, BETA0, THETA0, BETA,
+    BETA_hat, Y_HAT, THETA, D, my_values, legacy
+  )
 
   remove(invmat)
-  remove(my_values)
   remove(my_W_hat)
 
-  obj[1] <- obj[2]
+  loop_output$obj[1] <- loop_output$obj[2]
 
-  pred <- data.frame(Lambda = lam, nzero = n_main_terms, nzero_inter = non_zero_theta, OBJ_main = obj)
-  out <- list(beta0 = BETA0, beta = BETA, BETA_hat = BETA_hat, theta0 = THETA0, theta = THETA, path = pred, Lambdas = lam, non_zero = n_main_terms, LOSS = obj, Y_HAT = Y_HAT, gg = gg)
+  pred <- data.frame(
+    Lambda = lam,
+    nzero = loop_output$n_main_terms,
+    nzero_inter = loop_output$non_zero_theta,
+    OBJ_main = loop_output$obj
+  )
+  out <- list(
+    beta0 = loop_output$BETA0,
+    beta = loop_output$BETA,
+    BETA_hat = loop_output$BETA_hat,
+    theta0 = loop_output$THETA0,
+    theta = loop_output$THETA,
+    path = pred,
+    Lambdas = lam,
+    non_zero = loop_output$n_main_terms,
+    LOSS = loop_output$obj,
+    Y_HAT = loop_output$Y_HAT,
+    gg = gg
+  )
   class(out) <- "MADMMplasso"
   # Return results
   return(out)

R/RcppExports.R

@@ -37,6 +37,22 @@ admm_MADMMplasso_cpp <- function(beta0, theta0, beta, beta_hat, theta, rho1, X,
     .Call(`_MADMMplasso_admm_MADMMplasso_cpp`, beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, W_hat, XtY, y, N, e_abs, e_rel, alpha, lambda, alph, svd_w, tree, invmat, gg, my_print)
 }
 
+count_nonzero_a_cpp <- function(x) {
+    .Call(`_MADMMplasso_count_nonzero_a_cpp`, x)
+}
+
+count_nonzero_a_sp_mat <- function(x) {
+    .Call(`_MADMMplasso_count_nonzero_a_sp_mat`, x)
+}
+
+count_nonzero_a_cube <- function(x) {
+    .Call(`_MADMMplasso_count_nonzero_a_cube`, x)
+}
+
+hh_nlambda_loop_cpp <- function(lam, nlambda, beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, my_W_hat, XtY, y, N, e_abs, e_rel, alpha, alph, svd_w, tree, my_print, invmat, gg, tol, parallel, pal, BETA0, THETA0, BETA, BETA_hat, Y_HAT, THETA, D, my_values) {
+    .Call(`_MADMMplasso_hh_nlambda_loop_cpp`, lam, nlambda, beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, my_W_hat, XtY, y, N, e_abs, e_rel, alpha, alph, svd_w, tree, my_print, invmat, gg, tol, parallel, pal, BETA0, THETA0, BETA, BETA_hat, Y_HAT, THETA, D, my_values)
+}
+
 model_intercept <- function(beta0, theta0, beta, theta, X, Z) {
     .Call(`_MADMMplasso_model_intercept`, beta0, theta0, beta, theta, X, Z)
 }

R/admm_MADMMplasso.R

@@ -41,8 +41,8 @@ admm_MADMMplasso <- function(beta0, theta0, beta, beta_hat, theta, rho1, X, Z, m
     return(out)
   }
   warning(
-    "Using legacy R code for MADMMplasso.",
-    "This functionality will be removed in a future release.",
+    "Using legacy R code for MADMMplasso. ",
+    "This functionality will be removed in a future release. ",
     "Please consider using legacy = FALSE instead."
   )
   TT <- tree

R/cv_MADMMplasso.R

@@ -25,7 +25,9 @@ cv_MADMMplasso <- function(fit, nfolds, X, Z, y, alpha = 0.5, lambda = fit$Lambd
   nfolds <- length(table(foldid))
 
   for (ii in 1:nfolds) {
-    print(c("fold,", ii))
+    if (my_print) {
+      print(c("fold,", ii))
+    }
     oo <- foldid == ii
 
     ggg[[ii]] <- MADMMplasso(X = X[!oo, , drop = FALSE], Z = Z[!oo, , drop = FALSE], y = y[!oo, , drop = FALSE], alpha = alpha, my_lambda = lambda, lambda_min = 0.01, max_it = max_it, e.abs = e.abs, e.rel = e.rel, nlambda = length(lambda[, 1]), rho = rho, tree = TT, my_print = my_print, alph = alph, parallel = parallel, pal = pal, gg = gg, tol = tol, cl = cl, legacy)

R/hh_nlambda_loop.R

@@ -0,0 +1,102 @@
+hh_nlambda_loop <- function(
+  lam, nlambda, beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it,
+  my_W_hat, XtY, y, N, e.abs, e.rel, alpha, alph, svd.w, tree, my_print,
+  invmat, gg, tol, parallel, pal, BETA0, THETA0, BETA,
+  BETA_hat, Y_HAT, THETA, D, my_values, legacy = TRUE
+) {
+  if (legacy) {
+    obj <- NULL
+    non_zero_theta <- NULL
+    my_obj <- list()
+    n_main_terms <- NULL
+    lam_list <- list()
+    hh <- 1
+    while (hh <= nlambda) {
+      lambda <- lam[hh, ]
+
+      start_time <- Sys.time()
+      if (pal) {
+        my_values <- admm_MADMMplasso(
+          beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it, my_W_hat, XtY,
+          y, N, e.abs, e.rel, alpha, lambda, alph, svd.w, tree, my_print, invmat,
+          gg[hh, ], legacy
+        )
+
+        beta <- my_values$beta
+        theta <- my_values$theta
+        my_obj[[hh]] <- list(my_values$obj)
+        beta0 <- my_values$beta0
+        theta0 <- my_values$theta0 ### iteration
+        beta_hat <- my_values$beta_hat
+        y_hat <- my_values$y_hat
+      }
+      cost_time <- Sys.time() - start_time
+      if (my_print) {
+        print(cost_time)
+      }
+      if (parallel && !pal) {
+        beta <- my_values[hh, ]$beta
+        theta <- my_values[hh, ]$theta
+        my_obj[[hh]] <- list(my_values[hh, ]$obj)
+        beta0 <- my_values[hh, ]$beta0
+        theta0 <- my_values[hh, ]$theta0 ### iteration
+        beta_hat <- my_values[hh, ]$beta_hat
+        y_hat <- my_values[hh, ]$y_hat
+      } else if (!parallel && !pal) {
+        beta <- my_values[[hh]]$beta
+        theta <- my_values[[hh]]$theta
+        my_obj[[hh]] <- list(my_values[[hh]]$obj)
+        beta0 <- my_values[[hh]]$beta0
+        theta0 <- my_values[[hh]]$theta0 ### iteration
+        beta_hat <- my_values[[hh]]$beta_hat
+        y_hat <- my_values[[hh]]$y_hat
+      }
+      # Executed if par == TRUE, independent of parallel
+
+      beta1 <- as(beta * (abs(beta) > tol), "sparseMatrix")
+      theta1 <- as.sparse3Darray(theta * (abs(theta) > tol))
+      beta_hat1 <- as(beta_hat * (abs(beta_hat) > tol), "sparseMatrix")
+
+      n_interaction_terms <- count_nonzero_a((theta1))
+
+      n_main_terms <- (c(n_main_terms, count_nonzero_a((beta1))))
+
+      obj1 <- (sum(as.vector((y - y_hat)^2))) / (D * N)
+      obj <- c(obj, obj1)
+
+      non_zero_theta <- (c(non_zero_theta, n_interaction_terms))
+      lam_list <- (c(lam_list, lambda))
+
+      BETA0[[hh]] <- beta0
+      THETA0[[hh]] <- theta0
+      BETA[[hh]] <- as(beta1, "sparseMatrix")
+      BETA_hat[[hh]] <- as(beta_hat1, "sparseMatrix")
+
+      Y_HAT[[hh]] <- y_hat
+      THETA[[hh]] <- as.sparse3Darray(theta1)
+
+      if (my_print) {
+        if (hh == 1) {
+          print(c(hh, (n_main_terms[hh]), non_zero_theta[hh], obj1))
+        } else {
+          print(c(hh, (n_main_terms[hh]), non_zero_theta[hh], obj[hh - 1], obj1))
+        }
+      }
+
+      hh <- hh + 1
+    } ### lambda
+    out <- list(
+      obj = obj, n_main_terms = n_main_terms, non_zero_theta = non_zero_theta,
+      BETA0 = BETA0, THETA0 = THETA0, BETA = BETA, BETA_hat = BETA_hat,
+      Y_HAT = Y_HAT, THETA = THETA
+    )
+  } else {
+    out <- hh_nlambda_loop_cpp(
+      lam, nlambda, beta0, theta0, beta, beta_hat, theta, rho1, X, Z, max_it,
+      my_W_hat, XtY, y, N, e.abs, e.rel, alpha, alph, svd.w, tree, my_print,
+      invmat, gg, tol, parallel, pal, BETA0, THETA0, BETA,
+      BETA_hat, Y_HAT, THETA, D, my_values
+    )
+  }
+  return(out)
+}

src/MADMMplasso.h

@@ -0,0 +1,26 @@
+#include <RcppArmadillo.h>
+Rcpp::List admm_MADMMplasso_cpp(
+  const arma::vec beta0,
+  const arma::mat theta0,
+  arma::mat beta,
+  arma::mat beta_hat,
+  arma::cube theta,
+  const double rho1,
+  const arma::mat X,
+  const arma::mat Z,
+  const int max_it,
+  const arma::mat W_hat,
+  const arma::mat XtY,
+  const arma::mat y,
+  const int N,
+  const double e_abs,
+  const double e_rel,
+  const double alpha,
+  const arma::vec lambda,
+  const double alph,
+  const Rcpp::List svd_w,
+  const Rcpp::List tree,
+  const Rcpp::List invmat,
+  const arma::vec gg,
+  const bool my_print = true
+);