From: Benjamin Auder <benjamin.auder@somewhere> Date: Fri, 3 Nov 2017 09:54:04 +0000 (+0100) Subject: progress in debug: fix LLF/llh mismatch, and length + add adapter test in test/ X-Git-Url: https://git.auder.net/%7B%7B%20asset%28%27mixstore/css/user/doc/pieces/%3C?a=commitdiff_plain;h=6279ba8656582370e7242ff9e77d22c23fe8ca04;p=valse.git 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 --- 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), - 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") } diff --git a/pkg/R/EMGrank.R b/pkg/R/EMGrank.R index b85a0fa..f2bf58e 100644 --- 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 tau Seuil pour accepter la convergence +#' @param eps Seuil pour accepter la convergence #' @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 -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 - 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 @@ -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 - .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") } @@ -43,7 +43,7 @@ matricize <- function(X) } # 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) @@ -64,7 +64,7 @@ matricize <- function(X) # 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) diff --git a/pkg/R/computeGridLambda.R b/pkg/R/computeGridLambda.R index 8449d10..ac0788a 100644 --- 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 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, - maxi, tau, fast) + maxi, eps, fast) { 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, - X, Y, tau, fast) + X, Y, eps, fast) + 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 --- 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), ] - - + 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 } - + 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)) diff --git a/pkg/src/adapters/a.EMGLLF.c b/pkg/src/adapters/a.EMGLLF.c index e93b7ec..9b004c2 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, 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; @@ -56,10 +56,10 @@ SEXP EMGLLF( 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)); - 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); //////////////////// @@ -67,14 +67,14 @@ SEXP EMGLLF( //////////////////// 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)); - 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])); @@ -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, 3, LLF); + SET_VECTOR_ELT(listParams, 3, llh); 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 --- 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, - 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") @@ -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) - 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, - 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) diff --git a/test/generateRunSaveTest_EMGrank.R b/test/generateRunSaveTest_EMGrank.R index 935eada..9969620 100644 --- 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)) { - tau = 1e-6 + eps = 1e-6 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) - 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) - 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) diff --git a/test/test.EMGLLF.c b/test/test.EMGLLF.c index fa6e36c..68f73d9 100644 --- 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 tau = read_real("tau"); + Real eps = read_real("eps"); //////////// ///////////// @@ -38,7 +38,7 @@ int main(int argc, char** argv) //////////////////// // 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); //////////////////// diff --git a/test/test_EMGLLF.R b/test/test_EMGLLF.R new file mode 100644 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) )