From: Benjamin Auder <benjamin.auder@somewhere>
Date: Tue, 7 Apr 2020 13:58:56 +0000 (+0200)
Subject: Fix constructionModelesLassoMLE: phiInit was reshaped from array to matrix. Add examp... 
X-Git-Url: https://git.auder.net/game/current/%7B%7B%20asset%28%27mixstore/css/base.css%27%29%20%7D%7D?a=commitdiff_plain;h=fb3557f39487d9631ffde30f20b70938d2a6ab0c;p=valse.git

Fix constructionModelesLassoMLE: phiInit was reshaped from array to matrix. Add examples in main.R
---

diff --git a/TODO b/TODO
deleted file mode 100644
index c662bec..0000000
--- a/TODO
+++ /dev/null
@@ -1,29 +0,0 @@
-n = 50; m = 10; p = 5
-X = matrix(runif(n*p, -10, 10), nrow=n)
-Y = matrix(runif(n*m, -5, 15), nrow=n)
-beta = array(0, dim=c(p,m,2))
-beta[,,1] = 1
-beta[,,2] = 2
-data = generateXY(n, c(0.4,0.6), rep(0,p), beta, diag(0.5, p), diag(0.5, m))
-X = data$X
-Y = data$Y
-class = data$class
-V1 = runValse(X, Y, fast=FALSE)
-Error in while (!pi2AllPositive) { : 
-  missing value where TRUE/FALSE needed
-
-V2 = runValse(X, Y, fast=TRUE)
-list()
-Error in out[[ind_uniq[l]]] : 
-  attempt to select less than one element in get1index
-
-==> Error identified: line 61 in initSmallEM.R, division by 0
-It occurs also for smallers values of n and m, e.g.: n = 20; m = 20; p = 3
-
-=====
-
-Also:
-X <- matrix(runif(100), nrow=50)
-Y <- matrix(runif(100), nrow=50)
-(...)
-Error: cannot allocate vector of size 16.0 Gb
diff --git a/pkg/R/constructionModelesLassoMLE.R b/pkg/R/constructionModelesLassoMLE.R
index 2d04adb..0584382 100644
--- a/pkg/R/constructionModelesLassoMLE.R
+++ b/pkg/R/constructionModelesLassoMLE.R
@@ -21,7 +21,7 @@
 #'
 #' @export
 constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini,
-  maxi, gamma, X, Y, eps, S, ncores = 3, fast, verbose)
+  maxi, gamma, X, Y, eps, S, ncores, fast, verbose)
 {
   if (ncores > 1)
   {
@@ -51,8 +51,9 @@ constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini,
       return(NULL)
 
     # lambda == 0 because we compute the EMV: no penalization here
-    res <- EMGLLF(array(phiInit,dim=c(p,m,k))[col.sel, , ], rhoInit, piInit, gamInit,
-      mini, maxi, gamma, 0, as.matrix(X[, col.sel]), Y, eps, fast)
+    res <- EMGLLF(array(phiInit[col.sel, , ], dim=c(length(col.sel),m,k)),
+      rhoInit, piInit, gamInit, mini, maxi, gamma, 0,
+      as.matrix(X[, col.sel]), Y, eps, fast)
 
     # Eval dimension from the result + selected
     phiLambda2 <- res$phi
diff --git a/pkg/R/main.R b/pkg/R/main.R
index c74d7fb..13df89f 100644
--- a/pkg/R/main.R
+++ b/pkg/R/main.R
@@ -23,10 +23,22 @@
 #' @param verbose TRUE to show some execution traces
 #' @param plot TRUE to plot the selected models after run
 #'
-#' @return a list with estimators of parameters
+#' @return
+#' The selected model if enough data are available to estimate it,
+#' or a list of models otherwise.
 #'
 #' @examples
-#' #TODO: a few examples
+#' n = 50; m = 10; p = 5
+#' beta = array(0, dim=c(p,m,2))
+#' beta[,,1] = 1
+#' beta[,,2] = 2
+#' data = generateXY(n, c(0.4,0.6), rep(0,p), beta, diag(0.5, p), diag(0.5, m))
+#' X = data$X
+#' Y = data$Y
+#' res = runValse(X, Y)
+#' X <- matrix(runif(100), nrow=50)
+#' Y <- matrix(runif(100), nrow=50)
+#' res = runValse(X, Y)
 #'
 #' @export
 runValse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mini = 10,
@@ -125,30 +137,31 @@ runValse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1,
       complexity = sumPen, contrast = -LLH)
   }))
   tableauRecap <- tableauRecap[which(tableauRecap[, 4] != Inf), ]
-
-  if (verbose == TRUE)
+  if (verbose)
     print(tableauRecap)
-  modSel <- capushe::capushe(tableauRecap, n)
-  indModSel <- if (selecMod == "DDSE")
-  {
-    as.numeric(modSel@DDSE@model)
-  } else if (selecMod == "Djump")
-  {
-    as.numeric(modSel@Djump@model)
-  } else if (selecMod == "BIC")
-  {
-    modSel@BIC_capushe$model
-  } else if (selecMod == "AIC")
-  {
-    modSel@AIC_capushe$model
-  }
 
-  listMod <- as.integer(unlist(strsplit(as.character(indModSel), "[.]")))
-  modelSel <- models_list[[listMod[1]]][[listMod[2]]]
-  modelSel$tableau <- tableauRecap
-
-  if (plot)
-    print(plot_valse(X, Y, modelSel))
+  if (nrow(tableauRecap) > 10) {
+    modSel <- capushe::capushe(tableauRecap, n)
+    indModSel <- if (selecMod == "DDSE")
+    {
+      as.numeric(modSel@DDSE@model)
+    } else if (selecMod == "Djump")
+    {
+      as.numeric(modSel@Djump@model)
+    } else if (selecMod == "BIC")
+    {
+      modSel@BIC_capushe$model
+    } else if (selecMod == "AIC")
+    {
+      modSel@AIC_capushe$model
+    }
+    listMod <- as.integer(unlist(strsplit(as.character(indModSel), "[.]")))
+    modelSel <- models_list[[listMod[1]]][[listMod[2]]]
+    modelSel$models <- tableauRecap
 
-  return(modelSel)
+    if (plot)
+      print(plot_valse(X, Y, modelSel))
+    return(modelSel)
+  }
+  tableauRecap
 }
diff --git a/pkg/R/selectVariables.R b/pkg/R/selectVariables.R
index 99959ca..2d1c9b7 100644
--- a/pkg/R/selectVariables.R
+++ b/pkg/R/selectVariables.R
@@ -66,7 +66,6 @@ selectVariables <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma
   if (ncores > 1)
     parallel::stopCluster(cl)
 
-  print(out) #DEBUG TRACE
   # Suppress models which are computed twice
   # sha1_array <- lapply(out, digest::sha1) out[ duplicated(sha1_array) ]
   selec <- lapply(out, function(model) model$selected)