From: Benjamin Auder <benjamin.auder@somewhere>
Date: Sun, 2 Apr 2017 09:32:20 +0000 (+0200)
Subject: add affect[ations] in EMGLLF.c return
X-Git-Url: https://git.auder.net/variants/img/pieces/doc/doc/html/left.jpg?a=commitdiff_plain;h=8be79c465ecbbe849c9ee43e2a25c2760134e07a;p=valse.git

add affect[ations] in EMGLLF.c return
---

diff --git a/pkg/src/adapters/a.EMGLLF.c b/pkg/src/adapters/a.EMGLLF.c
index 4ef4396..f24cd2a 100644
--- a/pkg/src/adapters/a.EMGLLF.c
+++ b/pkg/src/adapters/a.EMGLLF.c
@@ -46,7 +46,7 @@ SEXP EMGLLF(
 	// OUTPUTS //
 	/////////////
 
-	SEXP phi, rho, pi, LLF, S, dimPhiS, dimRho;
+	SEXP phi, rho, pi, LLF, S, affec, dimPhiS, dimRho;
 	PROTECT(dimPhiS = allocVector(INTSXP, 3));
 	int* pDimPhiS = INTEGER(dimPhiS);
 	pDimPhiS[0] = p; pDimPhiS[1] = m; pDimPhiS[2] = k;
@@ -58,22 +58,25 @@ SEXP EMGLLF(
 	PROTECT(pi = allocVector(REALSXP, k));
 	PROTECT(LLF = allocVector(REALSXP, maxi-mini+1));
 	PROTECT(S = allocArray(REALSXP, dimPhiS));
+	PROTECT(affec = allocVector(INTSXP, n));
 	double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pLLF=REAL(LLF), *pS=REAL(S);
+	int *pAffec=INTEGER(affec);
 
 	////////////////////
 	// Call to EMGLLF //
 	////////////////////
 
 	EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau,
-		pPhi,pRho,pPi,pLLF,pS,
+		pPhi,pRho,pPi,pLLF,pS,pAffec,
 		n,p,m,k);
 
 	// Build list from OUT params and return it
 	SEXP listParams, listNames;
-	PROTECT(listParams = allocVector(VECSXP, 5));
-	char* lnames[5] = {"phi", "rho", "pi", "LLF", "S"}; //lists labels
-	PROTECT(listNames = allocVector(STRSXP,5));
-	for (int i=0; i<5; i++)
+	int nouts = 6;
+	PROTECT(listParams = allocVector(VECSXP, nouts));
+	char* lnames[6] = {"phi", "rho", "pi", "LLF", "S", "affec"}; //lists labels
+	PROTECT(listNames = allocVector(STRSXP,nouts));
+	for (int i=0; i<nouts; i++)
 		SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
 	setAttrib(listParams, R_NamesSymbol, listNames);
 	SET_VECTOR_ELT(listParams, 0, phi);
@@ -81,7 +84,8 @@ SEXP EMGLLF(
 	SET_VECTOR_ELT(listParams, 2, pi);
 	SET_VECTOR_ELT(listParams, 3, LLF);
 	SET_VECTOR_ELT(listParams, 4, S);
+	SET_VECTOR_ELT(listParams, 5, affec);
 
-	UNPROTECT(9);
+	UNPROTECT(10);
 	return listParams;
 }
diff --git a/pkg/src/sources/EMGLLF.c b/pkg/src/sources/EMGLLF.c
index 520d57b..e019588 100644
--- a/pkg/src/sources/EMGLLF.c
+++ b/pkg/src/sources/EMGLLF.c
@@ -22,6 +22,7 @@ void EMGLLF_core(
 	Real* pi, // parametre des proportions renormalisé, calculé par l'EM
 	Real* LLF, // log vraisemblance associée à cet échantillon, pour les valeurs estimées des paramètres
 	Real* S,
+	int* affec,
 	// additional size parameters
 	int n, // nombre d'echantillons
 	int p, // nombre de covariables
@@ -381,7 +382,20 @@ void EMGLLF_core(
 		ite++;
 	}
 
-	//TODO: affec = apply(gam, 1,which.max) à traduire, et retourner affec aussi.
+	//affec = apply(gam, 1, which.max)
+	for (int i=0; i<n; i++)
+	{
+		Real rowMax = 0.;
+		affec[i] = 0;
+		for (int j=0; j<k; j++)
+		{
+			if (gam[mi(i,j,n,k)] > rowMax)
+			{
+				affec[i] = j+1; //R indices start at 1
+				rowMax = gam[mi(i,j,n,k)];
+			}
+		}
+	}
 
 	//free memory
 	free(b);
diff --git a/pkg/src/sources/EMGLLF.h b/pkg/src/sources/EMGLLF.h
index 8f375ff..e15cb87 100644
--- a/pkg/src/sources/EMGLLF.h
+++ b/pkg/src/sources/EMGLLF.h
@@ -22,6 +22,7 @@ void EMGLLF_core(
 	Real* pi,
 	Real* LLF,
 	Real* S,
+	int* affec,
 	// additional size parameters
 	int n,
 	int p,
diff --git a/test/generate_test_data/EMGLLF.R b/test/generate_test_data/EMGLLF.R
index fce0a8f..039e291 100644
--- a/test/generate_test_data/EMGLLF.R
+++ b/test/generate_test_data/EMGLLF.R
@@ -151,6 +151,6 @@ EMGLLF_R = function(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,ta
     ite = ite+1
   }
   
-  affec = apply(gam, 1,which.max)
+  affec = apply(gam, 1, which.max)
   return(list("phi"=phi, "rho"=rho, "pi"=pi, "LLF"=LLF, "S"=S, "affec" = affec ))
 }
diff --git a/test/generate_test_data/generateRunSaveTest_EMGLLF.R b/test/generate_test_data/generateRunSaveTest_EMGLLF.R
index 77c0e71..e0bb3b8 100644
--- a/test/generate_test_data/generateRunSaveTest_EMGLLF.R
+++ b/test/generate_test_data/generateRunSaveTest_EMGLLF.R
@@ -45,4 +45,5 @@ generateRunSaveTest_EMGLLF = function(n=200, p=15, m=10, k=3, mini=5, maxi=10,
 	write.table(as.double(res$pi), paste(testFolder,"pi",sep=""), row.names=F, col.names=F)
 	write.table(as.double(res$LLF), paste(testFolder,"LLF",sep=""), row.names=F, col.names=F)
 	write.table(as.double(res$S), paste(testFolder,"S",sep=""), row.names=F, col.names=F)
+	write.table(as.integer(res$affec), paste(testFolder,"affec",sep=""), row.names=F, col.names=F)
 }
diff --git a/test/test.EMGLLF.c b/test/test.EMGLLF.c
index 168f81e..7eed301 100644
--- a/test/test.EMGLLF.c
+++ b/test/test.EMGLLF.c
@@ -33,12 +33,13 @@ int main(int argc, char** argv)
 	Real* pi = (Real*)malloc(k*sizeof(Real));
 	Real* LLF = (Real*)malloc(maxi*sizeof(Real));
 	Real* S = (Real*)malloc(p*m*k*sizeof(Real));
+	int* affec = (int*)malloc(n*sizeof(int));
 	/////////////
 
 	////////////////////
 	// Call to EMGLLF //
 	EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau,
-		phi,rho,pi,LLF,S,
+		phi,rho,pi,LLF,S,affec,
 		n,p,m,k);
 	////////////////////
 
@@ -75,5 +76,10 @@ int main(int argc, char** argv)
 	free(S);
 	free(ref_S);
 
+	int* ref_affec = readArray_int("affec");
+	compareArray_int("affec", affec, ref_affec, n);
+	free(affec);
+	free(ref_affec);
+
 	return 0;
 }