From 43d76c49d2f98490abc782c7e8a8b94baee40247 Mon Sep 17 00:00:00 2001
From: emilie <emilie@devijver.org>
Date: Thu, 13 Apr 2017 13:52:32 +0200
Subject: [PATCH] fix EMGRank.R, and add some lines in the roxygen code for
 some functions

---
 pkg/DESCRIPTION                      |   3 +-
 pkg/R/A_NAMESPACE.R                  |   2 +-
 pkg/R/EMGLLF.R                       |  24 ++---
 pkg/R/EMGrank.R                      |   8 +-
 pkg/R/constructionModelesLassoMLE.R  |  32 ++++--
 pkg/R/constructionModelesLassoRank.R | 145 ++++++++++++++-------------
 pkg/R/initSmallEM.R                  |   2 +-
 pkg/R/main.R                         |  21 ++--
 pkg/R/selectVariables.R              |  10 +-
 9 files changed, 136 insertions(+), 111 deletions(-)

diff --git a/pkg/DESCRIPTION b/pkg/DESCRIPTION
index c975ee9..e18fddb 100644
--- a/pkg/DESCRIPTION
+++ b/pkg/DESCRIPTION
@@ -28,7 +28,7 @@ URL: http://git.auder.net/?p=valse.git
 License: MIT + file LICENSE
 RoxygenNote: 5.0.1
 Collate:
-    'plot.R'
+    'plot_valse.R'
     'main.R'
     'selectVariables.R'
     'constructionModelesLassoRank.R'
@@ -39,4 +39,3 @@ Collate:
     'EMGLLF.R'
     'generateXY.R'
     'A_NAMESPACE.R'
-    'plot_valse.R'
diff --git a/pkg/R/A_NAMESPACE.R b/pkg/R/A_NAMESPACE.R
index 81e91ec..8e1783e 100644
--- a/pkg/R/A_NAMESPACE.R
+++ b/pkg/R/A_NAMESPACE.R
@@ -7,7 +7,7 @@
 #' @include constructionModelesLassoRank.R
 #' @include selectVariables.R
 #' @include main.R
-#' @include plot.R
+#' @include plot_valse.R
 #'
 #' @useDynLib valse
 #'
diff --git a/pkg/R/EMGLLF.R b/pkg/R/EMGLLF.R
index 5a69a52..13a08da 100644
--- a/pkg/R/EMGLLF.R
+++ b/pkg/R/EMGLLF.R
@@ -2,17 +2,17 @@
 #'
 #' Description de EMGLLF
 #'
-#' @param phiInit Parametre initial de moyenne renormalisé
-#' @param rhoInit Parametre initial de variance renormalisé
-#' @param piInit Parametre initial des proportions
-#' @param gamInit Paramètre initial des probabilités a posteriori de chaque échantillon
-#' @param mini Nombre minimal d'itérations dans l'algorithme EM
-#' @param maxi Nombre maximal d'itérations dans l'algorithme EM
-#' @param gamma Puissance des proportions dans la pénalisation pour un Lasso adaptatif
-#' @param lambda Valeur du paramètre de régularisation du Lasso
-#' @param X Régresseurs
-#' @param Y Réponse
-#' @param tau Seuil pour accepter la convergence
+#' @param phiInit an initialization for phi
+#' @param rhoInit an initialization for rho
+#' @param piInit an initialization for pi
+#' @param gamInit initialization for the a posteriori probabilities
+#' @param mini integer, minimum number of iterations in the EM algorithm, by default = 10
+#' @param maxi integer, maximum number of iterations in the EM algorithm, by default = 100
+#' @param gamma integer for the power in the penaly, by default = 1
+#' @param lambda regularization parameter in the Lasso estimation
+#' @param X matrix of covariates (of size n*p)
+#' @param Y matrix of responses (of size n*m)
+#' @param eps real, threshold to say the EM algorithm converges, by default = 1e-4
 #'
 #' @return A list ... phi,rho,pi,LLF,S,affec:
 #'   phi : parametre de moyenne renormalisé, calculé par l'EM
