From 09ab3c164abb566764e86a175b5973241e708fd6 Mon Sep 17 00:00:00 2001 From: Benjamin Auder Date: Sat, 17 Dec 2016 17:34:11 +0100 Subject: [PATCH] R package can now be installed (compilation OK) --- R/main.R | 65 ++++---- R/selectVariables.R | 147 +++++++++--------- src/adapters/a.EMGLLF.c | 10 +- src/adapters/a.EMGrank.c | 10 +- src/adapters/a.constructionModelesLassoMLE.c | 16 +- src/adapters/a.constructionModelesLassoRank.c | 6 +- src/adapters/a.selectiontotale.c | 37 +++-- src/sources/EMGLLF.c | 2 +- src/sources/EMGLLF.h | 2 +- src/sources/EMGrank.c | 2 +- src/sources/EMGrank.h | 2 +- src/sources/constructionModelesLassoMLE.c | 6 +- src/sources/constructionModelesLassoMLE.h | 2 +- src/sources/constructionModelesLassoRank.c | 6 +- src/sources/constructionModelesLassoRank.h | 2 +- src/sources/selectiontotale.c | 2 +- src/sources/selectiontotale.h | 2 +- 17 files changed, 164 insertions(+), 155 deletions(-) diff --git a/R/main.R b/R/main.R index b9b8d5b..059843f 100644 --- a/R/main.R +++ b/R/main.R @@ -14,13 +14,13 @@ Valse = setRefClass( # Optionally user defined (some default values) # power in the penalty - gamma = "double", + gamma = "numeric", # minimum number of iterations for EM algorithm mini = "integer", # maximum number of iterations for EM algorithm maxi = "integer", # threshold for stopping EM algorithm - eps = "double", + eps = "numeric", # minimum number of components in the mixture kmin = "integer", # maximum number of components in the mixture @@ -31,28 +31,28 @@ Valse = setRefClass( # Computed through the workflow # initialisation for the reparametrized conditional mean parameter - phiInit, + phiInit = "numeric", # initialisation for the reparametrized variance parameter - rhoInit, + rhoInit = "numeric", # initialisation for the proportions - piInit, + piInit = "numeric", # initialisation for the allocations probabilities in each component - tauInit, + tauInit = "numeric", # values for the regularization parameter grid - gridLambda = c(), + gridLambda = "numeric", # je ne crois pas vraiment qu'il faille les mettre en sortie, d'autant plus qu'on construit # une matrice A1 et A2 pour chaque k, et elles sont grandes, donc ca coute un peu cher ... - A1, - A2, + A1 = "integer", + A2 = "integer", # collection of estimations for the reparametrized conditional mean parameters - Phi, + Phi = "numeric", # collection of estimations for the reparametrized variance parameters - Rho, + Rho = "numeric", # collection of estimations for the proportions parameters - Pi, + Pi = "numeric", - #immutable - seuil = 1e-15 + #immutable (TODO:?) + seuil = "numeric" ), methods = list( @@ -75,6 +75,7 @@ Valse = setRefClass( kmax <<- ifelse (hasArg("kmax"), kmax, as.integer(3)) rangmin <<- ifelse (hasArg("rangmin"), rangmin, as.integer(2)) rangmax <<- ifelse (hasArg("rangmax"), rangmax, as.integer(3)) + seuil <<- 1e-15 #immutable (TODO:?) }, ################################## @@ -160,22 +161,22 @@ Valse = setRefClass( Pi2 = Pi p = ncol(X) m = ncol(Y) - if size(Phi2) == 0 + if (is.null(dim(Phi2))) #test was: size(Phi2) == 0 { - Phi[,,1:k] = r1$phi - Rho[,,1:k] = r1$rho - Pi[1:k,] = r1$pi + Phi[,,1:k] <<- r1$phi + Rho[,,1:k] <<- r1$rho + Pi[1:k,] <<- r1$pi } else { - Phi = array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(r1$phi)[4])) - Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] = Phi2 - Phi[,,1:k,dim(Phi2)[4]+1] = r1$phi - Rho = array(0., dim=c(m,m,kmax,dim(Rho2)[4]+dim(r1$rho)[4])) - Rho[,,1:(dim(Rho2)[3]),1:(dim(Rho2)[4])] = Rho2 - Rho[,,1:k,dim(Rho2)[4]+1] = r1$rho - Pi = array(0., dim=c(kmax,dim(Pi2)[2]+dim(r1$pi)[2])) - Pi[1:nrow(Pi2),1:ncol(Pi2)] = Pi2 - Pi[1:k,ncol(Pi2)+1] = r1$pi + Phi <<- array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(r1$phi)[4])) + Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] <<- Phi2 + Phi[,,1:k,dim(Phi2)[4]+1] <<- r1$phi + Rho <<- array(0., dim=c(m,m,kmax,dim(Rho2)[4]+dim(r1$rho)[4])) + Rho[,,1:(dim(Rho2)[3]),1:(dim(Rho2)[4])] <<- Rho2 + Rho[,,1:k,dim(Rho2)[4]+1] <<- r1$rho + Pi <<- array(0., dim=c(kmax,dim(Pi2)[2]+dim(r1$pi)[2])) + Pi[1:nrow(Pi2),1:ncol(Pi2)] <<- Pi2 + Pi[1:k,ncol(Pi2)+1] <<- r1$pi } } else { @@ -183,14 +184,12 @@ Valse = setRefClass( Phi2 = Phi if (dim(Phi2)[1] == 0) { - Phi(:,:,1:k,:) = phi + Phi[,,1:k,] <<- phi } else { - size(Phi2) - Phi = zeros(p,m,kmax,size(Phi2,4)+size(phi,4)) - size(Phi) - Phi(:,:,1:size(Phi2,3),1:size(Phi2,4)) = Phi2 - Phi(:,:,1:k,size(Phi2,4)+1:end) = phi + Phi <<- array(0., dim=c(p,m,kmax,dim(Phi2)[4]+dim(phi)[4])) + Phi[,,1:(dim(Phi2)[3]),1:(dim(Phi2)[4])] <<- Phi2 + Phi[,,1:k,-(1:(dim(Phi2)[4]))] <<- phi } } } diff --git a/R/selectVariables.R b/R/selectVariables.R index be53d85..e1a4e33 100644 --- a/R/selectVariables.R +++ b/R/selectVariables.R @@ -1,82 +1,89 @@ -#' selectVaribles -#' It is a function which construct, for a given lambda, the sets of +#' selectVaribles +#' It is a function which construct, for a given lambda, the sets of #' relevant variables and irrelevant variables. #' #' @param phiInit an initial estimator for phi (size: p*m*k) #' @param rhoInit an initial estimator for rho (size: m*m*k) -#' @param piInit an initial estimator for pi (size : k) +#' @param piInit an initial estimator for pi (size : k) #' @param gamInit an initial estimator for gamma -#' @param mini minimum number of iterations in EM algorithm -#' @param maxi maximum number of iterations in EM algorithm -#' @param gamma power in the penalty +#' @param mini minimum number of iterations in EM algorithm +#' @param maxi maximum number of iterations in EM algorithm +#' @param gamma power in the penalty #' @param glambda grid of regularization parameters -#' @param X matrix of regressors -#' @param Y matrix of responses -#' @param thres threshold to consider a coefficient to be equal to 0 -#' @param tau threshold to say that EM algorithm has converged +#' @param X matrix of regressors +#' @param Y matrix of responses +#' @param thres threshold to consider a coefficient to be equal to 0 +#' @param tau threshold to say that EM algorithm has converged #' #' @return #' @export #' #' @examples selectVariables <- function(phiInit,rhoInit,piInit,gamInit, - mini,maxi,gamma,glambda,X,Y,thres,tau){ - - dimphi <- dim(phiInit) - p <- dimPhi[1] - m <- dimPhi[2] - k <- dimPhi[3] - L <- length(glambda); - A1 <- array(0, dim <- c(p,m+1,L)) - A2 <- array(0, dim <- c(p,m+1,L)) - Rho <- array(0, dim <- c(m,m,k,L)) - Pi <- array(0, dim <- c(k,L)); - - # For every lambda in gridLambda, comutation of the coefficients - for (lambdaIndex in c(1:L)) { - Res <- EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi, - gamma,glambda[lambdaIndex],X,Y,tau); - phi <- Res$phi - rho <- Res$rho - pi <- Res$pi - - # If a coefficient is larger than the threshold, we keep it - selectedVariables <- array(0, dim = c(p,m)) - discardedVariables <- array(0, dim = c(p,m)) - atLeastOneSelectedVariable <- false - for (j in c(1:p)){ - cpt <- 1 - cpt2 <-1 - for (mm in c(1:m)){ - if (max(abs(phi[j,mm,])) > thres){ - selectedVariables[j,cpt] <- mm - cpt <- cpt+1 - atLeastOneSelectedVariable <- true - } else{ - discardedVariables[j,cpt2] <- mm - cpt2 <- cpt2+1 - } - } - } - - # If no coefficients have been selected, we provide the zero matrix - # We delete zero coefficients: vec = indices of zero values - if atLeastOneSelectedVariable{ - vec <- c() - for (j in c(1:p)){ - if (selectedVariables(j,1) =! 0){ - vec <- c(vec,j) - } - } - # Else, we provide the indices of relevant coefficients - A1[,1,lambdaIndex] <- c(vec,rep(0,p-length(vec))) - A1[1:length(vec),2:(m+1),lambdaIndex] <- selectedVariables[vec,] - A2[,1,lambdaIndex] <- 1:p - A2[,2:(m+1),lambdaIndex] <- discardedVariables - Rho[,,,lambdaIndex] <- rho - Pi[,lambdaIndex] <- pi - } - - } - return(res = list(A1 = A1, A2 = A2 , Rho = Rho, Pi = Pi)) -} \ No newline at end of file + mini,maxi,gamma,glambda,X,Y,thres,tau) +{ + dimphi <- dim(phiInit) + p <- dimPhi[1] + m <- dimPhi[2] + k <- dimPhi[3] + L <- length(glambda); + A1 <- array(0, dim <- c(p,m+1,L)) + A2 <- array(0, dim <- c(p,m+1,L)) + Rho <- array(0, dim <- c(m,m,k,L)) + Pi <- array(0, dim <- c(k,L)); + + # For every lambda in gridLambda, comutation of the coefficients + for (lambdaIndex in c(1:L)) + { + Res <- EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi, + gamma,glambda[lambdaIndex],X,Y,tau); + phi <- Res$phi + rho <- Res$rho + pi <- Res$pi + + # If a coefficient is larger than the threshold, we keep it + selectedVariables <- array(0, dim = c(p,m)) + discardedVariables <- array(0, dim = c(p,m)) + atLeastOneSelectedVariable <- false + for (j in c(1:p)) + { + cpt <- 1 + cpt2 <-1 + for (mm in c(1:m)) + { + if (max(abs(phi[j,mm,])) > thres) + { + selectedVariables[j,cpt] <- mm + cpt <- cpt+1 + atLeastOneSelectedVariable <- true + } else + { + discardedVariables[j,cpt2] <- mm + cpt2 <- cpt2+1 + } + } + } + + # If no coefficients have been selected, we provide the zero matrix + # We delete zero coefficients: vec = indices of zero values + if (atLeastOneSelectedVariable) + { + vec <- c() + for (j in c(1:p)) + { + if (selectedVariables(j,1) != 0) + vec <- c(vec,j) + # Else ( NOTE: [auder] else ?! TODO: explain? ) + # we provide the indices of relevant coefficients + A1[,1,lambdaIndex] <- c(vec,rep(0,p-length(vec))) + A1[1:length(vec),2:(m+1),lambdaIndex] <- selectedVariables[vec,] + A2[,1,lambdaIndex] <- 1:p + A2[,2:(m+1),lambdaIndex] <- discardedVariables + Rho[,,,lambdaIndex] <- rho + Pi[,lambdaIndex] <- pi + } + } + } + + return(res = list(A1 = A1, A2 = A2 , Rho = Rho, Pi = Pi)) +} diff --git a/src/adapters/a.EMGLLF.c b/src/adapters/a.EMGLLF.c index f70d69d..0df00bd 100644 --- a/src/adapters/a.EMGLLF.c +++ b/src/adapters/a.EMGLLF.c @@ -63,7 +63,7 @@ SEXP EMGLLF( // Call to EMGLLF // //////////////////// - EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau, + EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau, pPhi,pRho,pPi,pLLF,pS, n,p,m,k); @@ -75,11 +75,11 @@ SEXP EMGLLF( for (int i=0; i<5; i++) SET_STRING_ELT(listNames,i,mkChar(lnames[i])); setAttrib(listParams, R_NamesSymbol, listNames); - SET_ARRAY_ELT(listParams, 0, phi); - SET_ARRAY_ELT(listParams, 1, rho); - SET_MATRIX_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_ARRAY_ELT(listParams, 4, S); + SET_VECTOR_ELT(listParams, 4, S); UNPROTECT(9); return listParams; diff --git a/src/adapters/a.EMGrank.c b/src/adapters/a.EMGrank.c index 763ff4e..469349c 100644 --- a/src/adapters/a.EMGrank.c +++ b/src/adapters/a.EMGrank.c @@ -2,7 +2,7 @@ #include #include "EMGrank.h" -SEXP EMGLLF( +SEXP EMGrank( SEXP Pi_, SEXP Rho_, SEXP mini_, @@ -34,13 +34,13 @@ SEXP EMGLLF( double* Rho = REAL(Rho_); double* X = REAL(X_); double* Y = REAL(Y_); - double* rank = REAL(rank_); + int* rank = INTEGER(rank_); ///////////// // OUTPUTS // ///////////// - SEXP phi, LLF; + SEXP phi, LLF, dimPhi; PROTECT(dimPhi = allocVector(INTSXP, 3)); int* pDimPhi = INTEGER(dimPhi); pDimPhi[0] = p; pDimPhi[1] = m; pDimPhi[2] = k; @@ -52,7 +52,7 @@ SEXP EMGLLF( // Call to EMGrank // ///////////////////// - EMGrank(Pi, Rho, mini, maxi, X, Y, tau, rank, + EMGrank_core(Pi, Rho, mini, maxi, X, Y, tau, rank, pPhi,pLLF, n,p,m,k); @@ -64,7 +64,7 @@ SEXP EMGLLF( for (int i=0; i<2; i++) SET_STRING_ELT(listNames,i,mkChar(lnames[i])); setAttrib(listParams, R_NamesSymbol, listNames); - SET_ARRAY_ELT(listParams, 0, phi); + SET_VECTOR_ELT(listParams, 0, phi); SET_VECTOR_ELT(listParams, 1, LLF); UNPROTECT(5); diff --git a/src/adapters/a.constructionModelesLassoMLE.c b/src/adapters/a.constructionModelesLassoMLE.c index 8543658..72c3173 100644 --- a/src/adapters/a.constructionModelesLassoMLE.c +++ b/src/adapters/a.constructionModelesLassoMLE.c @@ -2,7 +2,7 @@ #include #include "constructionModelesLassoMLE.h" -SEXP EMGLLF( +SEXP constructionModelesLassoMLE( SEXP phiInit_, SEXP rhoInit_, SEXP piInit_, @@ -24,7 +24,7 @@ SEXP EMGLLF( int p = INTEGER(dim)[0]; int m = INTEGER(dim)[1]; int k = INTEGER(dim)[2]; - int L = INTEGER(getAttrib(glambda_, R_LengthSymbol))[0]; + int L = length(glambda_); //////////// // INPUTS // @@ -45,8 +45,8 @@ SEXP EMGLLF( double* glambda = REAL(glambda_); double* X = REAL(X_); double* Y = REAL(Y_); - double* A1 = REAL(A1_); - double* A2 = REAL(A2_); + int* A1 = INTEGER(A1_); + int* A2 = INTEGER(A2_); ///////////// // OUTPUTS // @@ -69,7 +69,7 @@ SEXP EMGLLF( // Call to constructionModelesLassoMLE // ///////////////////////////////////////// - constructionModelesLassoMLE( + constructionModelesLassoMLE_core( phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2, pPhi,pRho,pPi,pLvraisemblance, n,p,m,k,L); @@ -82,9 +82,9 @@ SEXP EMGLLF( for (int i=0; i<4; i++) SET_STRING_ELT(listNames,i,mkChar(lnames[i])); setAttrib(listParams, R_NamesSymbol, listNames); - SET_ARRAY_ELT(listParams, 0, phi); - SET_ARRAY_ELT(listParams, 1, rho); - SET_MATRIX_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, lvraisemblance); UNPROTECT(8); diff --git a/src/adapters/a.constructionModelesLassoRank.c b/src/adapters/a.constructionModelesLassoRank.c index 0d056c8..4833e19 100644 --- a/src/adapters/a.constructionModelesLassoRank.c +++ b/src/adapters/a.constructionModelesLassoRank.c @@ -2,7 +2,7 @@ #include #include "constructionModelesLassoRank.h" -SEXP EMGLLF( +SEXP constructionModelesLassoRank( SEXP Pi_, SEXP Rho_, SEXP mini_, @@ -58,7 +58,7 @@ SEXP EMGLLF( // Call to constructionModelesLassoRank // ////////////////////////////////////////// - constructionModelesLassoRank( + constructionModelesLassoRank_core( Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax, pPhi,pLvraisemblance, n,p,m,k,L); @@ -71,7 +71,7 @@ SEXP EMGLLF( for (int i=0; i<2; i++) SET_STRING_ELT(listNames,i,mkChar(lnames[i])); setAttrib(listParams, R_NamesSymbol, listNames); - SET_ARRAY_ELT(listParams, 0, phi); + SET_VECTOR_ELT(listParams, 0, phi); SET_VECTOR_ELT(listParams, 1, lvraisemblance); UNPROTECT(5); diff --git a/src/adapters/a.selectiontotale.c b/src/adapters/a.selectiontotale.c index 3d647fb..3bfab9f 100644 --- a/src/adapters/a.selectiontotale.c +++ b/src/adapters/a.selectiontotale.c @@ -2,7 +2,7 @@ #include #include "selectiontotale.h" -SEXP EMGLLF( +SEXP selectiontotale( SEXP phiInit_, SEXP rhoInit_, SEXP piInit_, @@ -20,10 +20,10 @@ SEXP EMGLLF( SEXP dimX = getAttrib(X_, R_DimSymbol); int n = INTEGER(dimX)[0]; int p = INTEGER(dimX)[1]; - SEXP dimRho = getAttrib(rhoInit_, R_DimSymbol); - int m = INTEGER(dimRho)[0]; - int k = INTEGER(dimRho)[2]; - int L = INTEGER(getAttrib(glambda_, R_LengthSymbol))[0]; + SEXP dimRhoInit = getAttrib(rhoInit_, R_DimSymbol); + int m = INTEGER(dimRhoInit)[0]; + int k = INTEGER(dimRhoInit)[2]; + int L = length(glambda_); //////////// // INPUTS // @@ -37,7 +37,7 @@ SEXP EMGLLF( double tau = NUMERIC_VALUE(tau_); // Get pointers from SEXP arrays ; WARNING: by columns ! - double* piInit = REAL(phiInit_); + double* phiInit = REAL(phiInit_); double* rhoInit = REAL(rhoInit_); double* piInit = REAL(piInit_); double* gamInit = REAL(gamInit_); @@ -49,22 +49,25 @@ SEXP EMGLLF( // OUTPUTS // ///////////// - int Size = pow(rangmax-rangmin+1,k); - SEXP A1, A2, rho, pi, dimA; + SEXP A1, A2, rho, pi, dimA, dimRho; PROTECT(dimA = allocVector(INTSXP, 3)); int* pDimA = INTEGER(dimA); pDimA[0] = p; pDimA[1] = m+1; pDimA[2] = L; - PROTECT(A1 = allocArray(REALSXP, dimA)); - PROTECT(A2 = allocArray(REALSXP, dimA)); - PROTECT(rho = allocArray(REALSXP, dimRho); + PROTECT(A1 = allocArray(INTSXP, dimA)); + PROTECT(A2 = allocArray(INTSXP, dimA)); + PROTECT(dimRho = allocVector(INTSXP, 4)); + int* pDimRho = INTEGER(dimRho); + pDimRho[0] = m; pDimRho[1] = m; pDimRho[2] = k; pDimRho[3] = L; + PROTECT(rho = allocArray(REALSXP, dimRho)); PROTECT(pi = allocMatrix(REALSXP, k, L)); - double *pA1=REAL(A1), *pA2=REAL(A2), *pRho=REAL(rho), *pPi=REAL(pi); + int *pA1=INTEGER(A1), *pA2=INTEGER(A2); + double *pRho=REAL(rho), *pPi=REAL(pi); ///////////////////////////// // Call to selectiontotale // ///////////////////////////// - selectiontotale( + selectiontotale_core( phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau, pA1,pA2,pRho,pPi, n,p,m,k,L); @@ -77,10 +80,10 @@ SEXP EMGLLF( for (int i=0; i<4; i++) SET_STRING_ELT(listNames,i,mkChar(lnames[i])); setAttrib(listParams, R_NamesSymbol, listNames); - SET_ARRAY_ELT(listParams, 0, A1); - SET_ARRAY_ELT(listParams, 1, A2); - SET_ARRAY_ELT(listParams, 2, rho); - SET_MATRIX_ELT(listParams, 3, pi); + SET_VECTOR_ELT(listParams, 0, A1); + SET_VECTOR_ELT(listParams, 1, A2); + SET_VECTOR_ELT(listParams, 2, rho); + SET_VECTOR_ELT(listParams, 3, pi); UNPROTECT(7); return listParams; diff --git a/src/sources/EMGLLF.c b/src/sources/EMGLLF.c index 6966e9c..b305641 100644 --- a/src/sources/EMGLLF.c +++ b/src/sources/EMGLLF.c @@ -3,7 +3,7 @@ #include // TODO: don't recompute indexes every time...... -void EMGLLF( +void EMGLLF_core( // IN parameters const double* phiInit, // parametre initial de moyenne renormalisé const double* rhoInit, // parametre initial de variance renormalisé diff --git a/src/sources/EMGLLF.h b/src/sources/EMGLLF.h index df4ae7c..d09d27c 100644 --- a/src/sources/EMGLLF.h +++ b/src/sources/EMGLLF.h @@ -1,7 +1,7 @@ #ifndef valse_EMGLLF_H #define valse_EMGLLF_H -void EMGLLF( +void EMGLLF_core( // IN parameters const double* phiInit, const double* rhoInit, diff --git a/src/sources/EMGrank.c b/src/sources/EMGrank.c index 5ee44e3..0791892 100644 --- a/src/sources/EMGrank.c +++ b/src/sources/EMGrank.c @@ -38,7 +38,7 @@ static double* pinv(const double* matrix, int dim) } // TODO: comment EMGrank purpose -void EMGrank( +void EMGrank_core( // IN parameters const double* Pi, // parametre de proportion const double* Rho, // parametre initial de variance renormalisé diff --git a/src/sources/EMGrank.h b/src/sources/EMGrank.h index 53a1d01..8123e68 100644 --- a/src/sources/EMGrank.h +++ b/src/sources/EMGrank.h @@ -1,7 +1,7 @@ #ifndef valse_EMGrank_H #define valse_EMGrank_H -void EMGrank( +void EMGrank_core( // IN parameters const double* Pi, const double* Rho, diff --git a/src/sources/constructionModelesLassoMLE.c b/src/sources/constructionModelesLassoMLE.c index 272f01b..4ab62ad 100644 --- a/src/sources/constructionModelesLassoMLE.c +++ b/src/sources/constructionModelesLassoMLE.c @@ -5,7 +5,7 @@ #include // TODO: comment on constructionModelesLassoMLE purpose -void constructionModelesLassoMLE( +void constructionModelesLassoMLE_core( // IN parameters const double* phiInit, // parametre initial de moyenne renormalisé const double* rhoInit, // parametre initial de variance renormalisé @@ -83,7 +83,7 @@ void constructionModelesLassoMLE( double* piLambda = (double*)malloc(k*sizeof(double)); double* LLF = (double*)malloc((maxi+1)*sizeof(double)); double* S = (double*)malloc(lengthA*m*k*sizeof(double)); - EMGLLF(phia,rhoInit,piInit,gamInit,mini,maxi,gamma,0.0,Xa,Y,tau, + EMGLLF_core(phia,rhoInit,piInit,gamInit,mini,maxi,gamma,0.0,Xa,Y,tau, phiLambda,rhoLambda,piLambda,LLF,S, n,lengthA,m,k); free(Xa); @@ -138,7 +138,7 @@ void constructionModelesLassoMLE( for (int mm=0; mm