progress in debug: fix LLF/llh mismatch, and length + add adapter test in test/

diff --git a/pkg/R/EMGLLF.R b/pkg/R/EMGLLF.R
index 08ff203..57638f9 100644 (file)
--- a/pkg/R/EMGLLF.R
+++ b/pkg/R/EMGLLF.R
@@ -40,7 +40,7 @@ EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
   k <- length(piInit)  #nombre de composantes dans le mélange
   .Call("EMGLLF", phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, 
     X, Y, eps, phi = double(p * m * k), rho = double(m * m * k), pi = double(k), 
   k <- length(piInit)  #nombre de composantes dans le mélange
   .Call("EMGLLF", phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda, 
     X, Y, eps, phi = double(p * m * k), rho = double(m * m * k), pi = double(k), 

-    LLF = double(maxi), S = double(p * m * k), affec = integer(n), n, p, m, k, 
+    llh = double(1), S = double(p * m * k), affec = integer(n), n, p, m, k, 

     PACKAGE = "valse")
 }
 
     PACKAGE = "valse")
 }
 

diff --git a/pkg/R/EMGrank.R b/pkg/R/EMGrank.R
index b85a0fa..f2bf58e 100644 (file)
--- a/pkg/R/EMGrank.R
+++ b/pkg/R/EMGrank.R
@@ -8,7 +8,7 @@
 #' @param maxi Nombre maximal d'itérations dans l'algorithme EM
 #' @param X Régresseurs
 #' @param Y Réponse
 #' @param maxi Nombre maximal d'itérations dans l'algorithme EM
 #' @param X Régresseurs
 #' @param Y Réponse

-#' @param tau Seuil pour accepter la convergence
+#' @param eps Seuil pour accepter la convergence

 #' @param rank Vecteur des rangs possibles
 #'
 #' @return A list ...
 #' @param rank Vecteur des rangs possibles
 #'
 #' @return A list ...


@@ -16,12 +16,12 @@
 #'   LLF : log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres
 #'
 #' @export
 #'   LLF : log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres
 #'
 #' @export

-EMGrank <- function(Pi, Rho, mini, maxi, X, Y, tau, rank, fast = TRUE)
+EMGrank <- function(Pi, Rho, mini, maxi, X, Y, eps, rank, fast = TRUE)

 {
   if (!fast)
   {
     # Function in R
 {
   if (!fast)
   {
     # Function in R

-    return(.EMGrank_R(Pi, Rho, mini, maxi, X, Y, tau, rank))
+    return(.EMGrank_R(Pi, Rho, mini, maxi, X, Y, eps, rank))

   }
 
   # Function in C
   }
 
   # Function in C


@@ -29,7 +29,7 @@ EMGrank <- function(Pi, Rho, mini, maxi, X, Y, tau, rank, fast = TRUE)
   p <- ncol(X)  #nombre de covariables
   m <- ncol(Y)  #taille de Y (multivarié)
   k <- length(Pi)  #nombre de composantes dans le mélange
   p <- ncol(X)  #nombre de covariables
   m <- ncol(Y)  #taille de Y (multivarié)
   k <- length(Pi)  #nombre de composantes dans le mélange

-  .Call("EMGrank", Pi, Rho, mini, maxi, X, Y, tau, rank, phi = double(p * m * k), 
+  .Call("EMGrank", Pi, Rho, mini, maxi, X, Y, eps, rank, phi = double(p * m * k), 

     LLF = double(1), n, p, m, k, PACKAGE = "valse")
 }
 
     LLF = double(1), n, p, m, k, PACKAGE = "valse")
 }
 


@@ -43,7 +43,7 @@ matricize <- function(X)
 }
 
 # R version - slow but easy to read
 }
 
 # R version - slow but easy to read

