From c3bc47052f3ccb659659c59a82e9a99ea842398d Mon Sep 17 00:00:00 2001 From: Benjamin Auder 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= 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 -#include - 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 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 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 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= 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