add affect[ations] in EMGLLF.c return
authorBenjamin Auder <benjamin.auder@somewhere>
Sun, 2 Apr 2017 09:32:20 +0000 (11:32 +0200)
committerBenjamin Auder <benjamin.auder@somewhere>
Sun, 2 Apr 2017 09:32:20 +0000 (11:32 +0200)
pkg/src/adapters/a.EMGLLF.c
pkg/src/sources/EMGLLF.c
pkg/src/sources/EMGLLF.h
test/generate_test_data/EMGLLF.R
test/generate_test_data/generateRunSaveTest_EMGLLF.R
test/test.EMGLLF.c

index 4ef4396..f24cd2a 100644 (file)
@@ -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;
 }
index 520d57b..e019588 100644 (file)
@@ -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);
index 8f375ff..e15cb87 100644 (file)
@@ -22,6 +22,7 @@ void EMGLLF_core(
        Real* pi,
        Real* LLF,
        Real* S,
+       int* affec,
        // additional size parameters
        int n,
        int p,
index fce0a8f..039e291 100644 (file)
@@ -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 ))
 }
index 77c0e71..e0bb3b8 100644 (file)
@@ -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)
 }
index 168f81e..7eed301 100644 (file)
@@ -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;
 }