remove extra prints
[valse.git] / pkg / R / constructionModelesLassoMLE.R
1 #' constructionModelesLassoMLE
2 #'
3 #' TODO: description
4 #'
5 #' @param ...
6 #'
7 #' @return ...
8 #'
9 #' export
10 constructionModelesLassoMLE = function(phiInit, rhoInit, piInit, gamInit, mini, maxi,
11 gamma, X, Y, thresh, tau, S, ncores=3, artefact = 1e3, fast=TRUE, verbose=FALSE)
12 {
13 if (ncores > 1)
14 {
15 cl = parallel::makeCluster(ncores, outfile='')
16 parallel::clusterExport( cl, envir=environment(),
17 varlist=c("phiInit","rhoInit","gamInit","mini","maxi","gamma","X","Y","thresh",
18 "tau","S","ncores","verbose") )
19 }
20
21 # Individual model computation
22 computeAtLambda <- function(lambda)
23 {
24 if (ncores > 1)
25 require("valse") #nodes start with an empty environment
26
27 if (verbose)
28 print(paste("Computations for lambda=",lambda))
29
30 n = dim(X)[1]
31 p = dim(phiInit)[1]
32 m = dim(phiInit)[2]
33 k = dim(phiInit)[3]
34
35 sel.lambda = S[[lambda]]$selected
36 # col.sel = which(colSums(sel.lambda)!=0) #if boolean matrix
37 col.sel <- which( sapply(sel.lambda,length) > 0 ) #if list of selected vars
38
39 if (length(col.sel) == 0)
40 return (NULL)
41
42 # lambda == 0 because we compute the EMV: no penalization here
43 res = EMGLLF(phiInit[col.sel,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0,
44 X[,col.sel], Y, tau, fast)
45
46 # Eval dimension from the result + selected
47 phiLambda2 = res$phi
48 rhoLambda = res$rho
49 piLambda = res$pi
50 phiLambda = array(0, dim = c(p,m,k))
51 for (j in seq_along(col.sel))
52 phiLambda[col.sel[j],,] = phiLambda2[j,,]
53 dimension = length(unlist(sel.lambda))
54
55 # Computation of the loglikelihood
56 densite = vector("double",n)
57 for (r in 1:k)
58 {
59 delta = (Y%*%rhoLambda[,,r] - (X[, col.sel]%*%phiLambda[col.sel,,r]))/artefact
60 densite = densite + piLambda[r] *
61 det(rhoLambda[,,r])/(sqrt(2*base::pi))^m * exp(-tcrossprod(delta)/2.0)
62 }
63 llhLambda = c( sum(artefact^2 * log(densite)), (dimension+m+1)*k-1 )
64 list("phi"= phiLambda, "rho"= rhoLambda, "pi"= piLambda, "llh" = llhLambda)
65 }
66
67 # For each lambda, computation of the parameters
68 out =
69 if (ncores > 1)
70 parLapply(cl, 1:length(S), computeAtLambda)
71 else
72 lapply(1:length(S), computeAtLambda)
73
74 if (ncores > 1)
75 parallel::stopCluster(cl)
76
77 out
78 }