-.EMGrank_R <- function(Pi, Rho, mini, maxi, X, Y, tau, rank)
+.EMGrank_R <- function(Pi, Rho, mini, maxi, X, Y, eps, rank)

 {
   # matrix dimensions
   n <- nrow(X)
 {
   # matrix dimensions
   n <- nrow(X)


@@ -64,7 +64,7 @@ matricize <- function(X)
   
   # main loop
   ite <- 1
   
   # main loop
   ite <- 1

-  while (ite <= mini || (ite <= maxi && sumDeltaPhi > tau))
+  while (ite <= mini || (ite <= maxi && sumDeltaPhi > eps))

   {
     # M step: update for Beta ( and then phi)
     for (r in 1:k)
   {
     # M step: update for Beta ( and then phi)
     for (r in 1:k)

diff --git a/pkg/R/computeGridLambda.R b/pkg/R/computeGridLambda.R
index 8449d10..ac0788a 100644 (file)
--- a/pkg/R/computeGridLambda.R
+++ b/pkg/R/computeGridLambda.R
@@ -11,13 +11,13 @@
 #' @param gamma power of weights in the penalty
 #' @param mini minimum number of iterations in EM algorithm
 #' @param maxi maximum number of iterations in EM algorithm
 #' @param gamma power of weights in the penalty
 #' @param mini minimum number of iterations in EM algorithm
 #' @param maxi maximum number of iterations in EM algorithm

-#' @param tau threshold to stop EM algorithm
+#' @param eps threshold to stop EM algorithm

 #'
 #' @return the grid of regularization parameters
 #'
 #' @export
 computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini, 
 #'
 #' @return the grid of regularization parameters
 #'
 #' @export
 computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini, 

-  maxi, tau, fast)
+  maxi, eps, fast)

 {
   n <- nrow(X)
   p <- ncol(X)
 {
   n <- nrow(X)
   p <- ncol(X)


@@ -25,7 +25,8 @@ computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mi
   k <- length(piInit)
 
   list_EMG <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda = 0, 
   k <- length(piInit)
 
   list_EMG <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda = 0, 

-    X, Y, tau, fast)
+    X, Y, eps, fast)
+

   grid <- array(0, dim = c(p, m, k))
   for (j in 1:p)
   {
   grid <- array(0, dim = c(p, m, k))
   for (j in 1:p)
   {

diff --git a/pkg/R/main.R b/pkg/R/main.R
index d29fe69..387d553 100644 (file)
--- a/pkg/R/main.R
+++ b/pkg/R/main.R
@@ -123,8 +123,7 @@ valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mi
       complexity = sumPen, contrast = -LLH)
   }))
   tableauRecap <- tableauRecap[which(tableauRecap[, 4] != Inf), ]
       complexity = sumPen, contrast = -LLH)
   }))
   tableauRecap <- tableauRecap[which(tableauRecap[, 4] != Inf), ]

-  
-  
+

   if (verbose == TRUE)
     print(tableauRecap)
   modSel <- capushe::capushe(tableauRecap, n)
   if (verbose == TRUE)
     print(tableauRecap)
   modSel <- capushe::capushe(tableauRecap, n)


@@ -141,11 +140,11 @@ valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mi
   {
     modSel@AIC_capushe$model
   }
   {
     modSel@AIC_capushe$model
   }

-    
+

   listMod <- as.integer(unlist(strsplit(as.character(indModSel), "[.]")))
   modelSel <- models_list[[listMod[1]]][[listMod[2]]]
   modelSel$tableau <- tableauRecap
   listMod <- as.integer(unlist(strsplit(as.character(indModSel), "[.]")))
   modelSel <- models_list[[listMod[1]]][[listMod[2]]]
   modelSel$tableau <- tableauRecap

-  
+

   if (plot)
     print(plot_valse(X, Y, modelSel, n))
 
   if (plot)
     print(plot_valse(X, Y, modelSel, n))
 

