pkg/R/constructionModelesLassoRank.R

   1 #' constructionModelesLassoRank
   2 #'
   3 #' Construct a collection of models with the Lasso-Rank procedure.
   4 #'
   5 #' @param S output of selectVariables.R
   6 #' @param k number of components
   7 #' @param mini integer, minimum number of iterations in the EM algorithm, by default = 10
   8 #' @param maxi integer, maximum number of iterations in the EM algorithm, by default = 100
   9 #' @param X matrix of covariates (of size n*p)
  10 #' @param Y matrix of responses (of size n*m)
  11 #' @param eps real, threshold to say the EM algorithm converges, by default = 1e-4
  12 #' @param rank.min integer, minimum rank in the low rank procedure, by default = 1
  13 #' @param rank.max integer, maximum rank in the low rank procedure, by default = 5
  14 #' @param ncores Number of cores, by default = 3
  15 #' @param fast TRUE to use compiled C code, FALSE for R code only
  16 #' @param verbose TRUE to show some execution traces
  17 #'
  18 #' @return a list with several models, defined by phi, rho, pi, llh
  19 #'
  20 #' @export
  21 constructionModelesLassoRank = function(S, k, mini, maxi, X, Y, eps, rank.min,
  22                                         rank.max, ncores, fast=TRUE, verbose=FALSE)
  23 {
  24   n = dim(X)[1]
  25   p = dim(X)[2]
  26   m = dim(Y)[2]
  27   L = length(S)
  28
  29   # Possible interesting ranks
  30   deltaRank = rank.max - rank.min + 1
  31   Size = deltaRank^k
  32   RankLambda = matrix(0, nrow=Size*L, ncol=k+1)
  33   for (r in 1:k)
  34   {
  35     # On veut le tableau de toutes les combinaisons de rangs possibles, et des lambdas
  36     # Dans la première colonne : on répète (rank.max-rank.min)^(k-1) chaque chiffre :
  37     #   ça remplit la colonne
  38     # Dans la deuxieme : on répète (rank.max-rank.min)^(k-2) chaque chiffre,
  39     #   et on fait ça (rank.max-rank.min)^2 fois
  40     # ...
  41     # Dans la dernière, on répète chaque chiffre une fois,
  42     #   et on fait ça (rank.min-rank.max)^(k-1) fois.
  43     RankLambda[,r] = rep(rank.min + rep(0:(deltaRank-1), deltaRank^(r-1), each=deltaRank^(k-r)), each = L)
  44   }
  45   RankLambda[,k+1] = rep(1:L, times = Size)
  46
  47   if (ncores > 1)
  48   {
  49     cl = parallel::makeCluster(ncores, outfile='')
  50     parallel::clusterExport( cl, envir=environment(),
  51                              varlist=c("A1","Size","Pi","Rho","mini","maxi","X","Y","eps",
  52                                        "Rank","m","phi","ncores","verbose") )
  53   }
  54
  55   computeAtLambda <- function(index)
  56   {
  57     lambdaIndex = RankLambda[index,k+1]
  58     rankIndex = RankLambda[index,1:k]
  59     if (ncores > 1)
  60       require("valse") #workers start with an empty environment
  61
  62     # 'relevant' will be the set of relevant columns
  63     selected = S[[lambdaIndex]]$selected
  64     relevant = c()
  65     for (j in 1:p){
  66       if (length(selected[[j]])>0){
  67         relevant = c(relevant,j)
  68       }
  69     }
  70     if (max(rankIndex)<length(relevant)){
  71       phi = array(0, dim=c(p,m,k))
  72       if (length(relevant) > 0)
  73       {
  74         res = EMGrank(S[[lambdaIndex]]$Pi, S[[lambdaIndex]]$Rho, mini, maxi,
  75                       X[,relevant], Y, eps, rankIndex, fast)
  76         llh = c( res$LLF, sum(rankIndex * (length(relevant)- rankIndex + m)) )
  77         phi[relevant,,] = res$phi
  78       }
  79       list("llh"=llh, "phi"=phi, "pi" = S[[lambdaIndex]]$Pi, "rho" = S[[lambdaIndex]]$Rho)
  80
  81     }
  82   }
  83
  84   #For each lambda in the grid we compute the estimators
  85   out =
  86     if (ncores > 1)
  87       parLapply(cl, seq_len(length(S)*Size), computeAtLambda)
  88   else
  89     lapply(seq_len(length(S)*Size), computeAtLambda)
  90
  91   if (ncores > 1)
  92     parallel::stopCluster(cl)
  93
  94   out
  95 }