pkg/R/main.R

   1 #' runValse
   2 #'
   3 #' Main function
   4 #'
   5 #' @param X matrix of covariates (of size n*p)
   6 #' @param Y matrix of responses (of size n*m)
   7 #' @param procedure among 'LassoMLE' or 'LassoRank'
   8 #' @param selecMod method to select a model among 'DDSE', 'DJump', 'BIC' or 'AIC'
   9 #' @param gamma integer for the power in the penaly, by default = 1
  10 #' @param mini integer, minimum number of iterations in the EM algorithm, by default = 10
  11 #' @param maxi integer, maximum number of iterations in the EM algorithm, by default = 100
  12 #' @param eps real, threshold to say the EM algorithm converges, by default = 1e-4
  13 #' @param kmin integer, minimum number of clusters, by default = 2
  14 #' @param kmax integer, maximum number of clusters, by default = 10
  15 #' @param rank.min integer, minimum rank in the low rank procedure, by default = 1
  16 #' @param rank.max integer, maximum rank in the low rank procedure, by default = 5
  17 #' @param ncores_outer Number of cores for the outer loop on k
  18 #' @param ncores_inner Number of cores for the inner loop on lambda
  19 #' @param thresh real, threshold to say a variable is relevant, by default = 1e-8
  20 #' @param grid_lambda, a vector with regularization parameters if known, by default numeric(0)
  21 #' @param size_coll_mod (Maximum) size of a collection of models
  22 #' @param fast TRUE to use compiled C code, FALSE for R code only
  23 #' @param verbose TRUE to show some execution traces
  24 #' @param plot TRUE to plot the selected models after run
  25 #'
  26 #' @return a list with estimators of parameters
  27 #'
  28 #' @examples
  29 #' #TODO: a few examples
  30 #'
  31 #' @export
  32 runValse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mini = 10,
  33   maxi = 50, eps = 1e-04, kmin = 2, kmax = 3, rank.min = 1, rank.max = 5, ncores_outer = 1,
  34   ncores_inner = 1, thresh = 1e-08, grid_lambda = numeric(0), size_coll_mod = 10,
  35   fast = TRUE, verbose = FALSE, plot = TRUE)
  36 {
  37   n <- nrow(X)
  38   p <- ncol(X)
  39   m <- ncol(Y)
  40
  41   if (verbose)
  42     print("main loop: over all k and all lambda")
  43
  44   if (ncores_outer > 1) {
  45     cl <- parallel::makeCluster(ncores_outer, outfile = "")
  46     parallel::clusterExport(cl = cl, envir = environment(), varlist = c("X",
  47       "Y", "procedure", "selecMod", "gamma", "mini", "maxi", "eps", "kmin",
  48       "kmax", "rank.min", "rank.max", "ncores_outer", "ncores_inner", "thresh",
  49       "size_coll_mod", "verbose", "p", "m"))
  50   }
  51
  52   # Compute models with k components
  53   computeModels <- function(k)
  54   {
  55     if (ncores_outer > 1)
  56       require("valse") #nodes start with an empty environment
  57
  58     if (verbose)
  59       print(paste("Parameters initialization for k =", k))
  60     # smallEM initializes parameters by k-means and regression model in each
  61     # component, doing this 20 times, and keeping the values maximizing the
  62     # likelihood after 10 iterations of the EM algorithm.
  63     P <- initSmallEM(k, X, Y, fast)
  64     if (length(grid_lambda) == 0)
  65     {
  66       grid_lambda <- computeGridLambda(P$phiInit, P$rhoInit, P$piInit, P$gamInit,
  67                                        X, Y, gamma, mini, maxi, eps, fast)
  68     }
  69     if (length(grid_lambda) > size_coll_mod)
  70       grid_lambda <- grid_lambda[seq(1, length(grid_lambda), length.out = size_coll_mod)]
  71
  72     if (verbose)
  73       print("Compute relevant parameters")
  74     # select variables according to each regularization parameter from the grid:
  75     # S$selected corresponding to selected variables
  76     S <- selectVariables(P$phiInit, P$rhoInit, P$piInit, P$gamInit, mini, maxi,
  77       gamma, grid_lambda, X, Y, thresh, eps, ncores_inner, fast)
  78
  79     if (procedure == "LassoMLE") {
  80       if (verbose)
  81         print("run the procedure Lasso-MLE")
  82       # compute parameter estimations, with the Maximum Likelihood Estimator,
  83       # restricted on selected variables.
  84       models <- constructionModelesLassoMLE(P$phiInit, P$rhoInit, P$piInit,
  85         P$gamInit, mini, maxi, gamma, X, Y, eps, S, ncores_inner, fast, verbose)
  86     } else {
  87       if (verbose)
  88         print("run the procedure Lasso-Rank")
  89       # compute parameter estimations, with the Low Rank Estimator, restricted on
  90       # selected variables.
  91       models <- constructionModelesLassoRank(S, k, mini, maxi, X, Y, eps, rank.min,
  92         rank.max, ncores_inner, fast, verbose)
  93     }
  94     # warning! Some models are NULL after running selectVariables
  95     models <- models[sapply(models, function(cell) !is.null(cell))]
  96     models
  97   }
  98
  99   # List (index k) of lists (index lambda) of models
 100   models_list <-
 101     if (ncores_outer > 1) {
 102       parLapply(cl, kmin:kmax, computeModels)
 103     } else {
 104       lapply(kmin:kmax, computeModels)
 105     }
 106   if (ncores_outer > 1)
 107     parallel::stopCluster(cl)
 108
 109   if (!requireNamespace("capushe", quietly = TRUE))
 110   {
 111     warning("'capushe' not available: returning all models")
 112     return(models_list)
 113   }
 114
 115   # Get summary 'tableauRecap' from models
 116   tableauRecap <- do.call(rbind, lapply(seq_along(models_list), function(i)
 117   {
 118     models <- models_list[[i]]
 119     # For a collection of models (same k, several lambda):
 120     LLH <- sapply(models, function(model) model$llh[1])
 121     k <- length(models[[1]]$pi)
 122     sumPen <- sapply(models, function(model) k * (dim(model$rho)[1] + sum(model$phi[,
 123       , 1] != 0) + 1) - 1)
 124     data.frame(model = paste(i, ".", seq_along(models), sep = ""), pen = sumPen/n,
 125       complexity = sumPen, contrast = -LLH)
 126   }))
 127   tableauRecap <- tableauRecap[which(tableauRecap[, 4] != Inf), ]
 128
 129   if (verbose == TRUE)
 130     print(tableauRecap)
 131   modSel <- capushe::capushe(tableauRecap, n)
 132   indModSel <- if (selecMod == "DDSE")
 133   {
 134     as.numeric(modSel@DDSE@model)
 135   } else if (selecMod == "Djump")
 136   {
 137     as.numeric(modSel@Djump@model)
 138   } else if (selecMod == "BIC")
 139   {
 140     modSel@BIC_capushe$model
 141   } else if (selecMod == "AIC")
 142   {
 143     modSel@AIC_capushe$model
 144   }
 145
 146   listMod <- as.integer(unlist(strsplit(as.character(indModSel), "[.]")))
 147   modelSel <- models_list[[listMod[1]]][[listMod[2]]]
 148   modelSel$tableau <- tableauRecap
 149
 150   if (plot)
 151     print(plot_valse(X, Y, modelSel, n))
 152
 153   return(modelSel)
 154 }