pkg/R/main.R

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