pkg/R/constructionModelesLassoRank.R

   1 #' constructionModelesLassoRank
   2 #'
   3 #' TODO: description
   4 #'
   5 #' @param ...
   6 #'
   7 #' @return ...
   8 #'
   9 #' export
  10 constructionModelesLassoRank = function(pi, rho, mini, maxi, X, Y, tau, A1, rangmin,
  11     rangmax, ncores, fast=TRUE, verbose=FALSE)
  12 {
  13   n = dim(X)[1]
  14   p = dim(X)[2]
  15   m = dim(rho)[2]
  16   k = dim(rho)[3]
  17   L = dim(A1)[2]
  18
  19     # On cherche les rangs possiblement intéressants
  20   deltaRank = rangmax - rangmin + 1
  21   Size = deltaRank^k
  22   Rank = matrix(0, nrow=Size, ncol=k)
  23   for (r in 1:k)
  24     {
  25         # On veut le tableau de toutes les combinaisons de rangs possibles
  26         # Dans la première colonne : on répète (rangmax-rangmin)^(k-1) chaque chiffre :
  27         #   ça remplit la colonne
  28         # Dans la deuxieme : on répète (rangmax-rangmin)^(k-2) chaque chiffre,
  29         #   et on fait ça (rangmax-rangmin)^2 fois
  30         # ...
  31         # Dans la dernière, on répète chaque chiffre une fois,
  32         #   et on fait ça (rangmin-rangmax)^(k-1) fois.
  33     Rank[,r] = rangmin + rep(0:(deltaRank-1), deltaRank^(r-1), each=deltaRank^(k-r))
  34   }
  35
  36   if (ncores > 1)
  37     {
  38     cl = parallel::makeCluster(ncores, outfile='')
  39     parallel::clusterExport( cl, envir=environment(),
  40             varlist=c("A1","Size","Pi","Rho","mini","maxi","X","Y","tau",
  41             "Rank","m","phi","ncores","verbose") )
  42     }
  43
  44     computeAtLambda <- function(lambdaIndex)
  45     {
  46         if (ncores > 1)
  47             require("valse") #workers start with an empty environment
  48
  49     # on ne garde que les colonnes actives
  50     # 'active' sera l'ensemble des variables informatives
  51     active = A1[,lambdaIndex]
  52     active = active[-(active==0)]
  53         phi = array(0, dim=c(p,m,k,Size))
  54         llh = matrix(0, Size, 2) #log-likelihood
  55     if (length(active) > 0)
  56         {
  57       for (j in 1:Size)
  58             {
  59         res = EMGrank(Pi[,lambdaIndex], Rho[,,,lambdaIndex], mini, maxi,
  60                     X[,active], Y, tau, Rank[j,], fast)
  61         llh = rbind(llh,
  62                     c( res$LLF, sum(Rank[j,] * (length(active)- Rank[j,] + m)) ) )
  63         phi[active,,,] = rbind(phi[active,,,], res$phi)
  64       }
  65     }
  66         list("llh"=llh, "phi"=phi)
  67     }
  68
  69     #Pour chaque lambda de la grille, on calcule les coefficients
  70   out =
  71         if (ncores > 1)
  72             parLapply(cl, seq_along(glambda), computeAtLambda)
  73         else
  74             lapply(seq_along(glambda), computeAtLambda)
  75
  76     if (ncores > 1)
  77     parallel::stopCluster(cl)
  78
  79     # TODO: this is a bit ugly. Better use bigmemory and fill llh/phi in-place
  80     # (but this also adds a dependency...)
  81     llh <- do.call( rbind, lapply(out, function(model) model$llh) )
  82     phi <- do.call( rbind, lapply(out, function(model) model$phi) )
  83     list("llh"=llh, "phi"=phi)
  84 }