diff --git a/pkg/src/adapters/a.EMGLLF.c b/pkg/src/adapters/a.EMGLLF.c
index e93b7ec..9b004c2 100644 (file)
--- a/pkg/src/adapters/a.EMGLLF.c
+++ b/pkg/src/adapters/a.EMGLLF.c
@@ -46,7 +46,7 @@ SEXP EMGLLF(
        // OUTPUTS //
        /////////////
 
        // OUTPUTS //
        /////////////
 

-       SEXP phi, rho, pi, LLF, S, affec, dimPhiS, dimRho;
+       SEXP phi, rho, pi, llh, S, affec, dimPhiS, dimRho;

        PROTECT(dimPhiS = allocVector(INTSXP, 3));
        int* pDimPhiS = INTEGER(dimPhiS);
        pDimPhiS[0] = p; pDimPhiS[1] = m; pDimPhiS[2] = k;
        PROTECT(dimPhiS = allocVector(INTSXP, 3));
        int* pDimPhiS = INTEGER(dimPhiS);
        pDimPhiS[0] = p; pDimPhiS[1] = m; pDimPhiS[2] = k;


@@ -56,10 +56,10 @@ SEXP EMGLLF(
        PROTECT(phi = allocArray(REALSXP, dimPhiS));
        PROTECT(rho = allocArray(REALSXP, dimRho));
        PROTECT(pi = allocVector(REALSXP, k));
        PROTECT(phi = allocArray(REALSXP, dimPhiS));
        PROTECT(rho = allocArray(REALSXP, dimRho));
        PROTECT(pi = allocVector(REALSXP, k));

-       PROTECT(LLF = allocVector(REALSXP, maxi));
+       PROTECT(llh = allocVector(REALSXP, 1));

        PROTECT(S = allocArray(REALSXP, dimPhiS));
        PROTECT(affec = allocVector(INTSXP, n));
        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);
+       double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pLlh=REAL(llh), *pS=REAL(S);

        int *pAffec=INTEGER(affec);
 
        ////////////////////
        int *pAffec=INTEGER(affec);
 
        ////////////////////


@@ -67,14 +67,14 @@ SEXP EMGLLF(
        ////////////////////
 
        EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,eps,
        ////////////////////
 
        EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,eps,

-               pPhi,pRho,pPi,pLLF,pS,pAffec,
+               pPhi,pRho,pPi,pLlh,pS,pAffec,

                n,p,m,k);
 
        // Build list from OUT params and return it
        SEXP listParams, listNames;
        int nouts = 6;
        PROTECT(listParams = allocVector(VECSXP, nouts));
                n,p,m,k);
 
        // Build list from OUT params and return it
        SEXP listParams, listNames;
        int nouts = 6;
        PROTECT(listParams = allocVector(VECSXP, nouts));

-       char* lnames[6] = {"phi", "rho", "pi", "LLF", "S", "affec"}; //lists labels
+       char* lnames[6] = {"phi", "rho", "pi", "llh", "S", "affec"}; //lists labels

        PROTECT(listNames = allocVector(STRSXP,nouts));
        for (int i=0; i<nouts; i++)
                SET_STRING_ELT(listNames,i,mkChar(lnames[i]));
        PROTECT(listNames = allocVector(STRSXP,nouts));
        for (int i=0; i<nouts; i++)
                SET_STRING_ELT(listNames,i,mkChar(lnames[i]));


