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