@@ -23,7 +23,7 @@
 #'
 #' @export
 EMGLLF <- function(phiInit, rhoInit, piInit, gamInit,
-                   mini, maxi, gamma, lambda, X, Y, tau, fast=TRUE)
+                   mini, maxi, gamma, lambda, X, Y, eps, fast=TRUE)
 {
   if (!fast)
   {
diff --git a/pkg/R/EMGrank.R b/pkg/R/EMGrank.R
index 0e68cb4..7c0d91f 100644
--- a/pkg/R/EMGrank.R
+++ b/pkg/R/EMGrank.R
@@ -2,7 +2,6 @@
 #'
 #' Description de EMGrank
 #'
-#' @param phiInit ...
 #' @param Pi Parametre de proportion
 #' @param Rho Parametre initial de variance renormalisé
 #' @param mini Nombre minimal d'itérations dans l'algorithme EM
@@ -49,6 +48,7 @@ matricize <- function(X)
 # R version - slow but easy to read
 .EMGrank_R = function(Pi, Rho, mini, maxi, X, Y, tau, rank)
 {
+  require(MASS)
   #matrix dimensions
   n = dim(X)[1]
   p = dim(X)[2]
@@ -70,10 +70,10 @@ matricize <- function(X)
   ite = 1
   while (ite<=mini || (ite<=maxi && sumDeltaPhi>tau))
 	{
-    #M step: Mise à jour de Beta (et donc phi)
+    #M step: update for Beta ( and then phi)
     for(r in 1:k)
 		{
-      Z_indice = seq_len(n)[Z==r] #indices où Z == r
+      Z_indice = seq_len(n)[Z==r] #indices where Z == r
       if (length(Z_indice) == 0)
         next
       #U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
@@ -87,7 +87,7 @@ matricize <- function(X)
       phi[,,r] = s$u %*% diag(S) %*% t(s$v) %*% Rho[,,r]
     }
 
-		#Etape E et calcul de LLF
+		#Step E and computation of the loglikelihood
 		sumLogLLF2 = 0
 		for(i in seq_len(n))
 		{
diff --git a/pkg/R/constructionModelesLassoMLE.R b/pkg/R/constructionModelesLassoMLE.R
index ac54319..ba6f125 100644
--- a/pkg/R/constructionModelesLassoMLE.R
+++ b/pkg/R/constructionModelesLassoMLE.R
@@ -2,21 +2,33 @@
 #'
 #' Construct a collection of models with the Lasso-MLE procedure.
 #' 
+#' @param phiInit an initialization for phi, get by initSmallEM.R
+#' @param rhoInit an initialization for rho, get by initSmallEM.R
+#' @param piInit an initialization for pi, get by initSmallEM.R
+#' @param gamInit an initialization for gam, get by initSmallEM.R
+#' @param mini integer, minimum number of iterations in the EM algorithm, by default = 10
+#' @param maxi integer, maximum number of iterations in the EM algorithm, by default = 100
+#' @param gamma integer for the power in the penaly, by default = 1
+#' @param X matrix of covariates (of size n*p)
+#' @param Y matrix of responses (of size n*m)
+#' @param eps real, threshold to say the EM algorithm converges, by default = 1e-4
+#' @param S output of selectVariables.R
+#' @param ncores Number of cores, by default = 3
+#' @param fast TRUE to use compiled C code, FALSE for R code only
+#' @param verbose TRUE to show some execution traces
+#' 
+#' @return a list with several models, defined by phi, rho, pi, llh
 #'
-#' @param ...
-#'
-#' @return ...
-#'
-#' export
-constructionModelesLassoMLE = function(phiInit, rhoInit, piInit, gamInit, mini, maxi,
-	gamma, X, Y, thresh, tau, S, ncores=3, fast=TRUE, verbose=FALSE)
+#' @export
+constructionModelesLassoMLE = function( phiInit, rhoInit, piInit, gamInit, mini, maxi,gamma, X, Y,
+	 eps, S, ncores=3, fast=TRUE, verbose=FALSE)
 {
 	if (ncores > 1)
 	{
 		cl = parallel::makeCluster(ncores, outfile='')
 		parallel::clusterExport( cl, envir=environment(),
-			varlist=c("phiInit","rhoInit","gamInit","mini","maxi","gamma","X","Y","thresh",
-			"tau","S","ncores","verbose") )
+			varlist=c("phiInit","rhoInit","gamInit","mini","maxi","gamma","X","Y","eps",
+			"S","ncores","fast","verbose") )
 	}
 
 	# Individual model computation
@@ -40,7 +52,7 @@ constructionModelesLassoMLE = function(phiInit, rhoInit, piInit, gamInit, mini,
 
 		# lambda == 0 because we compute the EMV: no penalization here
 		res = EMGLLF(phiInit[col.sel,,],rhoInit,piInit,gamInit,mini,maxi,gamma,0,
-			X[,col.sel], Y, tau, fast)
+			X[,col.sel], Y, eps, fast)
 		
 		# Eval dimension from the result + selected
 		phiLambda2 = res$phi
diff --git a/pkg/R/constructionModelesLassoRank.R b/pkg/R/constructionModelesLassoRank.R
index 339ba60..5da26e3 100644
--- a/pkg/R/constructionModelesLassoRank.R
+++ b/pkg/R/constructionModelesLassoRank.R
@@ -1,84 +1,95 @@
 #' constructionModelesLassoRank
 #'
-#' TODO: description
+#' Construct a collection of models with the Lasso-Rank procedure.
+#' 
+#' @param S output of selectVariables.R
+#' @param k number of components
+#' @param mini integer, minimum number of iterations in the EM algorithm, by default = 10
+#' @param maxi integer, maximum number of iterations in the EM algorithm, by default = 100
+#' @param X matrix of covariates (of size n*p)
+#' @param Y matrix of responses (of size n*m)
+#' @param eps real, threshold to say the EM algorithm converges, by default = 1e-4
+#' @param rank.min integer, minimum rank in the low rank procedure, by default = 1
+#' @param rank.max integer, maximum rank in the low rank procedure, by default = 5
+#' @param ncores Number of cores, by default = 3
+#' @param fast TRUE to use compiled C code, FALSE for R code only
+#' @param verbose TRUE to show some execution traces
+#' 
+#' @return a list with several models, defined by phi, rho, pi, llh
 #'
-#' @param ...
-#'
-#' @return ...
-#'
-#' export
-constructionModelesLassoRank = function(pi, rho, mini, maxi, X, Y, tau, A1, rangmin,
-	rangmax, ncores, fast=TRUE, verbose=FALSE)
+#' @export
+constructionModelesLassoRank = function(S, k, mini, maxi, X, Y, eps, rank.min,
+                                        rank.max, ncores, fast=TRUE, verbose=FALSE)
 {
   n = dim(X)[1]
   p = dim(X)[2]
-  m = dim(rho)[2]
-  k = dim(rho)[3]
-  L = dim(A1)[2]
-
-	# On cherche les rangs possiblement intéressants
-  deltaRank = rangmax - rangmin + 1
+  m = dim(Y)[2]
+  L = length(S)
+  
+  # Possible interesting ranks
+  deltaRank = rank.max - rank.min + 1
   Size = deltaRank^k
-  Rank = matrix(0, nrow=Size, ncol=k)
+  RankLambda = matrix(0, nrow=Size*L, ncol=k+1)
   for (r in 1:k)
-	{
-		# On veut le tableau de toutes les combinaisons de rangs possibles
-		# Dans la première colonne : on répète (rangmax-rangmin)^(k-1) chaque chiffre :
-		#   ça remplit la colonne
-		# Dans la deuxieme : on répète (rangmax-rangmin)^(k-2) chaque chiffre,
-		#   et on fait ça (rangmax-rangmin)^2 fois
-		# ...
-		# Dans la dernière, on répète chaque chiffre une fois,
-		#   et on fait ça (rangmin-rangmax)^(k-1) fois.
-    Rank[,r] = rangmin + rep(0:(deltaRank-1), deltaRank^(r-1), each=deltaRank^(k-r))
+  {
+    # On veut le tableau de toutes les combinaisons de rangs possibles, et des lambdas
+    # Dans la première colonne : on répète (rank.max-rank.min)^(k-1) chaque chiffre :
+    #   ça remplit la colonne
+    # Dans la deuxieme : on répète (rank.max-rank.min)^(k-2) chaque chiffre,
+    #   et on fait ça (rank.max-rank.min)^2 fois
+    # ...
+    # Dans la dernière, on répète chaque chiffre une fois,
+    #   et on fait ça (rank.min-rank.max)^(k-1) fois.
+    RankLambda[,r] = rep(rank.min + rep(0:(deltaRank-1), deltaRank^(r-1), each=deltaRank^(k-r)), each = L)
   }
-
+  RankLambda[,k+1] = rep(1:L, times = Size)
+  
   if (ncores > 1)
-	{
+  {
     cl = parallel::makeCluster(ncores, outfile='')
     parallel::clusterExport( cl, envir=environment(),
-			varlist=c("A1","Size","Pi","Rho","mini","maxi","X","Y","tau",
-			"Rank","m","phi","ncores","verbose") )
-	}
-
-	computeAtLambda <- function(lambdaIndex)
-	{
-		if (ncores > 1)
-			require("valse") #workers start with an empty environment
-
-    # on ne garde que les colonnes actives
-    # 'active' sera l'ensemble des variables informatives
-    active = A1[,lambdaIndex]
-    active = active[-(active==0)]
-		phi = array(0, dim=c(p,m,k,Size))
-		llh = matrix(0, Size, 2) #log-likelihood
-    if (length(active) > 0)
-		{
-      for (j in 1:Size)
-			{
-        res = EMGrank(Pi[,lambdaIndex], Rho[,,,lambdaIndex], mini, maxi,
-					X[,active], Y, tau, Rank[j,], fast)
-        llh = rbind(llh,
-					c( res$LLF, sum(Rank[j,] * (length(active)- Rank[j,] + m)) ) )
-        phi[active,,,] = rbind(phi[active,,,], res$phi)
+                             varlist=c("A1","Size","Pi","Rho","mini","maxi","X","Y","eps",
+                                       "Rank","m","phi","ncores","verbose") )
+  }
+  
+  computeAtLambda <- function(index)
+  {
+    lambdaIndex = RankLambda[index,k+1]
+    rankIndex = RankLambda[index,1:k]
+    if (ncores > 1)
+      require("valse") #workers start with an empty environment
+    
+    # 'relevant' will be the set of relevant columns
+    selected = S[[lambdaIndex]]$selected
+    relevant = c()
+    for (j in 1:p){
+      if (length(selected[[j]])>0){
+        relevant = c(relevant,j)
       }
     }
-		list("llh"=llh, "phi"=phi)
-	}
-
-	#Pour chaque lambda de la grille, on calcule les coefficients
+    if (max(rankIndex)<length(relevant)){
+      phi = array(0, dim=c(p,m,k))
+      if (length(relevant) > 0)
+      {
+        res = EMGrank(S[[lambdaIndex]]$Pi, S[[lambdaIndex]]$Rho, mini, maxi,
+                      X[,relevant], Y, eps, rankIndex, fast)
+        llh = c( res$LLF, sum(rankIndex * (length(relevant)- rankIndex + m)) ) 
+        phi[relevant,,] = res$phi
+      }
+      list("llh"=llh, "phi"=phi, "pi" = S[[lambdaIndex]]$Pi, "rho" = S[[lambdaIndex]]$Rho)
+      
+    }
+  }
+  
+  #For each lambda in the grid we compute the estimators
   out =
-		if (ncores > 1)
-			parLapply(cl, seq_along(glambda), computeAtLambda)
-		else
-			lapply(seq_along(glambda), computeAtLambda)
-
-	if (ncores > 1)
+    if (ncores > 1)
+      parLapply(cl, seq_len(length(S)*Size), computeAtLambda)
+  else
+    lapply(seq_len(length(S)*Size), computeAtLambda)
+  
+  if (ncores > 1)
     parallel::stopCluster(cl)
-
-	# TODO: this is a bit ugly. Better use bigmemory and fill llh/phi in-place
-	# (but this also adds a dependency...)
-	llh <- do.call( rbind, lapply(out, function(model) model$llh) )
-	phi <- do.call( rbind, lapply(out, function(model) model$phi) )
-	list("llh"=llh, "phi"=phi)
+  
+  out
 }
diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R
index 9b58a0c..5dcafb8 100644
--- a/pkg/R/initSmallEM.R
+++ b/pkg/R/initSmallEM.R
@@ -63,7 +63,7 @@ initSmallEM = function(k,X,Y, fast=TRUE)
 		maxiInit = 11
 		
 		new_EMG = EMGLLF(phiInit1[,,,repet], rhoInit1[,,,repet], piInit1[repet,],
-			gamInit1[,,repet], miniInit, maxiInit, gamma=1, lambda=0, X, Y, tau=1e-4, fast)
+			gamInit1[,,repet], miniInit, maxiInit, gamma=1, lambda=0, X, Y, eps=1e-4, fast)
 		LLFEessai = new_EMG$LLF
 		LLFinit1[repet] = LLFEessai[length(LLFEessai)]
 	}
diff --git a/pkg/R/main.R b/pkg/R/main.R
index 6d315cd..634c273 100644
--- a/pkg/R/main.R
+++ b/pkg/R/main.R
@@ -12,10 +12,11 @@
 #' @param eps real, threshold to say the EM algorithm converges, by default = 1e-4
 #' @param kmin integer, minimum number of clusters, by default = 2
 #' @param kmax integer, maximum number of clusters, by default = 10
-#' @param rang.min integer, minimum rank in the low rank procedure, by default = 1
-#' @param rang.max integer, maximum rank in the
+#' @param rank.min integer, minimum rank in the low rank procedure, by default = 1
+#' @param rank.max integer, maximum rank in the low rank procedure, by default = 5
 #' @param ncores_outer Number of cores for the outer loop on k
 #' @param ncores_inner Number of cores for the inner loop on lambda
+#' @param thresh real, threshold to say a variable is relevant, by default = 1e-8
 #' @param size_coll_mod (Maximum) size of a collection of models
 #' @param fast TRUE to use compiled C code, FALSE for R code only
 #' @param verbose TRUE to show some execution traces
@@ -26,7 +27,8 @@
 #' #TODO: a few examples
 #' @export
 valse = function(X, Y, procedure='LassoMLE', selecMod='DDSE', gamma=1, mini=10, maxi=50,
-                 eps=1e-4, kmin=2, kmax=2, rank.min=1, rank.max=10, ncores_outer=1, ncores_inner=1,
+                 eps=1e-4, kmin=2, kmax=3, rank.min=1, rank.max=5, ncores_outer=1, ncores_inner=1,
+                 thresh=1e-8,
                  size_coll_mod=10, fast=TRUE, verbose=FALSE, plot = TRUE)
 {
   p = dim(X)[2]
@@ -40,8 +42,8 @@ valse = function(X, Y, procedure='LassoMLE', selecMod='DDSE', gamma=1, mini=10,
   {
     cl = parallel::makeCluster(ncores_outer, outfile='')
     parallel::clusterExport( cl=cl, envir=environment(), varlist=c("X","Y","procedure",
-                                                                   "selecMod","gamma","mini","maxi","eps","kmin","kmax","rang.min","rang.max",
-                                                                   "ncores_outer","ncores_inner","verbose","p","m") )
+                                                                   "selecMod","gamma","mini","maxi","eps","kmin","kmax","rank.min","rank.max",
+                                                                   "ncores_outer","ncores_inner","thresh","size_coll_mod","verbose","p","m") )
   }
   
   # Compute models with k components
@@ -66,7 +68,7 @@ valse = function(X, Y, procedure='LassoMLE', selecMod='DDSE', gamma=1, mini=10,
     #select variables according to each regularization parameter
     #from the grid: S$selected corresponding to selected variables
     S = selectVariables(P$phiInit, P$rhoInit, P$piInit, P$gamInit, mini, maxi, gamma,
-                        grid_lambda, X, Y, 1e-8, eps, ncores_inner, fast) #TODO: 1e-8 as arg?! eps?
+                        grid_lambda, X, Y, thresh, eps, ncores_inner, fast) 
     
     if (procedure == 'LassoMLE')
     {
@@ -74,8 +76,9 @@ valse = function(X, Y, procedure='LassoMLE', selecMod='DDSE', gamma=1, mini=10,
         print('run the procedure Lasso-MLE')
       #compute parameter estimations, with the Maximum Likelihood
       #Estimator, restricted on selected variables.
-      models <- constructionModelesLassoMLE(P$phiInit, P$rhoInit, P$piInit, P$gamInit,
-                                            mini, maxi, gamma, X, Y, thresh, eps, S, ncores_inner, fast, verbose)
+      models <- constructionModelesLassoMLE( P$phiInit, P$rhoInit, P$piInit, P$gamInit, 
+                                            mini, maxi, gamma, X, Y, eps, S, ncores_inner, fast, verbose)
+      
     }
     else
     {
@@ -83,7 +86,7 @@ valse = function(X, Y, procedure='LassoMLE', selecMod='DDSE', gamma=1, mini=10,
         print('run the procedure Lasso-Rank')
       #compute parameter estimations, with the Low Rank
       #Estimator, restricted on selected variables.
-      models <- constructionModelesLassoRank(S$Pi, S$Rho, mini, maxi, X, Y, eps, S,
+      models <- constructionModelesLassoRank(S, k, mini, maxi, X, Y, eps,
                                              rank.min, rank.max, ncores_inner, fast, verbose)
     }
     #warning! Some models are NULL after running selectVariables
diff --git a/pkg/R/selectVariables.R b/pkg/R/selectVariables.R
index 65fbde5..0225287 100644
--- a/pkg/R/selectVariables.R
+++ b/pkg/R/selectVariables.R
@@ -12,8 +12,8 @@
 #' @param glambda grid of regularization parameters
 #' @param X			 matrix of regressors
 #' @param Y			 matrix of responses
-#' @param thres	 threshold to consider a coefficient to be equal to 0
-#' @param tau		 threshold to say that EM algorithm has converged
+#' @param thresh real, threshold to say a variable is relevant, by default = 1e-8
+#' @param eps		 threshold to say that EM algorithm has converged
 #' @param ncores Number or cores for parallel execution (1 to disable)
 #'
 #' @return a list of outputs, for each lambda in grid: selected,Rho,Pi
@@ -23,20 +23,20 @@
 #' @export
 #'
 selectVariables = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,
-                           X,Y,thresh,tau, ncores=3, fast=TRUE)
+                           X,Y,thresh=1e-8,eps, ncores=3, fast=TRUE)
 {
   if (ncores > 1)
   {
     cl = parallel::makeCluster(ncores, outfile='')
     parallel::clusterExport(cl=cl,
-                            varlist=c("phiInit","rhoInit","gamInit","mini","maxi","glambda","X","Y","thresh","tau"),
+                            varlist=c("phiInit","rhoInit","gamInit","mini","maxi","glambda","X","Y","thresh","eps"),
                             envir=environment())
   }
   
   # Computation for a fixed lambda
   computeCoefs <- function(lambda)
   {
-    params = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau,fast)
+    params = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,eps,fast)
     
     p = dim(phiInit)[1]
     m = dim(phiInit)[2]
-- 
2.44.0