@@ -82,7 +82,7 @@ SEXP EMGLLF(
        SET_VECTOR_ELT(listParams, 0, phi);
        SET_VECTOR_ELT(listParams, 1, rho);
        SET_VECTOR_ELT(listParams, 2, pi);
        SET_VECTOR_ELT(listParams, 0, phi);
        SET_VECTOR_ELT(listParams, 1, rho);
        SET_VECTOR_ELT(listParams, 2, pi);

-       SET_VECTOR_ELT(listParams, 3, LLF);
+       SET_VECTOR_ELT(listParams, 3, llh);

        SET_VECTOR_ELT(listParams, 4, S);
        SET_VECTOR_ELT(listParams, 5, affec);
 
        SET_VECTOR_ELT(listParams, 4, S);
        SET_VECTOR_ELT(listParams, 5, affec);
 

diff --git a/test/generateRunSaveTest_EMGLLF.R b/test/generateRunSaveTest_EMGLLF.R
index 8b3f4c1..674a726 100644 (file)
--- a/test/generateRunSaveTest_EMGLLF.R
+++ b/test/generateRunSaveTest_EMGLLF.R
@@ -2,7 +2,7 @@ source("helper.R")
 library(valse)
 
 generateRunSaveTest_EMGLLF = function(n=200, p=15, m=10, k=3, mini=5, maxi=10,
 library(valse)
 
 generateRunSaveTest_EMGLLF = function(n=200, p=15, m=10, k=3, mini=5, maxi=10,

-       gamma=1., lambda=0.5, tau=1e-6)
+       gamma=1., lambda=0.5, eps=1e-6)

 {
        testFolder = "data/"
        dir.create(testFolder, showWarnings=FALSE, mode="0755")
 {
        testFolder = "data/"
        dir.create(testFolder, showWarnings=FALSE, mode="0755")


@@ -31,13 +31,13 @@ generateRunSaveTest_EMGLLF = function(n=200, p=15, m=10, k=3, mini=5, maxi=10,
                row.names=F, col.names=F)
        write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
                row.names=F, col.names=F)
                row.names=F, col.names=F)
        write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
                row.names=F, col.names=F)

-       write.table(as.double(tau), paste(testFolder,"tau",sep=""),
+       write.table(as.double(eps), paste(testFolder,"eps",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 = valse::EMGLLF(params$phiInit,params$rhoInit,params$piInit,params$gamInit,mini,
                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 = valse::EMGLLF(params$phiInit,params$rhoInit,params$piInit,params$gamInit,mini,

-               maxi,gamma,lambda,xy$X,xy$Y,tau,fast=FALSE)
+               maxi,gamma,lambda,xy$X,xy$Y,eps,fast=FALSE)

 
        #save outputs
        write.table(as.double(res$phi),paste(testFolder,"phi",sep=""),row.names=F,col.names=F)
 
        #save outputs
        write.table(as.double(res$phi),paste(testFolder,"phi",sep=""),row.names=F,col.names=F)

diff --git a/test/generateRunSaveTest_EMGrank.R b/test/generateRunSaveTest_EMGrank.R
index 935eada..9969620 100644 (file)
--- a/test/generateRunSaveTest_EMGrank.R
+++ b/test/generateRunSaveTest_EMGrank.R
@@ -4,7 +4,7 @@ library(valse)
 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))
 {

-  tau = 1e-6
+  eps = 1e-6

   pi = rep(1.0/k, k)
   rho = array(dim=c(m,m,k))
   for(i in 1:k)
   pi = rep(1.0/k, k)
   rho = array(dim=c(m,m,k))
   for(i in 1:k)


@@ -26,14 +26,14 @@ generateRunSaveTest_EMGrank = function(n=200, p=15, m=10, k=3, mini=5, maxi=10,
                row.names=F, col.names=F)
   write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
                row.names=F, col.names=F)
                row.names=F, col.names=F)
   write.table(as.double(xy$Y), paste(testFolder,"Y",sep=""),
                row.names=F, col.names=F)

