projects / valse.git / blame
| Commit | Line | Data |
|---|---|---|
| 2279a641 BA |
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, | |
| aa480ac1 | 11 | rangmax, ncores, fast=TRUE, verbose=FALSE) |
| ef67d338 | 12 | { |
| 3f145e9a BG |
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] | |
| ef67d338 BA |
18 | |
| 19 | # On cherche les rangs possiblement intéressants | |
| 3f145e9a BG |
20 | deltaRank = rangmax - rangmin + 1 |
| 21 | Size = deltaRank^k | |
| ef67d338 | 22 | Rank = matrix(0, nrow=Size, ncol=k) |
| 0eb161e3 | 23 | for (r in 1:k) |
| ef67d338 BA |
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 | ||
| 0eb161e3 BA |
36 | if (ncores > 1) |
| 37 | { | |
| b9b0b72a | 38 | cl = parallel::makeCluster(ncores, outfile='') |
| 0eb161e3 BA |
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 | } | |
| 2279a641 | 43 | |
| 0eb161e3 | 44 | computeAtLambda <- function(lambdaIndex) |
| ef67d338 | 45 | { |
| 0eb161e3 BA |
46 | if (ncores > 1) |
| 47 | require("valse") #workers start with an empty environment | |
| 48 | ||
| ef67d338 BA |
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)] | |
| 0eb161e3 BA |
53 | phi = array(0, dim=c(p,m,k,Size)) |
| 54 | llh = matrix(0, Size, 2) #log-likelihood | |
| ef67d338 BA |
55 | if (length(active) > 0) |
| 56 | { | |
| 57 | for (j in 1:Size) | |
| 58 | { | |
| 59 | res = EMGrank(Pi[,lambdaIndex], Rho[,,,lambdaIndex], mini, maxi, | |
| aa480ac1 | 60 | X[,active], Y, tau, Rank[j,], fast) |
| 0eb161e3 BA |
61 | llh = rbind(llh, |
| 62 | c( res$LLF, sum(Rank[j,] * (length(active)- Rank[j,] + m)) ) ) | |
| 63 | phi[active,,,] = rbind(phi[active,,,], res$phi) | |
| 3f145e9a BG |
64 | } |
| 65 | } | |
| 0eb161e3 BA |
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) | |
| 9ade3f1b | 84 | } |