Update main.R
[valse.git] / pkg / R / constructionModelesLassoRank.R
CommitLineData
2279a641
BA
1#' constructionModelesLassoRank
2#'
3#' TODO: description
4#'
5#' @param ...
6#'
7#' @return ...
8#'
9#' export
10constructionModelesLassoRank = 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}