OOB issue on C++ code

[wleoncio]

Code for reproduction: https://gist.github.com/wleoncio/0318b15226128e6aef7898d4a2c009e9

Error message:

Error in eval(ei, envir) : field::operator(): index out of bounds

[wleoncio]

Problem seems to be malformation of my_values when legacy = FALSE. Here's what happens to str(my_values) when nlambda = 1L:

[ins] Browse[1]> str(my_values)
List of 1
 $ :List of 8
  ..$ beta0   : num [1:6] -0.6065 -0.5524 0.2878 0.3876 -0.0607 ...
  ..$ theta0  : num [1:4, 1:6] 1.054 0.173 0.546 0.863 0.976 ...
  ..$ beta    : num [1:50, 1:6] 0.355 0 0 0.274 0.359 ...
  ..$ theta   : num [1:50, 1:4, 1:6] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ converge: logi TRUE
  ..$ obj     : NULL
  ..$ beta_hat: num [1:250, 1:6] 0.355 0 0 0.274 0.359 ...
  ..$ y_hat   : num [1:100, 1:6] -0.927 -0.1251 -1.8735 -0.0907 0.645 ...

When nlamda = 2L, we get a list of lists of the same size, i.e.:

[ins] Browse[1]> str(my_values)
List of 2
 $ :List of 8
  ..$ beta0   : num [1:6] -0.6065 -0.5524 0.2878 0.3876 -0.0607 ...
  ..$ theta0  : num [1:4, 1:6] 1.054 0.173 0.546 0.863 0.976 ...
  ..$ beta    : num [1:50, 1:6] 0.355 0 0 0.274 0.359 ...
  ..$ theta   : num [1:50, 1:4, 1:6] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ converge: logi TRUE
  ..$ obj     : NULL
  ..$ beta_hat: num [1:250, 1:6] 0.355 0 0 0.274 0.359 ...
  ..$ y_hat   : num [1:100, 1:6] -0.927 -0.1251 -1.8735 -0.0907 0.645 ...
 $ :List of 8
  ..$ beta0   : num [1:6] -0.0799 -0.0809 0.1453 0.2028 -0.0286 ...
  ..$ theta0  : num [1:4, 1:6] -0.0473 -0.094 0.2274 0.139 -0.0811 ...
  ..$ beta    : num [1:50, 1:6] 1.92 1.71 1.71 1.83 1.85 ...
  ..$ theta   : num [1:50, 1:4, 1:6] 1.9265 0.0517 -0.2163 0.1983 0.3255 ...
  ..$ converge: logi TRUE
  ..$ obj     : NULL
  ..$ beta_hat: num [1:250, 1:6] 1.92 1.71 1.71 1.83 1.85 ...
  ..$ y_hat   : num [1:100, 1:6] 6.19 -1.42 -2.09 -5.25 2.33 ...

When legacy = FALSE, we get a proper my_values (size is 7 instead of 8 because C++ doesn't output obj, which is always NULL):

[ins] Browse[1]> str(my_values)
List of 1
 $ :List of 7
  ..$ : num [1:6, 1, 1] -0.607 -0.5532 0.2879 0.3877 -0.0606 ...
  ..$ : num [1:4, 1:6, 1] 1.049 0.161 0.54 0.862 0.975 ...
  ..$ : num [1:50, 1:6, 1] 0.393 0 0 0.28 0.4 ...
  ..$ : num [1:50, 1:4, 1:6] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ : num [1, 1, 1] 1
  ..$ : num [1:250, 1:6, 1] 0.393 0 0 0.28 0.4 ...
  ..$ : num [1:100, 1:6, 1] -0.969 -0.142 -1.881 -0.101 0.676 ...
  ..- attr(*, "dim")= int [1:2] 7 1

However, if nlambda = 2L, output is only beta0, leading to the OOB errors on these lines

[ins] Browse[1]> str(my_values)
List of 2
 $ :List of 1
  ..$ : num [1:6, 1, 1] -0.607 -0.5532 0.2879 0.3877 -0.0606 ...
 $ :List of 1
  ..$ : num [1:4, 1:6, 1] 1.049 0.161 0.54 0.862 0.975 ...

So the problem is possibly around here, I'll investigate further:

MADMMplasso/R/MADMMplasso.R

Lines 206 to 262 in 081383c

	if (parallel) {
	cl <- makeCluster(cl1, type = "FORK")
	doParallel::registerDoParallel(cl = cl)
	foreach::getDoParRegistered()
	if (legacy) {
	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, ]
	)
	}
	} else {
	my_values_matrix <- foreach(i = 1:nlambda, .packages = "MADMMplasso", .combine = rbind) %dopar% {
	admm_MADMMplasso_cpp(
	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_tu, svd_w_tv, svd_w_d,
	C, CW, gg[i, ], my_print
	)
	}
	}
	parallel::stopCluster(cl)
	# Converting to list so hh_nlambda_loop_cpp can handle it
	if (nlambda == 1) {
	my_values <- list(my_values_matrix)
	} else {
	my_values <- list()
	for (hh in seq_len(nlambda)) {
	my_values[[hh]] <- my_values_matrix[hh, ]
	}
	}
	} else if (!parallel && !pal) {
	if (legacy) {
	my_values <- lapply(
	seq_len(nlambda),
	function(g) {
	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[g, ], alph, svd.w, tree, my_print,
	invmat, gg[g, ]
	)
	}
	)
	} else {
	my_values <- lapply(
	seq_len(nlambda),
	function(i) {
	admm_MADMMplasso_cpp(
	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_tu, svd_w_tv, svd_w_d,
	C, CW, gg[i, ], my_print
	)
	}
	)
	}
	} else {

[Theo-qua]

Alright. Thank you for the update.

[wleoncio]

Almost done. Another question: is is OK for the output non_zero to be slightly different between engines (C++ gives {3, 17}, R gives {3, 25} for a run with nlambda = 2)?

[Theo-qua]

This is okay