From 4d9db27f0d1749e5577038dedbc5f4d0826f2772 Mon Sep 17 00:00:00 2001
From: emilie <emilie@devijver.org>
Date: Thu, 13 Apr 2017 16:24:26 +0200
Subject: [PATCH] fix initialization and made some update

---
 pkg/R/computeGridLambda.R |  2 +-
 pkg/R/initSmallEM.R       | 27 ++++++++---------
 pkg/R/main.R              | 63 ++++++++++++++++++++-------------------
 3 files changed, 47 insertions(+), 45 deletions(-)

diff --git a/pkg/R/computeGridLambda.R b/pkg/R/computeGridLambda.R
index 8ec4d66..9d06aed 100644
--- a/pkg/R/computeGridLambda.R
+++ b/pkg/R/computeGridLambda.R
@@ -33,5 +33,5 @@ computeGridLambda = function(phiInit, rhoInit, piInit, gamInit, X, Y,
 			grid[i,j,] = abs(list_EMG$S[i,j,]) / (n*list_EMG$pi^gamma)
 	}
 	grid = unique(grid)
-	grid
+	sort(grid)
 }
diff --git a/pkg/R/initSmallEM.R b/pkg/R/initSmallEM.R
index 5dcafb8..bd5ce17 100644
--- a/pkg/R/initSmallEM.R
+++ b/pkg/R/initSmallEM.R
@@ -13,21 +13,21 @@ initSmallEM = function(k,X,Y, fast=TRUE)
 	n = nrow(Y)
 	m = ncol(Y)
 	p = ncol(X)
-  
-	Zinit1 = array(0, dim=c(n,20))
-	betaInit1 = array(0, dim=c(p,m,k,20))
-	sigmaInit1 = array(0, dim = c(m,m,k,20))
-	phiInit1 = array(0, dim = c(p,m,k,20))
-	rhoInit1 = array(0, dim = c(m,m,k,20))
+  nIte = 20
+	Zinit1 = array(0, dim=c(n,nIte))
+	betaInit1 = array(0, dim=c(p,m,k,nIte))
+	sigmaInit1 = array(0, dim = c(m,m,k,nIte))
+	phiInit1 = array(0, dim = c(p,m,k,nIte))
+	rhoInit1 = array(0, dim = c(m,m,k,nIte))
 	Gam = matrix(0, n, k)
-	piInit1 = matrix(0,20,k)
-	gamInit1 = array(0, dim=c(n,k,20))
+	piInit1 = matrix(0,nIte,k)
+	gamInit1 = array(0, dim=c(n,k,nIte))
 	LLFinit1 = list()
 
 	#require(MASS) #Moore-Penrose generalized inverse of matrix
-	for(repet in 1:20)
+	for(repet in 1:nIte)
 	{
-	  distance_clus = dist(X)
+	  distance_clus = dist(cbind(X,Y))
 	  tree_hier = hclust(distance_clus)
 	  Zinit1[,repet] = cutree(tree_hier, k)
 
@@ -62,13 +62,12 @@ initSmallEM = function(k,X,Y, fast=TRUE)
 		miniInit = 10
 		maxiInit = 11
 		
-		new_EMG = EMGLLF(phiInit1[,,,repet], rhoInit1[,,,repet], piInit1[repet,],
+		init_EMG = EMGLLF(phiInit1[,,,repet], rhoInit1[,,,repet], piInit1[repet,],
 			gamInit1[,,repet], miniInit, maxiInit, gamma=1, lambda=0, X, Y, eps=1e-4, fast)
-		LLFEessai = new_EMG$LLF
+		LLFEessai = init_EMG$LLF
 		LLFinit1[repet] = LLFEessai[length(LLFEessai)]
 	}
-
-	b = which.max(LLFinit1)
+	b = which.min(LLFinit1)
 	phiInit = phiInit1[,,,b]
 	rhoInit = rhoInit1[,,,b]
 	piInit = piInit1[b,]
diff --git a/pkg/R/main.R b/pkg/R/main.R
index 634c273..238160c 100644
--- a/pkg/R/main.R
+++ b/pkg/R/main.R
@@ -123,36 +123,39 @@ valse = function(X, Y, procedure='LassoMLE', selecMod='DDSE', gamma=1, mini=10,
   
   print(tableauRecap)
   tableauRecap = tableauRecap[which(tableauRecap[,4]!= Inf),]
-  modSel = capushe::capushe(tableauRecap, n)
-  indModSel <-
-    if (selecMod == 'DDSE')
-      as.numeric(modSel@DDSE@model)
-  else if (selecMod == 'Djump')
-    as.numeric(modSel@Djump@model)
-  else if (selecMod == 'BIC')
-    modSel@BIC_capushe$model
-  else if (selecMod == 'AIC')
-    modSel@AIC_capushe$model
   
-  mod = as.character(tableauRecap[indModSel,1])
-  listMod = as.integer(unlist(strsplit(mod, "[.]")))
-  modelSel = models_list[[listMod[1]]][[listMod[2]]]
+  return(tableauRecap)
   
-  ##Affectations
-  Gam = matrix(0, ncol = length(modelSel$pi), nrow = n)
-  for (i in 1:n){
-    for (r in 1:length(modelSel$pi)){
-      sqNorm2 = sum( (Y[i,]%*%modelSel$rho[,,r]-X[i,]%*%modelSel$phi[,,r])^2 )
-      Gam[i,r] = modelSel$pi[r] * exp(-0.5*sqNorm2)* det(modelSel$rho[,,r])
-    }
-  }
-  Gam = Gam/rowSums(Gam)
-  modelSel$affec = apply(Gam, 1,which.max)
-  modelSel$proba = Gam
-  
-  if (plot){
-    print(plot_valse(X,Y,modelSel,n))
-  }
-  
-  return(modelSel)
+  # modSel = capushe::capushe(tableauRecap, n)
+  # indModSel <-
+  #   if (selecMod == 'DDSE')
+  #     as.numeric(modSel@DDSE@model)
+  # else if (selecMod == 'Djump')
+  #   as.numeric(modSel@Djump@model)
+  # else if (selecMod == 'BIC')
+  #   modSel@BIC_capushe$model
+  # else if (selecMod == 'AIC')
+  #   modSel@AIC_capushe$model
+  # 
+  # mod = as.character(tableauRecap[indModSel,1])
+  # listMod = as.integer(unlist(strsplit(mod, "[.]")))
+  # modelSel = models_list[[listMod[1]]][[listMod[2]]]
+  # 
+  # ##Affectations
+  # Gam = matrix(0, ncol = length(modelSel$pi), nrow = n)
+  # for (i in 1:n){
+  #   for (r in 1:length(modelSel$pi)){
+  #     sqNorm2 = sum( (Y[i,]%*%modelSel$rho[,,r]-X[i,]%*%modelSel$phi[,,r])^2 )
+  #     Gam[i,r] = modelSel$pi[r] * exp(-0.5*sqNorm2)* det(modelSel$rho[,,r])
+  #   }
+  # }
+  # Gam = Gam/rowSums(Gam)
+  # modelSel$affec = apply(Gam, 1,which.max)
+  # modelSel$proba = Gam
+  # 
+  # if (plot){
+  #   print(plot_valse(X,Y,modelSel,n))
+  # }
+  # 
+  # return(modelSel)
 }
-- 
2.44.0