| 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 | } |