From c3bc47052f3ccb659659c59a82e9a99ea842398d Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Sat, 11 Feb 2017 23:32:09 +0100
Subject: [PATCH] fix memory leaks on EMGLLF, test OK for EMGrank

---
 R/initSmallEM.R                               |   8 +-
 R/modelSelection.R                            |   2 +-
 R/vec_bin.R                                   |  23 ---
 src/adapters/a.constructionModelesLassoMLE.c  |  12 +-
 src/adapters/a.constructionModelesLassoRank.c |  12 +-
 src/sources/EMGrank.c                         |   4 +-
 src/sources/constructionModelesLassoMLE.c     |   6 +-
 src/sources/constructionModelesLassoMLE.h     |   2 +-
 src/sources/constructionModelesLassoRank.c    |   8 +-
 src/sources/constructionModelesLassoRank.h    |   2 +-
 src/sources/utils.h                           |   2 +-
 src/test/OUT                                  |  30 ----
 .../generateRunSaveTest_EMGLLF.R              |   4 +-
 .../generateRunSaveTest_EMGrank.R             |  52 ++++---
 ...eRunSaveTest_constructionModelesLassoMLE.R |  83 ++++++-----
 ...eRunSaveTest_constructionModelesLassoMLE.m |   4 +-
 ...RunSaveTest_constructionModelesLassoRank.R |  58 ++++----
 ...RunSaveTest_constructionModelesLassoRank.m |   4 +-
 .../generateRunSaveTest_selectiontotale.R     |  67 +++++----
 src/test/generate_test_data/helpers/EMGrank.R | 101 ++++++-------
 .../helpers/constructionModelesLassoMLE.R     |  10 +-
 .../helpers/constructionModelesLassoMLE.m     |   8 +-
 .../helpers/constructionModelesLassoRank.R    |   6 +-
 .../helpers/constructionModelesLassoRank.m    |   6 +-
 src/test/sourceAll.R                          |   6 +
 src/test/test.ConstructionModelesLassoMLE.c   | 135 +++++-------------
 src/test/test.EMGLLF.c                        |  46 +++---
 src/test/test.EMGrank.c                       |  32 ++---
 src/test/test.constructionModelesLassoRank.c  | 105 ++++----------
 src/test/test.selectiontotale.c               | 131 +++++------------
 src/test/test_utils.c                         |  21 +--
 src/test/test_utils.h                         |   4 +-
 32 files changed, 404 insertions(+), 590 deletions(-)
 delete mode 100644 R/vec_bin.R
 delete mode 100644 src/test/OUT
 create mode 100644 src/test/sourceAll.R

diff --git a/R/initSmallEM.R b/R/initSmallEM.R
index c836523..6dd7457 100644
--- a/R/initSmallEM.R
+++ b/R/initSmallEM.R
@@ -33,15 +33,13 @@ initSmallEM = function(k,X,Y,tau)
 		for(r in 1:k)
 		{
 			Z = Zinit1[,repet]
-			Z_bin = vec_bin(Z,r)
-			Z_vec = Z_bin$vec #vecteur 0 et 1 aux endroits o? Z==r
-			Z_indice = Z_bin$indice #renvoit les indices o? Z==r
+			Z_indice = seq_len(n)[Z == r] #renvoit les indices où Z==r
 			
-			betaInit1[,,r,repet] = ginv( crossprod(X[Z_indice,]) )   %*%   crossprod(X[Z_indice,], Y[Z_indice,]) 
+			betaInit1[,,r,repet] = ginv(crossprod(X[Z_indice,])) %*% crossprod(X[Z_indice,], Y[Z_indice,])
 			sigmaInit1[,,r,repet] = diag(m)
 			phiInit1[,,r,repet] = betaInit1[,,r,repet] #/ sigmaInit1[,,r,repet]
 			rhoInit1[,,r,repet] = solve(sigmaInit1[,,r,repet])
-			piInit1[repet,r] = sum(Z_vec)/n
+			piInit1[repet,r] = mean(Z == r)
 		}
 		
 		for(i in 1:n)
diff --git a/R/modelSelection.R b/R/modelSelection.R
index bc7eeae..81e832a 100644
--- a/R/modelSelection.R
+++ b/R/modelSelection.R
@@ -28,7 +28,7 @@ modelSelection = function(LLF)
 			}
 			b = max(a)
 			#indices[i] : first indices of the binary vector where u_i ==1
-			indices[i] = which.max(vec_bin(LLF,b)[[1]])
+			indices[i] = which.max(LLF == b)
 		}
 	}
 