-  write.table(as.double(tau), paste(testFolder,"tau",sep=""),
+  write.table(as.double(eps), paste(testFolder,"eps",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)
 
                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 = valse::EMGrank(pi,rho,mini,maxi,xy$X,xy$Y,tau,rank,fast=FALSE)
+  res = valse::EMGrank(pi,rho,mini,maxi,xy$X,xy$Y,eps,rank,fast=FALSE)

 
   #save output
   write.table(as.double(res$phi),paste(testFolder,"phi",sep=""),row.names=F,col.names=F)
 
   #save output
   write.table(as.double(res$phi),paste(testFolder,"phi",sep=""),row.names=F,col.names=F)

diff --git a/test/test.EMGLLF.c b/test/test.EMGLLF.c
index fa6e36c..68f73d9 100644 (file)
--- a/test/test.EMGLLF.c
+++ b/test/test.EMGLLF.c
@@ -23,7 +23,7 @@ int main(int argc, char** argv)
        Real lambda = read_real("lambda");
        Real* X = readArray_real("X");
        Real* Y = readArray_real("Y");
        Real lambda = read_real("lambda");
        Real* X = readArray_real("X");
        Real* Y = readArray_real("Y");

-       Real tau = read_real("tau");
+       Real eps = read_real("eps");

        ////////////
 
        /////////////
        ////////////
 
        /////////////


@@ -38,7 +38,7 @@ int main(int argc, char** argv)
 
        ////////////////////
        // Call to EMGLLF //
 
        ////////////////////
        // Call to EMGLLF //

-       EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau,
+       EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,eps,

                phi,rho,pi,&llh,S,affec,
                n,p,m,k);
        ////////////////////
                phi,rho,pi,&llh,S,affec,
                n,p,m,k);
        ////////////////////

diff --git a/test/test_EMGLLF.R b/test/test_EMGLLF.R
new file mode 100644 (file)
index 0000000..8ca7dea
--- /dev/null
+++ b/test/test_EMGLLF.R
@@ -0,0 +1,39 @@
+library(valse)
+testFolder = "data/"
+
+# NOTE: R typing is really terrible. as.double as.matrix ...and so on; don't remove.
+
+#get inputs
+npmk = as.matrix(read.table(paste(testFolder,"dimensions",sep="")))
+n = npmk[1]; p=npmk[2]; m=npmk[3]; k=npmk[4]
+phiInit = array(as.double(as.matrix(read.table(paste(testFolder,"phiInit",sep="")))), dim=c(p,m,k))
+rhoInit = array(as.double(as.matrix(read.table(paste(testFolder,"rhoInit",sep="")))), dim=c(m,m,k))
+piInit = as.double(as.matrix(read.table(paste(testFolder,"piInit",sep="")))[,])
+gamInit = matrix(as.double(as.matrix(read.table(paste(testFolder,"gamInit",sep="")))), n,k)
+mini = as.integer(as.matrix(read.table(paste(testFolder,"mini",sep="")))[1])
+maxi = as.integer(as.matrix(read.table(paste(testFolder,"maxi",sep="")))[1])
+gamma = as.double(as.matrix(read.table(paste(testFolder,"gamma",sep="")))[1])
+lambda = as.double(as.matrix(read.table(paste(testFolder,"lambda",sep="")))[1])
+X = matrix(as.double(as.matrix(read.table(paste(testFolder,"X",sep="")))), n,p)
+Y = matrix(as.double(as.matrix(read.table(paste(testFolder,"Y",sep="")))), n,m)
+eps = as.double(as.matrix(read.table(paste(testFolder,"eps",sep="")))[1])
+
+#get outputs
+phi = array(as.double(as.matrix(read.table(paste(testFolder,"phi",sep="")))), dim=c(p,m,k))
+rho = array(as.double(as.matrix(read.table(paste(testFolder,"rho",sep="")))), dim=c(m,m,k))
+pi = as.double(as.matrix(read.table(paste(testFolder,"pi",sep="")))[,])
+llh = as.double(as.matrix(read.table(paste(testFolder,"llh",sep="")))[1])
+S = array(as.double(as.matrix(read.table(paste(testFolder,"S",sep="")))), dim=c(p,m,k))
+affec = as.double(as.matrix(read.table(paste(testFolder,"affec",sep="")))[,])
+
+res = valse::EMGLLF(
+       phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,eps,fast=TRUE)
+
+#compare outputs
+nd=7 #number of digits
+print( all(round(phi,nd) == round(res$phi,nd)) )
+print( all(round(rho,nd) == round(res$rho,nd)) )
+print( all(round(pi,nd) == round(res$pi,nd)) )
+print( all(round(llh,nd) == round(res$llh,nd)) )
+print( all(round(S,nd) == round(res$S,nd)) )
+print( all(affec == res$affec) )