diff --git a/R/vec_bin.R b/R/vec_bin.R
deleted file mode 100644
index 27f771b..0000000
--- a/R/vec_bin.R
+++ /dev/null
@@ -1,23 +0,0 @@
-#' A function needed in initSmallEM
-#'
-#' @param X vector with integer values
-#' @param r integer
-#'
-#' @return a list with Z (a binary vector of size the size of X) and indices where Z is equal to 1
-vec_bin = function(X,r)
-{
-	Z = rep(0,length(X))
-	indice = c()
-	j = 1
-	for (i in 1:length(X))
-	{
-		if(X[i] == r)
-		{
-			Z[i] = 1
-			indice[j] = i
-			j=j+1
-		} else
-			Z[i] = 0
-	}
-	return (list(Z=Z,indice=indice))
-}
diff --git a/src/adapters/a.constructionModelesLassoMLE.c b/src/adapters/a.constructionModelesLassoMLE.c
index 72c3173..ec519a9 100644
--- a/src/adapters/a.constructionModelesLassoMLE.c
+++ b/src/adapters/a.constructionModelesLassoMLE.c
@@ -52,7 +52,7 @@ SEXP constructionModelesLassoMLE(
 	// OUTPUTS //
 	/////////////
 
-	SEXP phi, rho, pi, lvraisemblance, dimPhi, dimRho;
+	SEXP phi, rho, pi, llh, dimPhi, dimRho;
 	PROTECT(dimPhi = allocVector(INTSXP, 4));
 	int* pDimPhi = INTEGER(dimPhi);
 	pDimPhi[0] = p; pDimPhi[1] = m; pDimPhi[2] = k; pDimPhi[3] = L;
@@ -62,8 +62,8 @@ SEXP constructionModelesLassoMLE(
 	PROTECT(phi = allocArray(REALSXP, dimPhi));
 	PROTECT(rho = allocArray(REALSXP, dimRho));
 	PROTECT(pi = allocMatrix(REALSXP, k, L));
-	PROTECT(lvraisemblance = allocMatrix(REALSXP, L, 2));
-	double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pLvraisemblance=REAL(lvraisemblance);
+	PROTECT(llh = allocMatrix(REALSXP, L, 2));
+	double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pllh=REAL(llh);
 
 	/////////////////////////////////////////
 	// Call to constructionModelesLassoMLE //
@@ -71,13 +71,13 @@ SEXP constructionModelesLassoMLE(
 
 	constructionModelesLassoMLE_core(
 		phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2,
-		pPhi,pRho,pPi,pLvraisemblance,
+		pPhi,pRho,pPi,pllh,
 		n,p,m,k,L);
 
 	// Build list from OUT params and return it
 	SEXP listParams, listNames;
 	PROTECT(listParams = allocVector(VECSXP, 4));
-	char* lnames[4] = {"phi", "rho", "pi", "lvraisemblance"}; //lists labels
+	char* lnames[4] = {"phi", "rho", "pi", "llh"}; //lists labels
 	PROTECT(listNames = allocVector(STRSXP,4));
 	for (int i=0; i<4; i++)
 		SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
@@ -85,7 +85,7 @@ SEXP constructionModelesLassoMLE(
 	SET_VECTOR_ELT(listParams, 0, phi);
 	SET_VECTOR_ELT(listParams, 1, rho);
 	SET_VECTOR_ELT(listParams, 2, pi);
-	SET_VECTOR_ELT(listParams, 3, lvraisemblance);
+	SET_VECTOR_ELT(listParams, 3, llh);
 
 	UNPROTECT(8);
 	return listParams;
diff --git a/src/adapters/a.constructionModelesLassoRank.c b/src/adapters/a.constructionModelesLassoRank.c
index 4833e19..0e069d4 100644
--- a/src/adapters/a.constructionModelesLassoRank.c
+++ b/src/adapters/a.constructionModelesLassoRank.c
@@ -46,13 +46,13 @@ SEXP constructionModelesLassoRank(
 	/////////////
 
 	int Size = pow(rangmax-rangmin+1,k);
-	SEXP phi, lvraisemblance, dimPhi;
+	SEXP phi, llh, dimPhi;
 	PROTECT(dimPhi = allocVector(INTSXP, 4));
 	int* pDimPhi = INTEGER(dimPhi);
 	pDimPhi[0] = p; pDimPhi[1] = m; pDimPhi[2] = k; pDimPhi[3] = L*Size;
 	PROTECT(phi = allocArray(REALSXP, dimPhi));
-	PROTECT(lvraisemblance = allocMatrix(REALSXP, L*Size, 2));
-	double *pPhi=REAL(phi), *pLvraisemblance=REAL(lvraisemblance);
+	PROTECT(llh = allocMatrix(REALSXP, L*Size, 2));
+	double *pPhi=REAL(phi), *pllh=REAL(llh);
 
 	//////////////////////////////////////////
 	// Call to constructionModelesLassoRank //
@@ -60,19 +60,19 @@ SEXP constructionModelesLassoRank(
 
 	constructionModelesLassoRank_core(
 		Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax,
-		pPhi,pLvraisemblance,
+		pPhi,pllh,
 		n,p,m,k,L);
 
 	// Build list from OUT params and return it
 	SEXP listParams, listNames;
 	PROTECT(listParams = allocVector(VECSXP, 2));
-	char* lnames[2] = {"phi", "lvraisemblance"}; //lists labels
+	char* lnames[2] = {"phi", "llh"}; //lists labels
 	PROTECT(listNames = allocVector(STRSXP,2));
 	for (int i=0; i<2; i++)
 		SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
 	setAttrib(listParams, R_NamesSymbol, listNames);
 	SET_VECTOR_ELT(listParams, 0, phi);
-	SET_VECTOR_ELT(listParams, 1, lvraisemblance);
+	SET_VECTOR_ELT(listParams, 1, llh);
 
 	UNPROTECT(5);
 	return listParams;
diff --git a/src/sources/EMGrank.c b/src/sources/EMGrank.c
index 2422fc0..3a9bf94 100644
--- a/src/sources/EMGrank.c
+++ b/src/sources/EMGrank.c
@@ -132,7 +132,7 @@ void EMGrank_core(
 				{
 					Real dotProduct = 0.0;
 					for (int u=0; u<cardClustR; u++)
-						dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,j,n,m)];
+						dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,jj,n,m)];
 					tXrYr[mi(j,jj,p,m)] = dotProduct;
 				}
 			}
@@ -281,7 +281,7 @@ void EMGrank_core(
 			for (int jj=0; jj<p; jj++)
 			{
 				for (int r=0; r<k; r++)
-					Phi[ai(j,jj,r,p,m,k)] = phi[ai(j,jj,r,p,m,k)];
+					Phi[ai(jj,j,r,p,m,k)] = phi[ai(jj,j,r,p,m,k)];
 			}
 		}
 		ite++;
diff --git a/src/sources/constructionModelesLassoMLE.c b/src/sources/constructionModelesLassoMLE.c
index 6b92094..ad2a718 100644
--- a/src/sources/constructionModelesLassoMLE.c
+++ b/src/sources/constructionModelesLassoMLE.c
@@ -25,7 +25,7 @@ void constructionModelesLassoMLE_core(
 	Real* phi,// estimateur ainsi calculé par le Lasso
 	Real* rho,// estimateur ainsi calculé par le Lasso
 	Real* pi, // estimateur ainsi calculé par le Lasso
-	Real* lvraisemblance, // estimateur ainsi calculé par le Lasso
+	Real* llh, // estimateur ainsi calculé par le Lasso
 	// additional size parameters
 	int n, // taille de l'echantillon
 	int p, // nombre de covariables
@@ -204,8 +204,8 @@ void constructionModelesLassoMLE_core(
 			}
 			sumLogDensit += log(densite[lambdaIndex*n+i]);
 		}
-		lvraisemblance[mi(lambdaIndex,0,L,2)] = sumLogDensit;
-		lvraisemblance[mi(lambdaIndex,1,L,2)] = (dimension+m+1)*k-1;
+		llh[mi(lambdaIndex,0,L,2)] = sumLogDensit;
+		llh[mi(lambdaIndex,1,L,2)] = (dimension+m+1)*k-1;
 
 		free(a);
 		free(YiRhoR);
diff --git a/src/sources/constructionModelesLassoMLE.h b/src/sources/constructionModelesLassoMLE.h
index 992c2c7..06e500d 100644
--- a/src/sources/constructionModelesLassoMLE.h
+++ b/src/sources/constructionModelesLassoMLE.h
@@ -23,7 +23,7 @@ void constructionModelesLassoMLE_core(
 	Real* phi,
 	Real* rho,
 	Real* pi,
-	Real* lvraisemblance,
+	Real* llh,
 	// additional size parameters
 	int n,
 	int p,
diff --git a/src/sources/constructionModelesLassoRank.c b/src/sources/constructionModelesLassoRank.c
index 8eee0eb..2be1982 100644
--- a/src/sources/constructionModelesLassoRank.c
+++ b/src/sources/constructionModelesLassoRank.c
@@ -19,7 +19,7 @@ void constructionModelesLassoRank_core(
 	int rangmax,	//rang maximum autorisé
 	// OUT parameters (all pointers, to be modified)
 	Real* phi,// estimateur ainsi calculé par le Lasso
-	Real* lvraisemblance,// estimateur ainsi calculé par le Lasso
+	Real* llh,// estimateur ainsi calculé par le Lasso
 	// additional size parameters
 	int n,// taille de l'echantillon
 	int p,// nombre de covariables
@@ -106,13 +106,13 @@ for (int r=0; r<k; r++)
 			free(Xactive);
 			free(PiLambda);
 			free(RhoLambda);
-			//lvraisemblance((lambdaIndex-1)*Size+j,:) = [LLF, dot(Rank(j,:), length(active)-Rank(j,:)+m)];
-			lvraisemblance[mi(lambdaIndex*Size+j,0,L*Size,2)] = LLF;
+			//llh((lambdaIndex-1)*Size+j,:) = [LLF, dot(Rank(j,:), length(active)-Rank(j,:)+m)];
+			llh[mi(lambdaIndex*Size+j,0,L*Size,2)] = LLF;
 			//dot(Rank(j,:), length(active)-Rank(j,:)+m)
 			Real dotProduct = 0.0;
 			for (int r=0; r<k; r++)
 				dotProduct += Rank[mi(j,r,Size,k)] * (longueurActive-Rank[mi(j,r,Size,k)]+m);
-			lvraisemblance[mi(lambdaIndex*Size+j,1,Size*L,2)] = dotProduct;
+			llh[mi(lambdaIndex*Size+j,1,Size*L,2)] = dotProduct;
 			//phi(active,:,:,(lambdaIndex-1)*Size+j) = phiLambda;
 			for (int jj=0; jj<longueurActive; jj++)
 			{
diff --git a/src/sources/constructionModelesLassoRank.h b/src/sources/constructionModelesLassoRank.h
index 91a71d8..dbe8f59 100644
--- a/src/sources/constructionModelesLassoRank.h
+++ b/src/sources/constructionModelesLassoRank.h
@@ -18,7 +18,7 @@ void constructionModelesLassoRank_core(
 	int rangmax,
 	// OUT parameters
 	Real* phi,
-	Real* lvraisemblance,
+	Real* llh,
 	// additional size parameters
 	int n,
 	int p,
diff --git a/src/sources/utils.h b/src/sources/utils.h
index b33c69e..f8cf59e 100644
--- a/src/sources/utils.h
+++ b/src/sources/utils.h
@@ -7,7 +7,7 @@
  * Types
  *******/
 
-typedef float Real;
+typedef double Real;
 //typedef uint32_t UInt;
 //typedef int32_t Int;
 
diff --git a/src/test/OUT b/src/test/OUT
deleted file mode 100644
index 873fe4a..0000000
--- a/src/test/OUT
+++ /dev/null
@@ -1,30 +0,0 @@
-==12027== Memcheck, a memory error detector
-==12027== Copyright (C) 2002-2015, and GNU GPL'd, by Julian Seward et al.
-==12027== Using Valgrind-3.12.0 and LibVEX; rerun with -h for copyright info
-==12027== Command: ./test.EMGLLF
-==12027== 
-Checking phi
-    Inaccuracy: max(abs(error)) = 0.443413 >= 1e-05
-Checking rho
-    Inaccuracy: max(abs(error)) = 1.51588 >= 1e-05
-Checking pi
-    Inaccuracy: max(abs(error)) = 0.0839327 >= 1e-05
-Checking LLF
-    Inaccuracy: max(abs(error)) = 1.02552 >= 1e-05
-Checking S
-    Inaccuracy: max(abs(error)) = 90.3362 >= 1e-05
-==12027== 
-==12027== HEAP SUMMARY:
-==12027==     in use at exit: 523 bytes in 18 blocks
-==12027==   total heap usage: 168 allocs, 150 frees, 348,490 bytes allocated
-==12027== 
-==12027== LEAK SUMMARY:
-==12027==    definitely lost: 515 bytes in 17 blocks
-==12027==    indirectly lost: 0 bytes in 0 blocks
-==12027==      possibly lost: 0 bytes in 0 blocks
-==12027==    still reachable: 8 bytes in 1 blocks
-==12027==         suppressed: 0 bytes in 0 blocks
-==12027== Rerun with --leak-check=full to see details of leaked memory
-==12027== 
-==12027== For counts of detected and suppressed errors, rerun with: -v
-==12027== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)
diff --git a/src/test/generate_test_data/generateRunSaveTest_EMGLLF.R b/src/test/generate_test_data/generateRunSaveTest_EMGLLF.R
index f6059e5..1c16318 100644
--- a/src/test/generate_test_data/generateRunSaveTest_EMGLLF.R
+++ b/src/test/generate_test_data/generateRunSaveTest_EMGLLF.R
@@ -3,12 +3,12 @@ generateRunSaveTest_EMGLLF = function(n=200, p=15, m=10, k=3, mini=5, maxi=10,
 {
 	testFolder = "data/"
 	dir.create(testFolder, showWarnings=FALSE, mode="0755")
-	delimiter = " "
 
-	#save inputs
 	require(valse)
 	params = valse:::basic_Init_Parameters(n, p, m, k)
 	io = generateIOdefault(n, p, m, k)
+
+	#save inputs
 	write.table(as.double(params$phiInit), paste(testFolder,"phiInit",sep=""),
 		row.names=F, col.names=F)
 	write.table(as.double(params$rhoInit), paste(testFolder,"rhoInit",sep=""),
diff --git a/src/test/generate_test_data/generateRunSaveTest_EMGrank.R b/src/test/generate_test_data/generateRunSaveTest_EMGrank.R
index 0f97995..3740b58 100644
--- a/src/test/generate_test_data/generateRunSaveTest_EMGrank.R
+++ b/src/test/generate_test_data/generateRunSaveTest_EMGrank.R
@@ -1,35 +1,41 @@
-generateRunSaveTest_EMGrank = function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0, rank = c(1,2,4)){
+generateRunSaveTest_EMGrank = function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0,
+	rank = c(1,2,4))
+{
   testFolder = "data/"
   dir.create(testFolder, showWarnings=FALSE, mode="0755")
-  delimiter = " "
-  
+
   tau = 1e-6
-  
   pi = rep(1.0/k, k)
   rho = array(0, dim=c(m,m,k))
-  
   for(i in 1:k){
     rho[,,i] = diag(1,m)
   }
-
-  #Generate X and Y
 	require(valse)
-  generateIOdef = valse:::generateIOdefault(n, p, m, k)
-  
+  io = valse:::generateIOdefault(n, p, m, k)
+
   #save inputs
-  write.table(paste(testFolder,"rho",sep=""), rho, sep=delimiter)
-  write.table(paste(testFolder,"pi",sep=""), pi, sep=delimiter)
-  write.table(paste(testFolder,"mini",sep=""), mini, sep=delimiter)
-  write.table(paste(testFolder,"maxi",sep=""), maxi, sep=delimiter)
-  write.table(paste(testFolder,"X",sep=""), generateIOdef$X sep=delimiter)
-  write.table(paste(testFolder,"Y",sep=""), generateIOdef$Y, sep=delimiter)
-  write.table(paste(testFolder,"tau",sep=""), tau, sep=delimiter)
-  write.table(paste(testFolder,"rank",sep=""), rank, sep=delimiter)
-  write.table(paste(testFolder,"dimensions",sep=""), c(n,p,m,k), sep=delimiter)
-  
-  EMG_rank = EMG(pi,rho,mini,maxi,X,Y,tau,rank)
-  
+  write.table(as.double(rho), paste(testFolder,"rho",sep=""),
+		row.names=F, col.names=F)
+	write.table(as.double(pi), paste(testFolder,"pi",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(mini), paste(testFolder,"mini",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(maxi), paste(testFolder,"maxi",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(tau), paste(testFolder,"tau",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(rank), paste(testFolder,"rank",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensions",sep=""),
+		row.names=F, col.names=F)
+
+  res = EMGrank(pi,rho,mini,maxi,io$X,io$Y,tau,rank)
+
   #save output
-  write.table(paste(testFolder,"phi",sep=""), EMG_rank$phi, sep=delimiter)
-  write.table(paste(testFolder,"LLF",sep=""), EMG_rank$LLF, sep=delimiter)
+  write.table(as.double(res$phi), paste(testFolder,"phi",sep=""), row.names=F,col.names=F)
+  write.table(as.double(res$LLF), paste(testFolder,"LLF",sep=""), row.names=F,col.names=F)
 }
diff --git a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R
index b6ff68c..500ce2f 100644
--- a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R
+++ b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R
@@ -1,46 +1,55 @@
-generateRunSaveTest_constructionModelesLassoMLE = function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0, glambda=list(0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99)){
+generateRunSaveTest_constructionModelesLassoMLE = function(n=200, p=15, m=10, k=3, mini=5,
+	maxi=10, gamma=1.0, glambda=list(0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99))
+{
   testFolder = "data/"
   dir.create(testFolder, showWarnings=FALSE, mode="0755")
-  delimiter = " "
-  
-  #Generate phiInit,piInit,...
+
 	require(valse)
   params = valse:::basic_Init_Parameters(n, p, m, k)
-
-  #save inputs
-  write.table(paste(testFolder,"phiInit",sep=""), params$phiInit, sep=delimiter)
-  write.table(paste(testFolder,"rhoInit",sep=""), params$rhoInit, sep=delimiter)
-  write.table(paste(testFolder,"piInit",sep=""), params$piInit, sep=delimiter)
-  write.table(paste(testFolder,"gamInit",sep=""), params$gamInit, sep=delimiter)
-  write.table(paste(testFolder,"mini",sep=""), mini, sep=delimiter)
-  write.table(paste(testFolder,"maxi",sep=""), maxi, sep=delimiter)
-  write.table(paste(testFolder,"gamma",sep=""), gamma, sep=delimiter)
-  write.table(paste(testFolder,"lambda",sep=""), lambda, sep=delimiter)
-  write.table(paste(testFolder,"X",sep=""), io$X, sep=delimiter)
-  write.table(paste(testFolder,"Y",sep=""), io$Y, sep=delimiter)
-  write.table(paste(testFolder,"tau",sep=""), tau, sep=delimiter)
-  write.table(paste(testFolder,"dimensions",sep=""), c(n,p,m,k), sep=delimiter)
-  
   A2 = array(0, dim=c(p, m+1, L))
   A1 = array(0, dim=c(p, m+1, L))
-  for(i in 1:L){
-    A2[,1,i] = seq(1,p)
-    A2[,2,i] = seq(1,5)
-    A1[,1, i] = seq(1,p)
-    A1[,2,i] = seq(1,5)
+  for (i in 1:L)
+	{
+    A2[,1,i] = seq_len(p)
+    A2[,2,i] = seq_len(5)
+    A1[,1,i] = seq_len(p)
+    A1[,2,i] = seq_len(5)
   }
-  
-  #Generate X and Y
-  generateIOdef = generateIOdefault(n, p, m, k)
-  construct_LME = constructionModelesLassoMLE(params$phiInit,params$rhoInit,params$piInit,params$gamInit,mini,maxi,gamma,glambda,generateIOdef$X,generateIOdef$Y,seuil,tau,A1,A2)
-  phi = construct_LME$phi
-  rho = construct_LME$rho
-  pi = construct_LME$pi
-  lvraisemblance = construct_LME$lvraisemblance
-  
+  io = generateIOdefault(n, p, m, k)
+
+  #save inputs
+  write.table(as.double(params$phiInit), paste(testFolder,"phiInit",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(params$rhoInit), paste(testFolder,"rhoInit",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(params$piInit), paste(testFolder,"piInit",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(params$gamInit), paste(testFolder,"gamInit",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(mini), paste(testFolder,"mini",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(maxi), paste(testFolder,"maxi",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(gamma), paste(testFolder,"gamma",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(lambda), paste(testFolder,"lambda",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(tau), paste(testFolder,"tau",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensions",sep=""),
+		row.names=F, col.names=F)
+
+  res = constructionModelesLassoMLE(
+		params$phiInit,params$rhoInit,params$piInit,params$gamInit,
+		mini,maxi,gamma,glambda,io$X,io$Y,seuil,tau,A1,A2)
+
   #save output
-  write.table(paste(testFolder,"phi",sep=""), construct_LME$phi, sep=delimiter)
-  write.table(paste(testFolder,"rho",sep=""), construct_LME$rho, sep=delimiter)
-  write.table(paste(testFolder,"pi",sep=""), construct_LME$pi, sep=delimiter)
-  write.table(paste(testFolder,"lvraisemblance",sep=""), construct_LME$lvraisemblance, sep=delimiter)
+  write.table(as.double(res$phi), paste(testFolder,"phi",sep=""), row.names=F, col.names=F)
+  write.table(as.double(res$rho), paste(testFolder,"rho",sep=""), row.names=F, col.names=F)
+  write.table(as.double(res$pi), paste(testFolder,"pi",sep=""), row.names=F, col.names=F)
+  write.table(as.double(res$llh), paste(testFolder,"llh",sep=""), row.names=F, col.names=F)
 }
diff --git a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.m b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.m
index 0ba83c1..dde3daf 100644
--- a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.m
+++ b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.m
@@ -51,12 +51,12 @@ function[] = generateRunSaveTest_constructionModelesLassoMLE(n, p, m, k, mini, m
 	dlmwrite(strcat(testFolder,'A2'), reshape(A2,1,[]), delimiter);
 	dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k,L], delimiter);
 
-	[phi,rho,pi,lvraisemblance] = constructionModelesLassoMLE(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2);
+	[phi,rho,pi,llh] = constructionModelesLassoMLE(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2);
 
 	%save output
 	dlmwrite(strcat(testFolder,'phi'), reshape(phi,1,[]), delimiter);
 	dlmwrite(strcat(testFolder,'rho'), reshape(rho,1,[]), delimiter);
 	dlmwrite(strcat(testFolder,'pi'), reshape(pi,1,[]), delimiter);
-	dlmwrite(strcat(testFolder,'lvraisemblance'), reshape(lvraisemblance,1,[]), delimiter);
+	dlmwrite(strcat(testFolder,'llh'), reshape(llh,1,[]), delimiter);
 
 end
diff --git a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.R b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.R
index dfb2736..1a6076c 100644
--- a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.R
+++ b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.R
@@ -1,11 +1,10 @@
-generateRunSaveTest_constructionModelesLassoRank = function(n=200, p=15, m=10, L=12, mini=5, maxi=10, gamma=1.0, rangmin=3, rangmax=6){
+generateRunSaveTest_constructionModelesLassoRank = function(n=200, p=15, m=10, L=12, mini=5,
+	maxi=10, gamma=1.0, rangmin=3, rangmax=6)
+{
   testFolder = "data/"
   dir.create(testFolder, showWarnings=FALSE, mode="0755")
-  delimiter = " "
-  
+
   tau = 1e-6
-  
-  
   pi = matrix(0, k,L)
   for(i in 1:L){
     pi[,i] = rep(1.0/k, k)
@@ -16,31 +15,40 @@ generateRunSaveTest_constructionModelesLassoRank = function(n=200, p=15, m=10, L
       rho[,,r,l] = diag(1,m)
     }
   }
-  #Generate X and Y
 	require(valse)
-  generateIOdef = valse:::generateIOdefault(n, p, m, k)
-  
+  io = valse:::generateIOdefault(n, p, m, k)
   A1 = matrix(0,p,L)
   for(i in 1:L){
     A1[,i] = seq(1,p)
   }
+
   #save inputs
-  write.table(paste(testFolder,"rho",sep=""), rho, sep=delimiter)
-  write.table(paste(testFolder,"pi",sep=""), pi, sep=delimiter)
-  write.table(paste(testFolder,"mini",sep=""), mini, sep=delimiter)
-  write.table(paste(testFolder,"maxi",sep=""), maxi, sep=delimiter)
-  write.table(paste(testFolder,"X",sep=""), generateIOdef$X sep=delimiter)
-  write.table(paste(testFolder,"Y",sep=""), generateIOdef$Y, sep=delimiter)
-  write.table(paste(testFolder,"tau",sep=""), tau, sep=delimiter)
-  write.table(paste(testFolder,"A1",sep=""), A1, sep=delimiter)
-  write.table(paste(testFolder,"rangmin",sep=""), rangmin, sep=delimiter)
-  write.table(paste(testFolder,"rangmax",sep=""), rangmax, sep=delimiter)
-  write.table(paste(testFolder,"dimensions",sep=""), c(n,p,m,k), sep=delimiter)
-  
-  construct = constructionModelesLassoRank(pi,rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax))
-  
+  write.table(as.double(rho),paste(testFolder,"rho",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(pi), paste(testFolder,"pi",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(mini), paste(testFolder,"mini",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(maxi), paste(testFolder,"maxi",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(tau),paste(testFolder,"tau",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(A1),paste(testFolder,"A1",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(rangmin),paste(testFolder,"rangmin",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(rangmax),paste(testFolder,"rangmax",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(c(n,p,m,k)),paste(testFolder,"dimensions",sep=""),
+		row.names=F, col.names=F)
+
+  res = constructionModelesLassoRank(pi,rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax)
+
   #save output
-  write.table(paste(testFolder,"phi",sep=""), construct$phi, sep=delimiter)
-  write.table(paste(testFolder,"lvraisemblance",sep=""), construct$lvraisemblance, sep=delimiter)
-  
+  write.table(as.double(res$phi), paste(testFolder,"phi",sep=""), row.names=F, col.names=F)
+  write.table(as.double(res$llh), paste(testFolder,"llh",sep=""), row.names=F, col.names=F)
 }
diff --git a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.m b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.m
index 80d9db3..a599f19 100644
--- a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.m
+++ b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.m
@@ -49,11 +49,11 @@ function[] = generateRunSaveTest_constructionModelesLassoRank(n, p, m, k, L, min
 	dlmwrite(strcat(testFolder,'rangmax'), rangmax, delimiter);
 	dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k,L], delimiter);
 
-	[phi,lvraisemblance] = constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax);
+	[phi,llh] = constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax);
 
 	%save output
 	Size = (rangmax-rangmin+1)^k;
 	dlmwrite(strcat(testFolder,'phi'), reshape(phi,1,[]), delimiter);
-	dlmwrite(strcat(testFolder,'lvraisemblance'), reshape(lvraisemblance,1,[]), delimiter);
+	dlmwrite(strcat(testFolder,'llh'), reshape(llh,1,[]), delimiter);
 
 end
diff --git a/src/test/generate_test_data/generateRunSaveTest_selectiontotale.R b/src/test/generate_test_data/generateRunSaveTest_selectiontotale.R
index adfe26c..bdac898 100644
--- a/src/test/generate_test_data/generateRunSaveTest_selectiontotale.R
+++ b/src/test/generate_test_data/generateRunSaveTest_selectiontotale.R
@@ -1,36 +1,45 @@
-generateRunSaveTest_selectiontotale= function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0, glambda=list(0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99)){
+generateRunSaveTest_selectiontotale= function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0,
+	glambda=list(0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99))
+{
   testFolder = "data/"
   dir.create(testFolder, showWarnings=FALSE, mode="0755")
-  delimiter = " "
-  
-	require(valse)
 
-  #Generate phiInit,piInit,...
+	require(valse)
   params = valse:::basic_Init_Parameters(n, p, m, k)
-  
-  #Generate X and Y
-  generateIOdef = valse:::generateIOdefault(n, p, m, k)
-  
+  io = valse:::generateIOdefault(n, p, m, k)
+
   #save inputs
-  write.table(paste(testFolder,"phiInit",sep=""), params$phiInit, sep=delimiter)
-  write.table(paste(testFolder,"rhoInit",sep=""), params$rhoInit, sep=delimiter)
-  write.table(paste(testFolder,"piInit",sep=""), params$piInit, sep=delimiter)
-  write.table(paste(testFolder,"gamInit",sep=""), params$gamInit, sep=delimiter)
-  write.table(paste(testFolder,"mini",sep=""), mini, sep=delimiter)
-  write.table(paste(testFolder,"maxi",sep=""), maxi, sep=delimiter)
-  write.table(paste(testFolder,"gamma",sep=""), gamma, sep=delimiter)
-  write.table(paste(testFolder,"lambda",sep=""), glambda, sep=delimiter)
-  write.table(paste(testFolder,"X",sep=""), io$X, sep=delimiter)
-  write.table(paste(testFolder,"Y",sep=""), io$Y, sep=delimiter)
-  write.table(paste(testFolder,"tau",sep=""), tau, sep=delimiter)
-  write.table(paste(testFolder,"dimensions",sep=""), c(n,p,m,k), sep=delimiter)
-  
-  
-  selec = selectiontotale(params$phiInit,params$rhoInit,params$piInit,params$gamInit,mini,maxi,gamma,glambda,generateIOdef$X,generateIOdef$Y,seuil, tau)
-  
+  write.table(as.double(params$phiInit), paste(testFolder,"phiInit",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(params$rhoInit), paste(testFolder,"rhoInit",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(params$piInit), paste(testFolder,"piInit",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(params$gamInit), paste(testFolder,"gamInit",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(mini), paste(testFolder,"mini",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(maxi), paste(testFolder,"maxi",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(gamma), paste(testFolder,"gamma",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(glambda), paste(testFolder,"lambda",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(io$X), paste(testFolder,"X",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(io$Y), paste(testFolder,"Y",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.double(tau), paste(testFolder,"tau",sep=""),
+		row.names=F, col.names=F)
+  write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensions",sep=""),
+		row.names=F, col.names=F)
+
+  res = selectiontotale(params$phiInit,params$rhoInit,params$piInit,params$gamInit,
+		mini,maxi,gamma,glambda,io$X,io$Y,seuil, tau)
+
   #save output
-  write.table(paste(testFolder,"A1",sep=""), selec$A1, sep=delimiter)
-  write.table(paste(testFolder,"A2",sep=""), selec$A2, sep=delimiter)
-  write.table(paste(testFolder,"rho",sep=""), selec$rho, sep=delimiter)
-  write.table(paste(testFolder,"pi",sep=""), selec$pi, sep=delimiter)
+  write.table(as.double(res$A1), paste(testFolder,"A1",sep=""), row.names=F, col.names=F)
+  write.table(as.double(res$A2), paste(testFolder,"A2",sep=""), row.names=F, col.names=F)
+  write.table(as.double(res$rho), paste(testFolder,"rho",sep=""), row.names=F, col.names=F)
+  write.table(as.double(res$pi), paste(testFolder,"pi",sep=""), row.names=F, col.names=F)
 }
diff --git a/src/test/generate_test_data/helpers/EMGrank.R b/src/test/generate_test_data/helpers/EMGrank.R
index a1e3df7..d419813 100644
--- a/src/test/generate_test_data/helpers/EMGrank.R
+++ b/src/test/generate_test_data/helpers/EMGrank.R
@@ -1,3 +1,11 @@
+#helper to always have matrices as arg (TODO: put this elsewhere? improve?)
+matricize <- function(X)
+{
+	if (!is.matrix(X))
+		return (t(as.matrix(X)))
+	return (X)
+}
+
 require(MASS)
 EMGrank = function(Pi, Rho, mini, maxi, X, Y, tau, rank){
   #matrix dimensions
@@ -9,7 +17,7 @@ EMGrank = function(Pi, Rho, mini, maxi, X, Y, tau, rank){
   #init outputs
   phi = array(0, dim=c(p,m,k))
   Z = rep(1, n)
-  Pi = piInit
+#  Pi = piInit
   LLF = 0
   
   #local variables
@@ -20,65 +28,58 @@ EMGrank = function(Pi, Rho, mini, maxi, X, Y, tau, rank){
   
   #main loop
   ite = 1
-  while(ite<=mini || (ite<=maxi && sumDeltaPhi>tau)){
+  while(ite<=mini || (ite<=maxi && sumDeltaPhi>tau))
+	{
     #M step: Mise à jour de Beta (et donc phi)
-    for(r in 1:k){
-      Z_bin = valse:::vec_bin(Z,r)
-      Z_vec = Z_bin$vec #vecteur 0 et 1 aux endroits o? Z==r
-      Z_indice = Z_bin$indice 
-      if(sum(Z_indice) == 0){
+    for(r in 1:k)
+		{
+      Z_indice = seq_len(n)[Z==r] #indices où Z == r
+      if (length(Z_indice) == 0)
         next
-      }
       #U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
-      sv = svd(ginv( crossprod(X[Z_indice,]) )   %*%   crossprod(X[Z_indice,], Y[Z_indice,])     )
-      S = diag(sv$d)
-      U = sv$u
-      V = sv$v
+      s = svd( ginv(crossprod(matricize(X[Z_indice,]))) %*%
+				crossprod(matricize(X[Z_indice,]),matricize(Y[Z_indice,])) )
+      S = s$d
+      U = s$u
+      V = s$v
       #Set m-rank(r) singular values to zero, and recompose
       #best rank(r) approximation of the initial product
-      if(r==k){
-        j_r_1 = length(S)
-      }
-      else{
-        j_r_1 = c(rank[r]+1:length(S))
-      }
-      S[j_r_1] = 0
-      S = diag(S, nrow = ncol(U))
-      phi[,,r] = U %*% S %*% t(V) %*% Rho[,,r]
+      if(rank[r] < length(S))
+        S[(rank[r]+1):length(S)] = 0
+      phi[,,r] = U %*% diag(S) %*% t(V) %*% Rho[,,r]
     }
   
-  #Etape E et calcul de LLF
-  sumLogLLF2 = 0
-  for(i in 1:n){
-    sumLLF1 = 0
-    maxLogGamIR = -Inf
-    for(r in 1:k){
-      dotProduct = tcrossprod(Y[i,]%*%Rho[,,r]-X[i,]%*%phi[,,r])
-      logGamIR = log(Pi[r]) + log(det(Rho[,,r])) - 0.5*dotProduct
-      #Z[i] = index of max (gam[i,])
-      if(logGamIR > maxLogGamIR){
-        Z[i] = r
-        maxLogGamIR = logGamIR
-      }
-    sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2)
-    }
-    sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
-  }
+		#Etape E et calcul de LLF
+		sumLogLLF2 = 0
+		for(i in 1:n){
+			sumLLF1 = 0
+			maxLogGamIR = -Inf
+			for(r in 1:k){
+				dotProduct = tcrossprod(Y[i,]%*%Rho[,,r]-X[i,]%*%phi[,,r])
+				logGamIR = log(Pi[r]) + log(det(Rho[,,r])) - 0.5*dotProduct
+				#Z[i] = index of max (gam[i,])
+				if(logGamIR > maxLogGamIR){
+					Z[i] = r
+					maxLogGamIR = logGamIR
+				}
+			sumLLF1 = sumLLF1 + exp(logGamIR) / (2*pi)^(m/2)
+			}
+			sumLogLLF2 = sumLogLLF2 + log(sumLLF1)
+		}
   
-  LLF = -1/n * sumLogLLF2
-  
-  #update distance parameter to check algorithm convergence (delta(phi, Phi))
-  deltaPhi = c(deltaPhi, max(max(max((abs(phi-Phi))/(1+abs(phi))))) )
-  if(length(deltaPhi) > deltaPhiBufferSize){
-    l_1 = c(2:length(deltaPhi))
-    deltaPhi = deltaPhi[l_1]
-  }
-  sumDeltaPhi = sum(abs(deltaPhi))
+		LLF = -1/n * sumLogLLF2
   
-  #update other local variables
-  Phi = phi
-  ite = ite+1
+		#update distance parameter to check algorithm convergence (delta(phi, Phi))
+		deltaPhi = c(deltaPhi, max(max(max((abs(phi-Phi))/(1+abs(phi))))) )
+		if(length(deltaPhi) > deltaPhiBufferSize){
+		  l_1 = c(2:length(deltaPhi))
+		  deltaPhi = deltaPhi[l_1]
+		}
+		sumDeltaPhi = sum(abs(deltaPhi))
   
+		#update other local variables
+		Phi = phi
+		ite = ite+1
   }
   return(list(phi=phi, LLF=LLF))
 }
diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R b/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R
index 3eac5d1..359bada 100644
--- a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R
+++ b/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R
@@ -10,7 +10,7 @@ constructionModelesLassoMLE = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,
   phi = array(0, dim=c(p,m,k,L))
   rho = array(0, dim=c(m,m,k,L))
   Pi = matrix(0, k, L)
-  lvraisemblance = matrix(0, L, 2)
+  llh = matrix(0, L, 2) #log-likelihood
 
   for(lambdaIndex in 1:L){
     a = A1[, 1, lambdaIndex]
@@ -51,8 +51,8 @@ constructionModelesLassoMLE = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,
 				densite[i,lambdaIndex] = densite[i,lambdaIndex] +	Pi[r,lambdaIndex]*det(rho[,,r,lambdaIndex])/(sqrt(2*pi))^m*exp(-tcrossprod(delta)/2.0)
 			}
 		}
-		lvraisemblance[lambdaIndex,1] = sum(log(densite[,lambdaIndex]))
-		lvraisemblance[lambdaIndex,2] = (dimension+m+1)*k-1
+		llh[lambdaIndex,1] = sum(log(densite[,lambdaIndex]))
+		llh[lambdaIndex,2] = (dimension+m+1)*k-1
   }
-  return(list(phi=phi, rho=rho, Pi=Pi, lvraisemblance = lvraisemblance))
-}
\ No newline at end of file
+  return(list(phi=phi, rho=rho, Pi=Pi, llh = llh))
+}
diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.m b/src/test/generate_test_data/helpers/constructionModelesLassoMLE.m
index 5b7c65e..3e852a6 100644
--- a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.m
+++ b/src/test/generate_test_data/helpers/constructionModelesLassoMLE.m
@@ -1,4 +1,4 @@
-function[phi,rho,pi,lvraisemblance] = constructionModelesLassoMLE(...
+function[phi,rho,pi,llh] = constructionModelesLassoMLE(...
 	phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2)
 
 	PI = 4.0 * atan(1.0);
@@ -12,7 +12,7 @@ function[phi,rho,pi,lvraisemblance] = constructionModelesLassoMLE(...
 	phi = zeros(p,m,k,L);
 	rho = zeros(m,m,k,L);
 	pi = zeros(k,L);
-	lvraisemblance = zeros(L,2);
+	llh = zeros(L,2);
 
 	for lambdaIndex=1:L
 		% Procedure Lasso-MLE
@@ -51,8 +51,8 @@ function[phi,rho,pi,lvraisemblance] = constructionModelesLassoMLE(...
 					pi(r,lambdaIndex)*det(rho(:,:,r,lambdaIndex))/(sqrt(2*PI))^m*exp(-dot(delta,delta)/2.0);
 			end
 		end
-		lvraisemblance(lambdaIndex,1) = sum(log(densite(:,lambdaIndex)));
-		lvraisemblance(lambdaIndex,2) = (dimension+m+1)*k-1;
+		llh(lambdaIndex,1) = sum(log(densite(:,lambdaIndex)));
+		llh(lambdaIndex,2) = (dimension+m+1)*k-1;
 	end
 
 end
diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoRank.R b/src/test/generate_test_data/helpers/constructionModelesLassoRank.R
index 183a6d1..ad4f725 100644
--- a/src/test/generate_test_data/helpers/constructionModelesLassoRank.R
+++ b/src/test/generate_test_data/helpers/constructionModelesLassoRank.R
@@ -14,7 +14,7 @@ constructionModelesLassoRank = function(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rang
 #  }
   
   phi = array(0, dim=c(p,m,k,L*Size))
-  lvraisemblance = matrix(0, L*Size, 2)
+  llh = matrix(0, L*Size, 2) #log-likelihood
   for(lambdaIndex in 1:L){
     #on ne garde que les colonnes actives
     #active sera l'ensemble des variables informatives
@@ -25,10 +25,10 @@ constructionModelesLassoRank = function(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rang
         EMG_rank = EMGrank(Pi[,lambdaIndex], Rho[,,,lambdaIndex], mini, maxi, X[, active], Y, tau, Rank[j,])
         phiLambda = EMG_rank$phi
         LLF = EMG_rank$LLF
-        lvraisemblance[(lambdaIndex-1)*Size+j,] = c(LLF, sum(Rank[j,]^(length(active)- Rank[j,]+m)))
+        llh[(lambdaIndex-1)*Size+j,] = c(LLF, sum(Rank[j,]^(length(active)- Rank[j,]+m)))
         phi[active,,,(lambdaIndex-1)*Size+j] = phiLambda
       }
     }
   }
-  return(list(phi=phi, lvraisemblance = lvraisemblance))
+  return(list(phi=phi, llh = llh))
 }
diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoRank.m b/src/test/generate_test_data/helpers/constructionModelesLassoRank.m
index 415ab12..6279416 100644
--- a/src/test/generate_test_data/helpers/constructionModelesLassoRank.m
+++ b/src/test/generate_test_data/helpers/constructionModelesLassoRank.m
@@ -1,4 +1,4 @@
-function[phi,lvraisemblance] = constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax)
+function[phi,llh] = constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax)
 
 	PI = 4.0 * atan(1.0);
 
@@ -22,7 +22,7 @@ function[phi,lvraisemblance] = constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y
 
 	%output parameters
 	phi = zeros(p,m,k,L*Size);
-	lvraisemblance = zeros(L*Size,2);
+	llh = zeros(L*Size,2);
 	for lambdaIndex=1:L
 		%On ne garde que les colonnes actives
 		%active sera l'ensemble des variables informatives
@@ -31,7 +31,7 @@ function[phi,lvraisemblance] = constructionModelesLassoRank(Pi,Rho,mini,maxi,X,Y
 		if length(active) > 0
 			for j=1:Size
 				[phiLambda,LLF] = EMGrank(Pi(:,lambdaIndex),Rho(:,:,:,lambdaIndex),mini,maxi,X(:,active),Y,tau,Rank(j,:));
-				lvraisemblance((lambdaIndex-1)*Size+j,:) = [LLF, sum(Rank(j,:) .* (length(active)-Rank(j,:)+m))];
+				llh((lambdaIndex-1)*Size+j,:) = [LLF, sum(Rank(j,:) .* (length(active)-Rank(j,:)+m))];
 				phi(active,:,:,(lambdaIndex-1)*Size+j) = phiLambda;
 			end
 		end
diff --git a/src/test/sourceAll.R b/src/test/sourceAll.R
new file mode 100644
index 0000000..6e790a2
--- /dev/null
+++ b/src/test/sourceAll.R
@@ -0,0 +1,6 @@
+for (file in list.files("generate_test_data",full.names=TRUE,recursive=TRUE))
+{
+	file_name_length = nchar(file)
+	if (substr(file, file_name_length, file_name_length) == "R")
+		source(file)
+}
diff --git a/src/test/test.ConstructionModelesLassoMLE.c b/src/test/test.ConstructionModelesLassoMLE.c
index e9c7678..45402db 100644
--- a/src/test/test.ConstructionModelesLassoMLE.c
+++ b/src/test/test.ConstructionModelesLassoMLE.c
@@ -3,108 +3,48 @@
 
 int main(int argc, char** argv)
 {
-	// read dimensions
-	const int nbDims = 5;
-	int* dimensions = readArray_int("dimensions",&nbDims,1);
+	int* dimensions = readArray_int("dimensions");
 	int n = dimensions[0];
 	int p = dimensions[1];
 	int m = dimensions[2];
 	int k = dimensions[3];
 	int L = dimensions[4];
 	free(dimensions);
-	int lengthOne = 1;
 
 	////////////
 	// INPUTS //
+	Real* phiInit = readArray_real("phiInit");
+	Real* rhoInit = readArray_real("rhoInit");
+	Real* piInit = readArray_real("piInit");
+	Real* gamInit = readArray_real("gamInit");
+	int mini = read_int("mini");
+	int maxi = read_int("maxi");
+	Real gamma = read_real("gamma");
+	Real* glambda = readArray_real("glambda");
+	Real* X = readArray_real("X");
+	Real* Y = readArray_real("Y");
+	Real seuil = read_real("seuil");
+	Real tau = read_real("tau");
+	int* A1 = readArray_int("A1");
+	int* A2 = readArray_int("A2");
 	////////////
 
-	// phiInit
-	const int dimPhiInit[] = {p, m, k};
-	float* phiInit = readArray_real("phiInit",dimPhiInit,3);
-
-	// rhoInit
-	const int dimRhoInit[] = {m, m, k};
-	float* rhoInit = readArray_real("rhoInit",dimRhoInit,3);
-
-	// piInit
-	float* piInit = readArray_real("piInit",&k,1);
-
-	// gamInit
-	const int dimGamInit[] = {n, k};
-	float* gamInit = readArray_real("gamInit",dimGamInit,2);
-
-	// min number of iterations
-	int* pmini = readArray_int("mini",&lengthOne,1);
-	int mini = *pmini;
-	free(pmini);
-	
-	// max number of iterations
-	int* pmaxi = readArray_int("maxi",&lengthOne,1);
-	int maxi = *pmaxi;
-	free(pmaxi);
-
-	// gamma
-	float* pgamma = readArray_real("gamma",&lengthOne,1);
-	float gamma = *pgamma;
-	free(pgamma);
-
-	// lambda
-	float* glambda = readArray_real("glambda",&L,1);
-
-	// X
-	const int dimX[] = {n, p};
-	float* X = readArray_real("X",dimX,2);
-
-	// Y
-	const int dimY[] = {n, m};
-	float* Y = readArray_real("Y",dimY,2);
-
-	// seuil
-	float* pseuil = readArray_real("seuil",&lengthOne,1);
-	float seuil = *pseuil;
-	free(pseuil);
-
-	// tau
-	float* ptau = readArray_real("tau",&lengthOne,1);
-	float tau = *ptau;
-	free(ptau);
-	
-	// A1
-	const int dimA[] = {p, m+1, L};
-	int* A1 = readArray_int("A1",dimA,3);
-	
-	// A2
-	int* A2 = readArray_int("A2",dimA,3);
-	
 	/////////////
 	// OUTPUTS //
+	Real* phi = (Real*)malloc(p*m*k*L*sizeof(Real));
+	Real* rho = (Real*)malloc(m*m*k*L*sizeof(Real));
+	Real* pi = (Real*)malloc(k*L*sizeof(Real));
+	Real* llh = (Real*)malloc(L*2*sizeof(Real));
 	/////////////
 
-	// phi
-	const int dimPhi[] = {dimPhiInit[0], dimPhiInit[1], dimPhiInit[2], L};
-	float* phi = (float*)malloc(dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]*sizeof(float));
-
-	// rho
-	const int dimRho[] = {dimRhoInit[0], dimRhoInit[1], dimRhoInit[2], L};
-	float* rho = (float*)malloc(dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]*sizeof(float));
-
-	// pi
-	const int dimPi[] = {k, L};
-	float* pi = (float*)malloc(dimPi[0]*dimPi[1]*sizeof(float));
-
-	// lvraisemblance
-	const int dimLvraisemblance[] = {L, 2};
-	float* lvraisemblance = (float*)malloc(dimLvraisemblance[0]*dimLvraisemblance[1]*sizeof(float));
-
 	/////////////////////////////////////////
 	// Call to constructionModelesLassoMLE //
-	/////////////////////////////////////////
-
 	constructionModelesLassoMLE(
 		phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2,
-		phi,rho,pi,lvraisemblance,
+		phi,rho,pi,llh,
 		n,p,m,k,L);
-	
+	/////////////////////////////////////////
+
 	free(phiInit);
 	free(rhoInit);
 	free(piInit);
@@ -114,30 +54,27 @@ int main(int argc, char** argv)
 	free(A1);
 	free(A2);
 	free(glambda);
-	
+
 	// Compare to reference outputs
-	float* ref_phi = readArray_real("phi",dimPhi,4);
-	compareArray_real("phi", phi, ref_phi, dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]);
+	Real* ref_phi = readArray_real("phi",dimPhi,4);
+	compareArray_real("phi", phi, ref_phi, p*m*k*L);
 	free(phi);
 	free(ref_phi);
-	
-	// rho
-	float* ref_rho = readArray_real("rho",dimRho,4);
-	compareArray_real("rho", rho, ref_rho, dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]);
+
+	Real* ref_rho = readArray_real("rho",dimRho,4);
+	compareArray_real("rho", rho, ref_rho, m*m*k*L);
 	free(rho);
 	free(ref_rho);
-	
-	// pi
-	float* ref_pi = readArray_real("pi",dimPi,2);
-	compareArray_real("pi", pi, ref_pi, dimPi[0]*dimPi[1]);
+
+	Real* ref_pi = readArray_real("pi",dimPi,2);
+	compareArray_real("pi", pi, ref_pi, k*L);
 	free(pi);
 	free(ref_pi);
-	
-	// lvraisemblance
-	float* ref_lvraisemblance = readArray_real("lvraisemblance",dimLvraisemblance,2);
-	compareArray_real("lvraisemblance", lvraisemblance, ref_lvraisemblance, dimLvraisemblance[0]*dimLvraisemblance[1]);
-	free(lvraisemblance);
-	free(ref_lvraisemblance);
+
+	Real* ref_llh = readArray_real("llh",dimllh,2);
+	compareArray_real("llh", llh, ref_llh, L*2);
+	free(llh);
+	free(ref_llh);
 
 	return 0;
 }
diff --git a/src/test/test.EMGLLF.c b/src/test/test.EMGLLF.c
index db38f14..168f81e 100644
--- a/src/test/test.EMGLLF.c
+++ b/src/test/test.EMGLLF.c
@@ -2,8 +2,6 @@
 #include "test_utils.h"
 #include <stdlib.h>
 
-#include <stdio.h>
-
 int main(int argc, char** argv)
 {
 	int* dimensions = readArray_int("dimensions");
@@ -15,36 +13,34 @@ int main(int argc, char** argv)
 
 	////////////
 	// INPUTS //
-	////////////
-
-	float* phiInit = readArray_real("phiInit");
-	float* rhoInit = readArray_real("rhoInit");
-	float* piInit = readArray_real("piInit");
-	float* gamInit = readArray_real("gamInit");
+	Real* phiInit = readArray_real("phiInit");
+	Real* rhoInit = readArray_real("rhoInit");
+	Real* piInit = readArray_real("piInit");
+	Real* gamInit = readArray_real("gamInit");
 	int mini = read_int("mini");
 	int maxi = read_int("maxi");
-	float gamma = read_real("gamma");
-	float lambda = read_real("lambda");
-	float* X = readArray_real("X");
-	float* Y = readArray_real("Y");
-	float tau = read_real("tau");
+	Real gamma = read_real("gamma");
+	Real lambda = read_real("lambda");
+	Real* X = readArray_real("X");
+	Real* Y = readArray_real("Y");
+	Real tau = read_real("tau");
+	////////////
 
 	/////////////
 	// OUTPUTS //
+	Real* phi = (Real*)malloc(p*m*k*sizeof(Real));
+	Real* rho = (Real*)malloc(m*m*k*sizeof(Real));
+	Real* pi = (Real*)malloc(k*sizeof(Real));
+	Real* LLF = (Real*)malloc(maxi*sizeof(Real));
+	Real* S = (Real*)malloc(p*m*k*sizeof(Real));
 	/////////////
 
-	float* phi = (float*)malloc(p*m*k*sizeof(float));
-	float* rho = (float*)malloc(m*m*k*sizeof(float));
-	float* pi = (float*)malloc(k*sizeof(float));
-	float* LLF = (float*)malloc(maxi*sizeof(float));
-	float* S = (float*)malloc(p*m*k*sizeof(float));
-
 	////////////////////
 	// Call to EMGLLF //
-	////////////////////
 	EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau,
 		phi,rho,pi,LLF,S,
 		n,p,m,k);
+	////////////////////
 
 	free(phiInit);
 	free(rhoInit);
@@ -54,27 +50,27 @@ int main(int argc, char** argv)
 	free(Y);
 
 	// Compare to reference outputs
-	float* ref_phi = readArray_real("phi");
+	Real* ref_phi = readArray_real("phi");
 	compareArray_real("phi", phi, ref_phi, p*m*k);
 	free(phi);
 	free(ref_phi);
 
-	float* ref_rho = readArray_real("rho");
+	Real* ref_rho = readArray_real("rho");
 	compareArray_real("rho", rho, ref_rho, m*m*k);
 	free(rho);
 	free(ref_rho);
 
-	float* ref_pi = readArray_real("pi");
+	Real* ref_pi = readArray_real("pi");
 	compareArray_real("pi", pi, ref_pi, k);
 	free(pi);
 	free(ref_pi);
 
-	float* ref_LLF = readArray_real("LLF");
+	Real* ref_LLF = readArray_real("LLF");
 	compareArray_real("LLF", LLF, ref_LLF, maxi);
 	free(LLF);
 	free(ref_LLF);
 
-	float* ref_S = readArray_real("S");
+	Real* ref_S = readArray_real("S");
 	compareArray_real("S", S, ref_S, p*m*k);
 	free(S);
 	free(ref_S);
diff --git a/src/test/test.EMGrank.c b/src/test/test.EMGrank.c
index 80263a0..8374b2f 100644
--- a/src/test/test.EMGrank.c
+++ b/src/test/test.EMGrank.c
@@ -1,5 +1,6 @@
 #include "EMGrank.h"
 #include "test_utils.h"
+#include <stdlib.h>
 
 int main(int argc, char** argv)
 {
@@ -12,46 +13,43 @@ int main(int argc, char** argv)
 
 	////////////
 	// INPUTS //
-	////////////
-
-	float* Rho = readArray_real("Rho");
-	float* Pi = readArray_real("Pi");
+	Real* rho = readArray_real("rho");
+	Real* pi = readArray_real("pi");
 	int mini = read_int("mini");
 	int maxi = read_int("maxi");
-	float* X = readArray_real("X");
-	float* Y = readArray_real("Y");
-	float tau = read_real("tau");
+	Real* X = readArray_real("X");
+	Real* Y = readArray_real("Y");
+	Real tau = read_real("tau");
 	int* rank = readArray_int("rank");
+	////////////
 
 	/////////////
 	// OUTPUTS //
+	Real* phi = (Real*)malloc(p*m*k*sizeof(Real));
+	Real* LLF = (Real*)malloc(1*sizeof(Real));
 	/////////////
 
-	float* phi = (float*)malloc(p*m*k*sizeof(float));
-	float* LLF = (float*)malloc(1*sizeof(float));
-
 	//////////////////////////
 	// Main call to EMGrank //
-	//////////////////////////
-
-	EMGrank(Pi,Rho,mini,maxi,X,Y,tau,rank,
+	EMGrank_core(pi,rho,mini,maxi,X,Y,tau,rank,
 		phi,LLF,
 		n,p,m,k);
+	//////////////////////////
 
-	free(Rho);
-	free(Pi);
+	free(rho);
+	free(pi);
 	free(X);
 	free(Y);
 	free(rank);
 
 	// Compare to reference outputs
-	float* ref_phi = readArray_real("phi");
+	Real* ref_phi = readArray_real("phi");
 	compareArray_real("phi", phi, ref_phi, p*m*k);
 	free(phi);
 	free(ref_phi);
 
 	// LLF
-	float* ref_LLF = readArray_real("LLF");
+	Real* ref_LLF = readArray_real("LLF");
 	compareArray_real("LLF", LLF, ref_LLF, 1);
 	free(LLF);
 	free(ref_LLF);
diff --git a/src/test/test.constructionModelesLassoRank.c b/src/test/test.constructionModelesLassoRank.c
index cf95bc4..b3ed34a 100644
--- a/src/test/test.constructionModelesLassoRank.c
+++ b/src/test/test.constructionModelesLassoRank.c
@@ -1,107 +1,62 @@
 #include "constructionModelesLassoRank.h"
 #include "test_utils.h"
+#include <stdlib.h>
 
 int main(int argc, char** argv)
 {
-	// read dimensions
-	const int nbDims = 5;
-	int* dimensions = readArray_int("dimensions",&nbDims,1);
+	int* dimensions = readArray_int("dimensions");
 	int n = dimensions[0];
 	int p = dimensions[1];
 	int m = dimensions[2];
 	int k = dimensions[3];
 	int L = dimensions[4];
 	free(dimensions);
-	int lengthOne = 1;
 
 	////////////
 	// INPUTS //
+	Real* Pi = readArray_real("Pi");
+	Real* Rho = readArray_real("Rho");
+	int mini = read_int("mini");
+	int maxi = read_int("maxi");
+	Real* X = readArray_real("X");
+	Real* Y = readArray_real("Y");
+	Real tau = read_real("tau");
+	int* A1 = readArray_int("A1");
+	int rangmin = read_int("rangmin");
+	int rangmax = read_int("rangmax");
 	////////////
 
-	// piInit
-	const int dimPi[] = {k, L};
-	float* Pi = readArray_real("Pi",dimPi,2);
-
-	// rhoInit
-	const int dimRho[] = {m, m, k, L};
-	float* Rho = readArray_real("Rho",dimRho,4);
-
-	// min number of iterations
-	int* pmini = readArray_int("mini",&lengthOne,1);
-	int mini = *pmini;
-	free(pmini);
-
-	// max number of iterations
-	int* pmaxi = readArray_int("maxi",&lengthOne,1);
-	int maxi = *pmaxi;
-	free(pmaxi);
-
-	// X
-	const int dimX[] = {n, p};
-	float* X = readArray_real("X",dimX,2);
-
-	// Y
-	const int dimY[] = {n, m};
-	float* Y = readArray_real("Y",dimY,2);
-
-	// tau
-	float* ptau = readArray_real("tau",&lengthOne,1);
-	float tau = *ptau;
-	free(ptau);
-
-	// A1
-	const int dimA[] = {p, L};
-	int* A1 = readArray_int("A1",dimA,2);
-
-	// rangmin
-	int* prangmin = readArray_int("rangmin",&lengthOne,1);
-	int rangmin = *prangmin;
-	free(prangmin);
-
-	// rangmax
-	int* prangmax = readArray_int("rangmax",&lengthOne,1);
-	int rangmax = *prangmax;
-	free(prangmax);
-	
 	/////////////
 	// OUTPUTS //
-	/////////////
-
-	// phi
 	int Size = (int)pow(rangmax-rangmin+1, k);
-	const int dimPhi[] = {p, m, k, L*Size};
-	float* phi = (float*)malloc(dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]*sizeof(float));
-
-	// lvraisemblance
-	const int dimLvraisemblance[] = {L*Size, 2};
-	float* lvraisemblance = (float*)malloc(dimLvraisemblance[0]*dimLvraisemblance[1]*sizeof(float));
-
-	//////////////////////////////////////////////
-	// Main call to constructionModelesLassoMLE //
-	//////////////////////////////////////////////
+	Real* phi = (Real*)malloc(p*m*k*L*Size*sizeof(Real));
+	Real* llh = (Real*)malloc(L*Size*2*sizeof(Real));
+	/////////////
 
-	constructionModelesLassoRank(
+	/////////////////////////////////////////
+	// Call to constructionModelesLassoMLE //
+	constructionModelesLassoRank_core(
 		Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax,
-		phi,lvraisemblance,
+		phi,llh,
 		n,p,m,k,L);
-	
+	/////////////////////////////////////////
+
 	free(Rho);
 	free(Pi);
 	free(X);
 	free(Y);
 	free(A1);
-	
+
 	// Compare to reference outputs
-	float* ref_phi = readArray_real("phi",dimPhi, 4);
-	compareArray_real("phi", phi, ref_phi, dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]);
+	Real* ref_phi = readArray_real("phi");
+	compareArray_real("phi", phi, ref_phi, p*m*k*L*Size);
 	free(phi);
 	free(ref_phi);
-	
-	// lvraisemblance
-	float* ref_lvraisemblance = readArray_real("lvraisemblance",dimLvraisemblance,2);
-	compareArray_real("lvraisemblance", lvraisemblance, ref_lvraisemblance, dimLvraisemblance[0]*dimLvraisemblance[1]);
-	free(lvraisemblance);
-	free(ref_lvraisemblance);
-	
+
+	Real* ref_llh = readArray_real("llh");
+	compareArray_real("llh", llh, ref_llh, L*Size*2);
+	free(llh);
+	free(ref_llh);
+
 	return 0;
 }
diff --git a/src/test/test.selectiontotale.c b/src/test/test.selectiontotale.c
index 2553d05..843a04d 100644
--- a/src/test/test.selectiontotale.c
+++ b/src/test/test.selectiontotale.c
@@ -1,103 +1,49 @@
 #include "selectiontotale.h"
 #include "test_utils.h"
+#include <stdlib.h>
 
 int main(int argc, char** argv)
 {
-	// read dimensions
-	const int nbDims = 5;
-	int* dimensions = readArray_int("dimensions",&nbDims,1);
+	int* dimensions = readArray_int("dimensions");
 	int n = dimensions[0];
 	int p = dimensions[1];
 	int m = dimensions[2];
 	int k = dimensions[3];
 	int L = dimensions[4];
 	free(dimensions);
-	int lengthOne = 1;
-	
+
 	////////////
 	// INPUTS //
+	Real* phiInit = readArray_real("phiInit");
+	Real* rhoInit = readArray_real("rhoInit");
+	Real* piInit = readArray_real("piInit");
+	Real* gamInit = readArray_real("gamInit");
+	int mini = read_int("mini");
+	int maxi = read_int("maxi");
+	Real gamma = read_real("gamma");
+	Real* glambda = readArray_real("glambda");
+	Real* X = readArray_real("X");
+	Real* Y = readArray_real("Y");
+	Real seuil = read_real("seuil");
+	Real tau = read_real("tau");
 	////////////
 
-	// phiInit
-	const int dimPhiInit[] = {p, m, k};
-	float* phiInit = readArray_real("phiInit",dimPhiInit,3);
-
-	// rhoInit
-	const int dimRhoInit[] = {m, m, k};
-	float* rhoInit = readArray_real("rhoInit",dimRhoInit,3);
-
-	// piInit
-	float* piInit = readArray_real("piInit",&k,1);
-
-	// gamInit
-	const int dimGamInit[] = {n, k};
-	float* gamInit = readArray_real("gamInit",dimGamInit, 2);
-
-	// min number of iterations
-	int* pmini = readArray_int("mini",&lengthOne,1);
-	int mini = *pmini;
-	free(pmini);
-	
-	// max number of iterations
-	int* pmaxi = readArray_int("maxi",&lengthOne,1);
-	int maxi = *pmaxi;
-	free(pmaxi);
-	
-	// gamma
-	float* pgamma = readArray_real("gamma",&lengthOne,1);
-	float gamma = *pgamma;
-	free(pgamma);
-	
-	// lambda
-	float* glambda = readArray_real("glambda",&L,1);
-
-	// X
-	const int dimX[] = {n, p};
-	float* X = readArray_real("X",dimX,2);
-
-	// Y
-	const int dimY[] = {n, m};
-	float* Y = readArray_real("Y",dimY,2);
-
-	// seuil
-	float* pseuil = readArray_real("seuil",&lengthOne,1);
-	float seuil = *pseuil;
-	free(pseuil);
-	
-	// tau
-	float* ptau = readArray_real("tau",&lengthOne,1);
-	float tau = *ptau;
-	free(ptau);
-	
 	/////////////
 	// OUTPUTS //
+	int* A1 = (int*)malloc(p*(m+1)*L*sizeof(int));
+	int* A2 = (int*)malloc(p*(m+1)*L*sizeof(int));
+	Real* Rho = (Real*)malloc(m*m*k*L*sizeof(Real));
+	Real* Pi = (Real*)malloc(k*L*sizeof(Real));
 	/////////////
 
-	// A1
-	const int dimA[] = {p, m+1, L};
-	int* A1 = (int*)malloc(dimA[0]*dimA[1]*dimA[2]*sizeof(int));
-
-	// A2
-	int* A2 = (int*)malloc(dimA[0]*dimA[1]*dimA[2]*sizeof(int));
-
-	// rho
-	const int dimRho[] = {m, m, k, L};
-	float* Rho = (float*)malloc(dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]*sizeof(float));
-
-	// pi
-	const int dimPi[] = {k, L};
-	float* Pi = (float*)malloc(dimPi[0]*dimPi[1]*sizeof(float));
-
-	//////////////////////////////////////////////
-	// Main call to constructionModelesLassoMLE //
-	//////////////////////////////////////////////
-
-	selectiontotale(
+	/////////////////////////////////////////
+	// Call to constructionModelesLassoMLE //
+	selectiontotale_core(
 		phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,
 		A1,A2,Rho,Pi,
 		n,p,m,k,L);
-	
-	// free input pointers
+	/////////////////////////////////////////
+
 	free(phiInit);
 	free(rhoInit);
 	free(piInit);
@@ -105,30 +51,27 @@ int main(int argc, char** argv)
 	free(glambda);
 	free(X);
 	free(Y);
-	
+
 	// Compare to reference outputs
-	int* ref_A1 = readArray_int("A1",dimA, 3);
-	compareArray_int("A1", A1, ref_A1, dimA[0]*dimA[1]*dimA[2]);
+	int* ref_A1 = readArray_int("A1");
+	compareArray_int("A1", A1, ref_A1, p*(m+1)*L);
 	free(A1);
 	free(ref_A1);
-	
-	// A2
-	int* ref_A2 = readArray_int("A2",dimA, 3);
-	compareArray_int("A2", A2, ref_A2, dimA[0]*dimA[1]*dimA[2]);
+
+	int* ref_A2 = readArray_int("A2");
+	compareArray_int("A2", A2, ref_A2, p*(m+1)*L);
 	free(A2);
 	free(ref_A2);
-	
-	// Rho
-	float* ref_Rho = readArray_real("Rho",dimRho, 4);
-	compareArray_real("Rho", Rho, ref_Rho, dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]);
+
+	Real* ref_Rho = readArray_real("Rho");
+	compareArray_real("Rho", Rho, ref_Rho, m*m*k*L);
 	free(Rho);
 	free(ref_Rho);
-	
-	// Pi
-	float* ref_Pi = readArray_real("Pi",dimPi, 2);
-	compareArray_real("Pi", Pi, ref_Pi, dimPi[0]*dimPi[1]);
+
+	Real* ref_Pi = readArray_real("Pi");
+	compareArray_real("Pi", Pi, ref_Pi, k*L);
 	free(Pi);
 	free(ref_Pi);
-	
+
 	return 0;
 }
diff --git a/src/test/test_utils.c b/src/test/test_utils.c
index 7fc240e..96972ee 100644
--- a/src/test/test_utils.c
+++ b/src/test/test_utils.c
@@ -8,14 +8,14 @@
 void compareArray(const char* ID, const void* array, const void* refArray, int size,
 	int isinteger)
 {
-	float EPS = 1e-5; //precision
+	Real EPS = 1e-5; //precision
 	printf("Checking %s\n",ID);
-	float maxError = 0.0;
+	Real maxError = 0.0;
 	for (int i=0; i<size; i++)
 	{
-		float error = isinteger
+		Real error = isinteger
 			? fabs(((int*)array)[i] - ((int*)refArray)[i])
-			: fabs(((float*)array)[i] - ((float*)refArray)[i]);
+			: fabs(((Real*)array)[i] - ((Real*)refArray)[i]);
 		if (error >= maxError)
 			maxError = error;
 	}
@@ -48,6 +48,7 @@ void* readArray(const char* fileName, int isinteger)
 	strcat(command, "wc -l ");
 	strcat(command, fullFileName);
 	FILE *arraySize = popen(command, "r");
+	free(command);
 	char* bufferNum = (char*)calloc(64, sizeof(char));
 	fgets(bufferNum, sizeof(bufferNum), arraySize);
 	int n = atoi(bufferNum);
@@ -58,16 +59,16 @@ void* readArray(const char* fileName, int isinteger)
 	free(fullFileName);
 
 	// read all values, and convert them to by-rows matrices format
-	size_t elementSize = isinteger ? sizeof(int) : sizeof(float);
+	size_t elementSize = isinteger ? sizeof(int) : sizeof(Real);
 	void* array = malloc(n*elementSize);
 	for (int i=0; i<n; i++)
 	{
 		fgets(bufferNum, 64, arrayFile);
-		// transform buffer content into float or int, and store it at appropriate location
+		// transform buffer content into Real or int, and store it at appropriate location
 		if (isinteger)
 			((int*)array)[i] = atoi(bufferNum);
 		else
-			((float*)array)[i] = atof(bufferNum);
+			((Real*)array)[i] = atof(bufferNum);
 	}
 	fclose(arrayFile);
 	free(bufferNum);
@@ -80,9 +81,9 @@ int* readArray_int(const char* fileName)
 	return (int*)readArray(fileName, 1);
 }
 
-float* readArray_real(const char* fileName)
+Real* readArray_real(const char* fileName)
 {
-	return (float*)readArray(fileName, 0);
+	return (Real*)readArray(fileName, 0);
 }
 
 int read_int(const char* fileName)
@@ -93,7 +94,7 @@ int read_int(const char* fileName)
 	return res;
 }
 
-float read_real(const char* fileName)
+Real read_real(const char* fileName)
 {
 	Real* array = readArray_real(fileName);
 	Real res = array[0];
diff --git a/src/test/test_utils.h b/src/test/test_utils.h
index cbffc24..e0aec96 100644
--- a/src/test/test_utils.h
+++ b/src/test/test_utils.h
@@ -10,8 +10,8 @@ void* readArray(const char* fileName, int isInteger);
 
 int* readArray_int(const char* fileName);
 
-float* readArray_real(const char* fileName);
+Real* readArray_real(const char* fileName);
 
 int read_int(const char* fileName);
 
-float read_real(const char* fileName);
+Real read_real(const char* fileName);
-- 
2.44.0